This script, derived from the 2016/17 BTO scripts, incorporates the code updates required for 2021. The script

  • Calculates and plots the relative abundance and six yearly mean relative abundances trends for each wintering non-breeding species, across each OSPAR region and subregion
  • Plots the species latest six year mean relative abundance value compared to its threshold value and compares this latest threshold status to a previous years threshold status, to produce a traffic light plot for the species within the OSPAR regions and birds subdivisions.

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

  • set the target threshold, depending on number of eggs species lays
  • combine the abundance counts across all regions in to new_dat dataframe, converting to wide format with years as columns and exporting as csv file, eg Arctic_Skua.csv
  • export the species abundance data eg Arctic Skua_non_breeding.csv
  • calculate the yearly and six yearly mean relative abundances and plotting relative abundance trends by running calculate_non_breeding_relative_abundance function for each region.
  • calculate the yearly and six yearly mean relative abundances and plotting relative abundance trends by running calculate_non_breeding_relative_abundance function for each sub region.
  • yearly and six year mean relative abundance values for each species+region/sub region added to dat_unsm data frame.
# 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.