wintering_seabird_trend_2021.Rmd
This script, derived from the 2016/17 BTO scripts, incorporates the code updates required for 2021. The script
library(b1indicator)
library(reshape)
library(dplyr)
library(purrr)
library(stringr)
library(tidyr)
library(fs)
library(luzlogr)
Set the year range, add the regions and sub region labels into vectors and create an empty dat_sum list to store the species+regions/subregions yearly and six yearly mean relative abundances.
# Set year ranges
startYearWintering <- list(
In = 1991,
Io = 1991,
I = 1991,
IIa = 1991,
IIb = 1991,
IIc = 1991,
IId = 1991,
IIdw = 1991, # Missing Wadden sea counts have effect on species trend
IIe = 1991,
IIf = 1991,
II = 1991,
IIw = 1991, # Missing Wadden sea counts have effect on species trend
III = 1991
)
endYearWintering <- list(
In = 2020,
Io = 2020,
I = 2020,
IIa = 2020,
IIb = 2020,
IIc = 2020,
IId = 2020,
IIdw = 2016, # Missing Wadden sea counts have effect on species trend
IIe = 2020,
IIf = 2020,
II = 2020,
IIw = 2016, # Missing Wadden sea counts have effect on species trend
III = 2020
)
# Previous assessment year used to compare threshold category changes in traffic light plots
previousYearWintering <- 2014
# species where reduced Wadden year range 1991 - 2016 should be used instead of 1991 -2020
wadden_species <- c("Bar_tailed_godwit", "Barnacle_goose", "Black_headed_gull", "Brent_goose", "Common_gull", "Curlew_sandpiper", "Eurasian_curlew", "Eurasian_spoonbill", "European_golden_plover", "Great_black_backed_gull", "Great_cormorant", "Great_crested_grebe", "Common_greenshank", "Grey_plover", "European_herring_gull", "Kentish_plover", "Northern_lapwing", "Lesser_black_backed_gull", "Mallard", "Mute_swan", "Eurasian_oystercatcher", "Pied_avocet", "Northern_pintail", "Red_breasted_merganser", "Red_knot", "Common_redshank", "Common_ringed_plover", "Ruff", "Sanderling", "Common_shelduck", "Northern_shoveler", "Slavonian_grebe", "Smew", "Spotted_redshank", "Eurasian_teal", "Ruddy_turnstone", "Eurasian_whimbrel", "Whooper_swan", "Eurasian_wigeon")
# list of OSPAR regions and bird subdivision
region_labels <- c("I", "II", "III")
sub_region_labels <- c("In", "Io", "IIa", "IIb", "IIc", "IId", "IIe", "IIf")
region_sub_region_labels <- c(region_labels, sub_region_labels)
# empty list to store species relative abundance values for regions and sub regions
dat_unsm <-list()
Create the directory structure for the results
# directory to place wintering waterbird data
dir.create(path("results"))
dir.create(path("results", "non_breeding_out"))
dir.create(path("results", "non_breeding_baseline_model"))
# and directory for EMECO Validation datasets
dir.create(path("results", "EMECO checks"))
dir.create(path("results", "EMECO checks", "non_breeding_out"))
dir.create(path("results", "EMECO checks", "species data"))
Import the species and sites information
# get .csv file containing species names & number of eggs
add.info = read.table(path("data", "number of eggs.csv"), header = TRUE, sep = ",",
stringsAsFactors = FALSE)
# standardise common names with underscores and not spaces or dashes in name
add.info$Common_name = str_replace_all(add.info$Common_name, "[:space:]|-", "_")
# bring in sites list to get ospar regions/sub-regions
sites = read.csv(path("data", "BreedingSites ABUNDANCE.csv"), header = TRUE)
Import the species yearly relative abundance counts into a single data frame, data
#read data
data <- read.csv(path("data", "wintering", "WinteringWaterbirdAbundance.csv"),
stringsAsFactors = FALSE)
data_I <- data[data$OSPAR_region == "I" & data$Year %in% startYearWintering$I:endYearWintering$I, ]
data_II <- data[data$OSPAR_region == "II" & data$Year %in% startYearWintering$II:endYearWintering$II, ]
data_III <- data[data$OSPAR_region == "III" & data$Year %in% startYearWintering$III:endYearWintering$III, ]
data <- bind_rows(data_I, data_II, data_III)
# standardise common names with underscores and not spaces or dashes in name
data$CommonName <- str_replace_all(data$CommonName, "[:space:]|-", "_")
Create Region_Subdivision column containing unique Region_Subdivisions labels over which to sum data
data$Region_Subdivision <- do.call(paste,data[c("OSPAR_region","OSPAR_sub.division")])
data$Region_Subdivision <- do.call(str_replace_all, list(data$Region_Subdivision, " NA", ""))
data$Region_Subdivision <- do.call(str_replace_all, list(data$Region_Subdivision, " ", ""))
Sum the data for each species across Regional Subdivisions, into data frame data.agg
data.agg <- aggregate(data$Count, list(data$Region_Subdivision, data$CommonName,
data$SpeciesID, data$Year), FUN = "sum")
names(data.agg) <- c("Region_Subdivision","CommonName", "SpeciesID", "Year", "SumCount")
data.agg <- data.agg[order(data.agg$CommonName, data.agg$Region_Subdivision, data.agg$Year),]
Get list of species in vector species
species <- unique(data.agg$CommonName)
species <- as.character(species)
For each species
# log setting baseline regression results
luzlogr::openlog(path("results", "non_breeding_baseline_model", "Wintering Seabirds setting baseline.csv"))
luzlogr::printlog("species,region,count_flag,p-value,method,baseline_count,first_relative_abundance", ts = FALSE)
for(i in 1:length(species)){
# create species directory
dir.create(path("results", "non_breeding_out", species[i]))
dir.create(path("results", "non_breeding_baseline_model", species[i]))
# get additional info
species.info = subset(add.info, Common_name == species[i])
# set target 0.7 if species lays >1 egg, 0.8 if species lays 1 egg
ifelse(species.info$Eggs[1] == 1, target <- 0.8, target <- 0.7)
# get region and subregion abundance data for species
new_dat <- data.agg[data.agg$CommonName == species[i],]
new_dat <- new_dat[, c("Year", "Region_Subdivision", "SumCount")]
names(new_dat) <- c("Year", "Colony", "Count")
new_dat$Colony <- paste("OSPAR", new_dat$Colony, sep = "")
# convert abundance count data to wide format with each year a column and add
## site information columns
new_dat = cast(new_dat, Colony ~ Year, value.var = Count, fun.aggregate = mean)
new_dat = merge(new_dat, sites, all.x = TRUE)
# export species data
write.csv(new_dat, path("results", "EMECO checks", "species data",
str_glue("{species[i]}_non_breeding.csv")), row.names = FALSE)
## calculate and plot relative abundance values for regions
for(j in as.numeric(as.roman(region_labels))) {
# Update end year for each OSPAR REGION
if (j == 1) {
endYear = endYearWintering$I
startYear = startYearWintering$I
} else if (j == 2 & !species[i] %in% wadden_species) {
endYear = endYearWintering$II
startYear = startYearWintering$II
} else if (j == 2 & species[i] %in% wadden_species) {
endYear = endYearWintering$IIw
startYear = startYearWintering$IIw
} else if (j == 3) {
endYear = endYearWintering$III
startYear = startYearWintering$III
} else if (j == 4) {
endYear = endYearWintering$IV
startYear = startYearWintering$IV
} else if (j == 5) {
endYear = endYearWintering$V
startYear = startYearWintering$V
}
dat_unsm[[str_glue("{species[i]}_{as.roman(j)}")]] <-
calculate_non_breeding_relative_abundance(ospar = subset(new_dat, OSPAR_REGION == j),
label = as.roman(j),
species_name = species[i],
target_threshold = target,
start_year = startYear,
end_year = endYear)
}
# ## calculate and plot relative abundance values for sub regions
for(j in sub_region_labels) {
# Update start and end year for each OSPAR SUB REGION
if (j == "Im") {
endYear = endYearWintering$Im
startYear = startYearWintering$Im
} else if (j == "In") {
endYear = endYearWintering$In
startYear = startYearWintering$In
} else if (j == "Io") {
endYear = endYearWintering$Io
startYear = startYearWintering$Io
} else if (j == "IIa") {
endYear = endYearWintering$IIa
startYear = startYearWintering$IIa
} else if (j == "IIb") {
endYear = endYearWintering$IIb
startYear = startYearWintering$IIb
} else if (j == "IIc") {
endYear = endYearWintering$IIc
startYear = startYearWintering$IIc
} else if (j == "IId" & !species[i] %in% wadden_species) {
endYear = endYearWintering$IId
startYear = startYearWintering$IId
} else if (j == "IId" & species[i] %in% wadden_species) {
endYear = endYearWintering$IIdw
startYear = startYearWintering$IIdw
} else if (j == "IIe") {
endYear = endYearWintering$IIe
startYear = startYearWintering$IIe
} else if (j == "IIf") {
endYear = endYearWintering$IIf
startYear = startYearWintering$IIf
}
dat_unsm[[str_glue("{species[i]}_{j}")]] <-
calculate_non_breeding_relative_abundance(ospar = subset(new_dat, SUBADMIN == j),
label = j,
species_name = species[i],
target_threshold = target,
start_year = startYear,
end_year = endYear)
}
}
dat_unsm <- bind_rows(dat_unsm)
# close baseline regression results log
luzlogr::closelog(sessionInfo = FALSE)
Format dat_unsm dataframe into same format used by BTO code in the subsequent analyses.
names(dat_unsm)[names(dat_unsm) == "species"] <- "Common_name"
dat_unsm <- merge(dat_unsm, add.info, all.x = TRUE)
dat_unsm$Region_Subdivision <- str_remove(dat_unsm$region, "OSPAR")
dat_unsm$Reg_sub_SpID <- do.call(paste, dat_unsm[c("Region_Subdivision","SpeciesID")])
dat_unsm <- dat_unsm[, c("Region_Subdivision", "Common_name", "SpeciesID", "years", "Reg_sub_SpID", "index", "six_year_mean")]
names(dat_unsm) <- c("Region_Subdivision", "CommonName", "SpeciesID", "Year", "Reg_sub_SpID", "Index", "SixYrMean")
Export dat_unsm data frame with yearly regional abundance values, WinteringWaterbird_RegionalIndices.csv
non_breeding_regional_indices <- dat_unsm[, c("Region_Subdivision", "CommonName", "SpeciesID", "Year", "Reg_sub_SpID", "Index")]
write.table(non_breeding_regional_indices, path("results", "EMECO checks", "non_breeding_out", "WinteringWaterbird_RegionalIndices.csv"), row.names = FALSE, sep=",")
Export dat_unsm data frame with six year regional abundance mean values, WinteringWaterbird_RegionalIndices_mean.csv
non_breeding_regional_indices_mean <- dat_unsm[, c("Region_Subdivision", "CommonName", "SpeciesID", "Year", "Reg_sub_SpID", "SixYrMean")]
write.table(non_breeding_regional_indices_mean, path("results", "EMECO checks", "non_breeding_out", "WinteringWaterbird_RegionalIndices_mean.csv"), row.names = FALSE, sep=",")
Tidy up removing all objects from global environment except object required in the subsequent analyses
required_data <- c("species", "add.info", "endYearWintering", "previousYearWintering",
"region_labels", "sub_region_labels", "region_sub_region_labels",
"wadden_species")
rm(list=ls()[!ls() %in% c("non_breeding_regional_indices_mean", required_data)])
From non_breeding_regional_indices_mean get the latest six year relative abundance mean for each species as rows and OSPAR region / sub-region as columns. including species number of eggs and functional group. Sorted by functional group and species scientific name.
# select the latest year
ByRegion.dat.latest.mean.I <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "I" & Year == endYearWintering$I)
ByRegion.dat.latest.mean.In <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "In" & Year == endYearWintering$In)
ByRegion.dat.latest.mean.Io <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "Io" & Year == endYearWintering$Io)
ByRegion.dat.latest.mean.II <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "II" & Year == endYearWintering$II &
!CommonName %in% wadden_species)
ByRegion.dat.latest.mean.IIw <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "II" & Year == endYearWintering$IIw &
CommonName %in% wadden_species)
ByRegion.dat.latest.mean.IIa <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIa" & Year == endYearWintering$IIa)
ByRegion.dat.latest.mean.IIb <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIb" & Year == endYearWintering$IIb)
ByRegion.dat.latest.mean.IIc <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIc" & Year == endYearWintering$IIc)
ByRegion.dat.latest.mean.IId <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IId" & Year == endYearWintering$IId &
!CommonName %in% wadden_species)
ByRegion.dat.latest.mean.IIdw <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IId" & Year == endYearWintering$IIdw &
CommonName %in% wadden_species)
ByRegion.dat.latest.mean.IIe <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIe" & Year == endYearWintering$IIe)
ByRegion.dat.latest.mean.IIf <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIf" & Year == endYearWintering$IIf)
ByRegion.dat.latest.mean.III <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "III" & Year == endYearWintering$III)
ByRegion.dat.latest.mean <- bind_rows(ByRegion.dat.latest.mean.I,
ByRegion.dat.latest.mean.In,
ByRegion.dat.latest.mean.Io,
ByRegion.dat.latest.mean.II,
ByRegion.dat.latest.mean.IIw,
ByRegion.dat.latest.mean.IIa,
ByRegion.dat.latest.mean.IIb,
ByRegion.dat.latest.mean.IIc,
ByRegion.dat.latest.mean.IId,
ByRegion.dat.latest.mean.IIdw,
ByRegion.dat.latest.mean.IIe,
ByRegion.dat.latest.mean.IIf,
ByRegion.dat.latest.mean.III)
# format region
ByRegion.dat.latest.mean$region <- str_glue("OSPAR{ByRegion.dat.latest.mean$Region_Subdivision}")
# select and rename columns
ByRegion.dat.latest.mean <- ByRegion.dat.latest.mean[, c("CommonName", "region", "SixYrMean")]
names(ByRegion.dat.latest.mean) <- c("species", "region", "SixYrMean")
# pivot data frame so species in row and regions+sub_regions in columns
ByRegion.dat.latest.mean <- ByRegion.dat.latest.mean %>%
pivot_wider(names_from = "region", values_from = "SixYrMean")
# add species information and sort by functional group and species name
ByRegion.dat.latest.mean <- merge(ByRegion.dat.latest.mean, add.info, by.x = "species", by.y = "Common_name", all.x = TRUE)
ByRegion.dat.latest.mean <- ByRegion.dat.latest.mean[order(ByRegion.dat.latest.mean$Fun_Group, ByRegion.dat.latest.mean$Species),]
Export species+region/subregion latest six year relative abundance means
write.csv(ByRegion.dat.latest.mean,
path("results", "EMECO checks", "non_breeding_out", "Wintering Seabirds Indicator Latest Six Year Mean.csv"), row.names = FALSE)
From non_breeding_regional_indices_mean get the previous six year relative abundance mean for each species as rows and OSPAR region / sub-region as columns. including species number of eggs and functional group. Sorted by functional group and species scientific name.
# select the previous year
ByRegion.dat.previous.mean.I <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "I" & Year == previousYearWintering)
ByRegion.dat.previous.mean.In <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "In" & Year == previousYearWintering)
ByRegion.dat.previous.mean.Io <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "Io" & Year == previousYearWintering)
ByRegion.dat.previous.mean.II <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "II" & Year == previousYearWintering &
!CommonName %in% wadden_species)
ByRegion.dat.previous.mean.IIw <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "II" & Year == previousYearWintering &
CommonName %in% wadden_species)
ByRegion.dat.previous.mean.IIa <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIa" & Year == previousYearWintering)
ByRegion.dat.previous.mean.IIb <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIb" & Year == previousYearWintering)
ByRegion.dat.previous.mean.IIc <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIc" & Year == previousYearWintering)
ByRegion.dat.previous.mean.IId <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IId" & Year == previousYearWintering &
!CommonName %in% wadden_species)
ByRegion.dat.previous.mean.IIdw <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IId" & Year == previousYearWintering &
CommonName %in% wadden_species)
ByRegion.dat.previous.mean.IIe <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIe" & Year == previousYearWintering)
ByRegion.dat.previous.mean.IIf <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "IIf" & Year == previousYearWintering)
ByRegion.dat.previous.mean.III <- subset(non_breeding_regional_indices_mean,
Region_Subdivision == "III" & Year == previousYearWintering)
ByRegion.dat.previous.mean <- bind_rows(ByRegion.dat.previous.mean.I,
ByRegion.dat.previous.mean.In,
ByRegion.dat.previous.mean.Io,
ByRegion.dat.previous.mean.II,
ByRegion.dat.previous.mean.IIw,
ByRegion.dat.previous.mean.IIa,
ByRegion.dat.previous.mean.IIb,
ByRegion.dat.previous.mean.IIc,
ByRegion.dat.previous.mean.IId,
ByRegion.dat.previous.mean.IIdw,
ByRegion.dat.previous.mean.IIe,
ByRegion.dat.previous.mean.IIf,
ByRegion.dat.previous.mean.III)
# format region
ByRegion.dat.previous.mean$region <- str_glue("OSPAR{ByRegion.dat.previous.mean$Region_Subdivision}")
# select and rename columns
ByRegion.dat.previous.mean <- ByRegion.dat.previous.mean[, c("CommonName", "region", "SixYrMean")]
names(ByRegion.dat.previous.mean) <- c("species", "region", "SixYrMean")
# pivot data frame so species in row and regions+sub_regions in columns
ByRegion.dat.previous.mean <- ByRegion.dat.previous.mean %>%
pivot_wider(names_from = "region", values_from = "SixYrMean")
# add species information and sort by functional group and species name
ByRegion.dat.previous.mean <- merge(ByRegion.dat.previous.mean, add.info, by.x = "species", by.y = "Common_name", all.x = TRUE)
ByRegion.dat.previous.mean <- ByRegion.dat.previous.mean[order(ByRegion.dat.previous.mean$Fun_Group, ByRegion.dat.previous.mean$Species),]
Export species+region/subregion latest six year relative abundance means
write.csv(ByRegion.dat.previous.mean,
path("results", "EMECO checks", "non_breeding_out",
str_glue("Wintering Seabirds Indicator {previousYearWintering} Six Year Mean.csv")), row.names = FALSE)
Create latest non-breeding assessment traffic light plots, with colour indicating if latest six year mean is above or below the relative abundance threshold for the species and arrows indicating indicating if there has been a change in threshold status between the latest year and previous year (= 2014) assessments.