This script plotting and outputting the proportion of species within regions and subdivisions and functional groups with latest six year relative abundance mean above threshold, was derived from the 2016/17 BTO scripts, incorporating updates required for 2021.

startYearBreeding <- list(
  Im = 1991,
  In = 1991,
  Io = 1991,
  I = 1991,
  IIa = 1991,
  IIb = 1991,
  IIc = 1991,
  IId = 1991,
  IIe = 1991,
  IIf = 1991,
  II = 1991,
  III = 1991,
  IV = 1991
)
endYearBreeding <- list(
  Im = 2020,
  In = 2020,
  Io = 2020,
  I = 2020,
  IIa = 2019,
  IIb = 2020,
  IIc = 2020,
  IId = 2017,
  IIe = 2019,
  IIf = 2019,
  II = 2017,
  III = 2019,
  IV = 2016
)
startYearWintering <- list(
  In = 1991,
  Io = 1991,
  I = 1991,
  IIa = 1991,
  IIb = 1991,
  IIc = 1991,
  IId = 1991,
  IIe = 1991,
  IIf = 1991,
  II = 1991,
  III = 1991
)
endYearWintering <- list(
  In = 2020,
  Io = 2020,
  I = 2020,
  IIa = 2020,
  IIb = 2020,
  IIc = 2020,
  IId = 2016,
  IIe = 2020,
  IIf = 2020,
  II = 2016,
  III = 2020
)

# list of OSPAR regions and functional groups
region_labels_breeding <- c("I", "Im", "In", "Io",
                            "II", "IIa", "IIb", "IIc", "IId", "IIe", "IIf",
                            "III",
                            "IV")
region_labels_wintering <- c("I", "In", "Io",
                             "II", "IIa", "IIb", "IIc", "IId", "IIe", "IIf",
                             "III")
functional_groups = c("All Feeders", "Wading Feeder", "Surface Feeder", 
                     "Water Column Feeder", "Benthic Feeder", "Grazing Feeder")

# set theme for plot
theme_set(
  theme_classic() +
  theme(panel.background = element_rect(colour = "black"),
        axis.text = element_text(face = "bold"),
        panel.grid.major.y = element_line(colour = "grey", size = 0.5),
        legend.box.background = element_rect(colour = "black"),
        legend.title = NULL,
        legend.spacing.y = unit(0, "mm"),
        plot.margin = unit(c(10, 10, 10, 10), "mm"))
)

Create the directory structure for the results

dir.create(path("results", "summary out"))
dir.create(path("results", "EMECO checks", "summary out"))

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:]|-", "_")

Import breeding six year relative abundance means

get_breeding_relative_abundance_means <- function(label, start_year, end_year) {
  
  # import six year relative abundance means
  unsm_all <- read.table(path("results", "EMECO checks", "breeding_out", 
                              str_glue("OSPAR{label}_unsm_all_mean.csv")), 
                         header = TRUE, sep = ",", stringsAsFactors = FALSE)
  
  # start and end year columns
  start_column <- str_glue("X{start_year}")
  end_column <- str_glue("X{end_year}")
  
  # format data, create dataframe with species, N eggs laid and indicator status in each year for OSPAR region / sub_region
  unsm_dat = t(unsm_all)
  colnames(unsm_dat) = paste("X", as.numeric(unsm_dat[1,]), sep = "")
  unsm_dat = unsm_dat[-1,]
  unsm_dat = data.frame(unsm_dat)
  unsm_dat$spp = row.names(unsm_dat)
  unsm_dat$spp = str_replace_all(unsm_dat$spp, "//.", "_")
  unsm_dat = merge(unsm_dat, add.info, by.x = "spp", by.y = "Common_name", all.x = TRUE)
  unsm_dat <- unsm_dat %>%
    select(Species, spp, Eggs, Fun_Group, {{start_column}}:{{end_column}})
  
  # export relative_abundance_means
  write.csv(unsm_dat, path("results", "EMECO checks", "breeding_out", 
                                str_glue("OSPAR{label} Seabird Six Year Mean Abundance Values.csv")), 
          row.names = FALSE)
  
  # return relative_abundance_means
  return(unsm_dat)
}
unsm_dat_OSPARI <- get_breeding_relative_abundance_means("I", startYearBreeding$I,  endYearBreeding$I)
unsm_dat_OSPARIm <- get_breeding_relative_abundance_means("Im", startYearBreeding$Im,  endYearBreeding$Im)
unsm_dat_OSPARIn <- get_breeding_relative_abundance_means("In", startYearBreeding$In,  endYearBreeding$In)
unsm_dat_OSPARIo <- get_breeding_relative_abundance_means("Io", startYearBreeding$Io,  endYearBreeding$Io)
unsm_dat_OSPARII <- get_breeding_relative_abundance_means("II", startYearBreeding$II,  endYearBreeding$II)
unsm_dat_OSPARIIa <- get_breeding_relative_abundance_means("IIa", startYearBreeding$IIa,  endYearBreeding$IIa)
unsm_dat_OSPARIIb <- get_breeding_relative_abundance_means("IIb", startYearBreeding$IIb,  endYearBreeding$IIb)
unsm_dat_OSPARIIc <- get_breeding_relative_abundance_means("IIc", startYearBreeding$IIc,  endYearBreeding$IIc)
unsm_dat_OSPARIId <- get_breeding_relative_abundance_means("IId", startYearBreeding$IId,  endYearBreeding$IId)
unsm_dat_OSPARIIe <- get_breeding_relative_abundance_means("IIe", startYearBreeding$IIe,  endYearBreeding$IIe)
unsm_dat_OSPARIIf <- get_breeding_relative_abundance_means("IIf", startYearBreeding$IIf,  endYearBreeding$IIf)
unsm_dat_OSPARIII <- get_breeding_relative_abundance_means("III", startYearBreeding$III,  endYearBreeding$III)
unsm_dat_OSPARIV <- get_breeding_relative_abundance_means("IV", startYearBreeding$IV,  endYearBreeding$IV)

Import wintering six year relative abundance means

OSPAR I non-breeding

OSPAR.dat <- read.table(path("results", "EMECO checks", "non_breeding_out",  "WinteringWaterbird_RegionalIndices_mean.csv"), 
                         header = TRUE, sep = ",", stringsAsFactors = FALSE) 
get_non_breeding_relative_abundance_means <- function(label, start_year, end_year) {
  
  # start and end year columns
  start_column <- as.character(start_year)
  end_column <- as.character(end_year)
  
  # format data, create dataframe with species, N eggs laid and indicator status in each year for OSPAR region / sub_region
  dat <- OSPAR.dat[OSPAR.dat$Region_Subdivision == label, ]
  dat <- dat[, c("CommonName", "Year", "SixYrMean")]
  names(dat) <- c("Common_name", "Year", "SixYrMean")
  dat <- dat %>% 
    pivot_wider(names_from = "Year", values_from = "SixYrMean")
  dat$Common_name <- str_replace_all(dat$Common_name, " ", "_")
  dat = merge(dat, add.info, by = "Common_name", all.x = TRUE)
  dat <- dat %>%
     select(Species, Common_name, Eggs, Fun_Group, {{start_column}}:{{end_column}})
  
  # export relative_abundance_means
  write.csv(dat, path("results", "EMECO checks", "non_breeding_out",
                      str_glue("OSPAR{label} Wintering Waterbirds Six Year Mean Abundance Values.csv")))
            
    # return relative_abundance_means
  return(dat) 
}
OSPARI.dat <- get_non_breeding_relative_abundance_means("I", startYearWintering$I, endYearWintering$I)
OSPARIn.dat <- get_non_breeding_relative_abundance_means("In", startYearWintering$In, endYearWintering$In)
OSPARIo.dat <- get_non_breeding_relative_abundance_means("Io", startYearWintering$Io, endYearWintering$Io)
OSPARII.dat <- get_non_breeding_relative_abundance_means("II", startYearWintering$II, endYearWintering$II)
OSPARIIa.dat <- get_non_breeding_relative_abundance_means("IIa", startYearWintering$IIa, endYearWintering$IIa)
OSPARIIb.dat <- get_non_breeding_relative_abundance_means("IIb", startYearWintering$IIb, endYearWintering$IIb)
OSPARIIc.dat <- get_non_breeding_relative_abundance_means("IIc", startYearWintering$IIc, endYearWintering$IIc)
OSPARIId.dat <- get_non_breeding_relative_abundance_means("IId", startYearWintering$IId, endYearWintering$IId)
OSPARIIe.dat <- get_non_breeding_relative_abundance_means("IIe", startYearWintering$IIe, endYearWintering$IIe)
OSPARIIf.dat <- get_non_breeding_relative_abundance_means("IIf", startYearWintering$IIf, endYearWintering$IIf)
OSPARIII.dat <- get_non_breeding_relative_abundance_means("III", startYearWintering$III, endYearWintering$III)

Plot graphs of yearly trend of percentage of species achieving GES

Plot yearly trend of percentage of breeding and wintering species exceeding relative abundance threshold for each region

# OSPAR I ----------------------------------------------------------------------

# create empty vector to store yearly percentage of breeding and wintering species exceeding threshold
OSPARI.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:min(ncol(unsm_dat_OSPARI), ncol(OSPARI.dat))) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species + 
  ## total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARI.total = length(which(!is.na(unsm_dat_OSPARI[, i]))) + 
    length(which(!is.na(OSPARI.dat[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg +
  # the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARI.pass = length(which(unsm_dat_OSPARI[, i] > 0.7 & unsm_dat_OSPARI$Eggs > 1)) +
    length(which(unsm_dat_OSPARI[, i] > 0.8 & unsm_dat_OSPARI$Eggs == 1)) + 
    length(which(OSPARI.dat[, i] > 0.7 & OSPARI.dat$Eggs > 1)) +
    length(which(OSPARI.dat[, i] > 0.8 & OSPARI.dat$Eggs == 1))
  
  # add percentage of breeding and wintering species that exceed the threshold for that year
  OSPARI.prop.pass[i] = OSPARI.pass / OSPARI.total * 100
}

dat_I <- data.frame(Year = seq(max(startYearBreeding$I, startYearWintering$I) + 5, 
                               min(endYearBreeding$I, endYearWintering$I), 1),
                    Region = "OSPARI",
                    Proportion = na.omit(OSPARI.prop.pass))

# OSPAR II ---------------------------------------------------------------------

# create empty vector to store yearly percentage of breeding and wintering species exceeding threshold
OSPARII.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:min(ncol(unsm_dat_OSPARII), ncol(OSPARII.dat))) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species + 
  ## total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARII.total = length(which(!is.na(unsm_dat_OSPARII[, i]))) + 
    length(which(!is.na(OSPARII.dat[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg +
  # the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARII.pass = length(which(unsm_dat_OSPARII[, i] > 0.7 & unsm_dat_OSPARII$Eggs > 1)) +
    length(which(unsm_dat_OSPARII[, i] > 0.8 & unsm_dat_OSPARII$Eggs == 1)) + 
    length(which(OSPARII.dat[, i] > 0.7 & OSPARII.dat$Eggs > 1)) +
    length(which(OSPARII.dat[, i] > 0.8 & OSPARII.dat$Eggs == 1))
  
  # add percentage of breeding and wintering species that exceed the threshold for that year
  OSPARII.prop.pass[i] = OSPARII.pass / OSPARII.total * 100
}

dat_II <- data.frame(Year = seq(max(startYearBreeding$II, startYearWintering$II) + 5, 
                                min(endYearBreeding$II, endYearWintering$II), 1),
                     Region = "OSPARII",
                     Proportion = na.omit(OSPARII.prop.pass))

# OSPAR III ---------------------------------------------------------------------

# create empty vector to store yearly percentage of breeding and wintering species exceeding threshold
OSPARIII.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:min(ncol(unsm_dat_OSPARIII), ncol(OSPARIII.dat))) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species + 
  ## total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARIII.total = length(which(!is.na(unsm_dat_OSPARIII[, i]))) + 
    length(which(!is.na(OSPARIII.dat[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg +
  # the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARIII.pass = length(which(unsm_dat_OSPARIII[, i] > 0.7 & unsm_dat_OSPARIII$Eggs > 1)) +
    length(which(unsm_dat_OSPARIII[, i] > 0.8 & unsm_dat_OSPARIII$Eggs == 1)) + 
    length(which(OSPARIII.dat[, i] > 0.7 & OSPARIII.dat$Eggs > 1)) +
    length(which(OSPARIII.dat[, i] > 0.8 & OSPARIII.dat$Eggs == 1))
  
  # add percentage of breeding and wintering species that exceed the threshold for that year
  OSPARIII.prop.pass[i] = OSPARIII.pass / OSPARIII.total * 100
}

dat_III <- data.frame(Year = seq(max(startYearBreeding$III, startYearWintering$III) + 5, 
                                 min(endYearBreeding$III, endYearWintering$III), 1),
                      Region = "OSPARIII",
                      Proportion = na.omit(OSPARIII.prop.pass))

Plot

# set the x axis year range expanded to 5 year intervals
xAxisMinYear <- plyr::round_any(reduce(append(startYearBreeding, startYearWintering), min), 5)
xAxisMaxYear <- plyr::round_any(reduce(append(endYearBreeding, endYearWintering), max), 5, f = ceiling)

# create plot
plot <- ggplot() +
  theme(legend.position = c(0.15, 0.1)) +
  coord_cartesian(xlim = c(xAxisMinYear, xAxisMaxYear),
                  ylim = c(40, 100)) +
  scale_x_continuous(breaks = seq(xAxisMinYear, xAxisMaxYear, 5),
                     name = NULL) +
  scale_y_continuous(breaks = seq(40, 100, 10),
                     name = "Percentage of species with relative abundance above threshold value") +
  scale_colour_manual(name = NULL,
                      breaks = c("I", "II", "III"),
                      values = c("#fde725",  "#440154", "#21918c"),
                      labels = c("Arctic Waters", "Greater North Sea", "Celtic Seas")) + 
  # 75% cut off line
  geom_hline(yintercept = 75, linetype = "solid", colour = "black", size = 1.2) +
   # add percentage of breeding and wintering species exceeding threshold line +
  geom_line(data = dat_I,
            mapping = aes(Year, Proportion, colour = "I"), size = 1.5) +
  geom_line(data = dat_II,
            mapping = aes(Year, Proportion, colour = "II"), size = 1.5) +
  geom_line(data = dat_III,
            mapping = aes(Year, Proportion, colour = "III"), size = 1.5) +
  annotate(geom = "text", x = xAxisMinYear, y = 74, 
           label = "Threshold value", color = "black", size = 3, hjust = 0.1)

# save plot
ggsave(plot = plot,
       filename = path("results", "summary out", "Trends_All.png"),
       width = 2100, height = 2100, units = "px")

Export results

### output data to check with EMECO
dat <- bind_rows(dat_I, dat_II, dat_III) %>%
  pivot_wider(names_from = Region, values_from = Proportion)
write.csv(dat, path("results", "EMECO checks", "summary out", "Proportion All Species Achieving GES.csv"),
          row.names = FALSE)

Plot yearly trend of percentage of breeding species exceeding relative abundance threshold for each region

# OSPAR I ----------------------------------------------------------------------

# create empty vector to store yearly percentage of breeding species exceeding threshold
OSPARI.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARI)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species 
  ### with six year relative abundance mean scores for the year
  OSPARI.total = length(which(!is.na(unsm_dat_OSPARI[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg
  OSPARI.pass = length(which(unsm_dat_OSPARI[, i] > 0.7 & unsm_dat_OSPARI$Eggs > 1)) +
    length(which(unsm_dat_OSPARI[, i] > 0.8 & unsm_dat_OSPARI$Eggs == 1))
  
  # add percentage of breeding species that exceed the threshold for that year
  OSPARI.prop.pass[i] = OSPARI.pass / OSPARI.total * 100
}

dat_I <- data.frame(Year = seq(startYearBreeding$I + 5, endYearBreeding$I, 1),
                    Region = "OSPARI",
                    Proportion = na.omit(OSPARI.prop.pass))

# OSPAR II ---------------------------------------------------------------------

# create empty vector to store yearly percentage of breeding species exceeding threshold
OSPARII.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARII)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species
  ### with six year relative abundance mean scores for the year
  OSPARII.total = length(which(!is.na(unsm_dat_OSPARII[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg
  OSPARII.pass = length(which(unsm_dat_OSPARII[, i] > 0.7 & unsm_dat_OSPARII$Eggs > 1)) +
    length(which(unsm_dat_OSPARII[, i] > 0.8 & unsm_dat_OSPARII$Eggs == 1))
  
  # add percentage of breeding species that exceed the threshold for that year
  OSPARII.prop.pass[i] = OSPARII.pass / OSPARII.total * 100
}

dat_II <- data.frame(Year = seq(startYearBreeding$II + 5, endYearBreeding$II, 1),
                     Region = "OSPARII",
                     Proportion = na.omit(OSPARII.prop.pass))

# OSPAR III ---------------------------------------------------------------------

# create empty vector to store yearly percentage of breeding species exceeding threshold
OSPARIII.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARIII)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species
  ### with six year relative abundance mean scores for the year
  OSPARIII.total = length(which(!is.na(unsm_dat_OSPARIII[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg
  OSPARIII.pass = length(which(unsm_dat_OSPARIII[, i] > 0.7 & unsm_dat_OSPARIII$Eggs > 1)) +
    length(which(unsm_dat_OSPARIII[, i] > 0.8 & unsm_dat_OSPARIII$Eggs == 1))
  
  # add percentage of breeding species that exceed the threshold for that year
  OSPARIII.prop.pass[i] = OSPARIII.pass / OSPARIII.total * 100
}

dat_III <- data.frame(Year = seq(startYearBreeding$III + 5, endYearBreeding$III, 1),
                      Region = "OSPARIII",
                      Proportion = na.omit(OSPARIII.prop.pass))

# OSPAR IV ----------------------------------------------------------------------

# create empty vector to store yearly percentage of breeding species exceeding threshold
OSPARIV.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARIV)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species
  ### with six year relative abundance mean scores for the year
  OSPARIV.total = length(which(!is.na(unsm_dat_OSPARIV[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg
  OSPARIV.pass = length(which(unsm_dat_OSPARIV[, i] > 0.7 & unsm_dat_OSPARIV$Eggs > 1)) +
    length(which(unsm_dat_OSPARIV[, i] > 0.8 & unsm_dat_OSPARIV$Eggs == 1))
  
  # add percentage of breeding species that exceed the threshold for that year
  OSPARIV.prop.pass[i] = OSPARIV.pass / OSPARIV.total * 100
}

dat_IV <- data.frame(Year =seq(startYearBreeding$IV + 5, endYearBreeding$IV, 1),
                     Region = "OSPARIV",
                     Proportion = na.omit(OSPARIV.prop.pass))

Plot

# set the x axis year range expanded to 5 year intervals
xAxisMinYear <- plyr::round_any(reduce(startYearBreeding, min), 5)
xAxisMaxYear <- plyr::round_any(reduce(endYearBreeding, max), 5, f = ceiling)

# create plot
plot <- ggplot() +
  theme(legend.position = c(0.22, 0.12)) +
  coord_cartesian(xlim = c(xAxisMinYear, xAxisMaxYear),
                  ylim = c(40, 100)) +
  scale_x_continuous(breaks = seq(xAxisMinYear, xAxisMaxYear, 5),
                     name = NULL) +
  scale_y_continuous(breaks = seq(40, 100, 10),
                     name = "Percentage of species with relative abundance above threshold value") +
  scale_colour_manual(name = NULL,
                      breaks = c("I", "II", "III", "IV"),
                      values = c("#fde725",  "#440154", "#21918c", "#5ec962"),
                      labels = c("Arctic Waters", "Greater North Sea", "Celtic Seas", "Bay of Biscay and Iberian Coast")) + # "Wider Atlantic"
  # 75% cut off line
  geom_hline(yintercept = 75, linetype = "solid", colour = "black", size = 1.2) +
   # add percentage of breeding and wintering species exceeding threshold line +
  geom_line(data = dat_I,
            mapping = aes(Year, Proportion, colour = "I"), size = 1.5) +
  geom_line(data = dat_II,
            mapping = aes(Year, Proportion, colour = "II"), size = 1.5) +
  geom_line(data = dat_III,
            mapping = aes(Year, Proportion, colour = "III"), size = 1.5) +
  geom_line(data = dat_IV,
            mapping = aes(Year, Proportion, colour = "IV"), size = 1.5) +
  annotate(geom = "text", x = xAxisMinYear, y = 74, 
           label = "Threshold value", color = "black", size = 3, hjust = 0.1)

# save plot
ggsave(plot = plot,
       filename = path("results", "summary out", "Trends_Breeding.png"),
       width = 2100, height = 2100, units = "px")

Export results

### output data to check with EMECO
dat <- bind_rows(dat_I, dat_II, dat_III, dat_IV) %>%
  pivot_wider(names_from = Region, values_from = Proportion)
write.csv(dat, path("results", "EMECO checks", "summary out", "Proportion Breeding Species Achieving GES.csv"),
          row.names = FALSE)

Plot yearly trend of percentage of wintering species exceeding relative abundance threshold for each region

# OSPAR I ----------------------------------------------------------------------

# create empty vector to store yearly percentage of wintering species exceeding threshold
OSPARI.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARI.dat)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARI.total = length(which(!is.na(OSPARI.dat[, i])))
  
  # count the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARI.pass =  length(which(OSPARI.dat[, i] > 0.7 & OSPARI.dat$Eggs > 1)) +
    length(which(OSPARI.dat[, i] > 0.8 & OSPARI.dat$Eggs == 1))
  
  # add percentage of wintering species that exceed the threshold for that year
  OSPARI.prop.pass[i] = OSPARI.pass / OSPARI.total * 100
}

dat_I <- data.frame(Year = seq(startYearWintering$I + 5, endYearWintering$I, 1),
                    Region = "OSPARI",
                    Proportion = na.omit(OSPARI.prop.pass))

# OSPAR II ---------------------------------------------------------------------

# create empty vector to store yearly percentage of wintering species exceeding threshold
OSPARII.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARII.dat)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARII.total = length(which(!is.na(OSPARII.dat[, i])))
  
  # count the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARII.pass = length(which(OSPARII.dat[, i] > 0.7 & OSPARII.dat$Eggs > 1)) +
    length(which(OSPARII.dat[, i] > 0.8 & OSPARII.dat$Eggs == 1))
  
  # add percentage of wintering species that exceed the threshold for that year
  OSPARII.prop.pass[i] = OSPARII.pass / OSPARII.total * 100
}

dat_II <- data.frame(Year = seq(startYearWintering$II + 5, endYearWintering$II, 1),
                     Region = "OSPARII",
                     Proportion = na.omit(OSPARII.prop.pass))

# OSPAR III ---------------------------------------------------------------------

# create empty vector to store yearly percentage of wintering species exceeding threshold
OSPARIII.prop.pass = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARIII.dat)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARIII.total = length(which(!is.na(OSPARIII.dat[, i])))
  
  # count the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARIII.pass = length(which(OSPARIII.dat[, i] > 0.7 & OSPARIII.dat$Eggs > 1)) +
    length(which(OSPARIII.dat[, i] > 0.8 & OSPARIII.dat$Eggs == 1))
  
  # add percentage of wintering species that exceed the threshold for that year
  OSPARIII.prop.pass[i] = OSPARIII.pass / OSPARIII.total * 100
}

dat_III <- data.frame(Year = seq(startYearWintering$III + 5, endYearWintering$III, 1),
                      Region = "OSPARIII",
                      Proportion = na.omit(OSPARIII.prop.pass))

Plot

# set the x axis year range expanded to 5 year intervals
xAxisMinYear <- plyr::round_any(reduce(startYearWintering, min), 5)
xAxisMaxYear <- plyr::round_any(reduce(endYearWintering, max), 5, f = ceiling)

# create plot
plot <- ggplot() +
  theme(legend.position = c(0.15, 0.1)) +
  coord_cartesian(xlim = c(xAxisMinYear, xAxisMaxYear),
                  ylim = c(40, 100)) +
  scale_x_continuous(breaks = seq(xAxisMinYear, xAxisMaxYear, 5),
                     name = NULL) +
  scale_y_continuous(breaks = seq(40, 100, 10),
                     name = "Percentage of species with relative abundance above threshold value") +
  scale_colour_manual(name = NULL,
                      breaks = c("I", "II", "III"),
                      values = c("#fde725",  "#440154", "#21918c"),
                      labels = c("Arctic Waters", "Greater North Sea", "Celtic Seas")) +
  # 75% cut off line
  geom_hline(yintercept = 75, linetype = "solid", colour = "black", size = 1.2) +
   # add percentage of breeding and wintering species exceeding threshold line +
  geom_line(data = dat_I,
            mapping = aes(Year, Proportion, colour = "I"), size = 1.5) +
  geom_line(data = dat_II,
            mapping = aes(Year, Proportion, colour = "II"), size = 1.5) +
  geom_line(data = dat_III,
            mapping = aes(Year, Proportion, colour = "III"), size = 1.5) +
  annotate(geom = "text", x = xAxisMinYear, y = 74, 
           label = "Threshold value", color = "black", size = 3, hjust = 0.1)

# save plot
ggsave(plot = plot,
       filename = path("results", "summary out", "Trends_Wintering.png"),
       width = 2100, height = 2100, units = "px")

Export results

### output data to check with EMECO
dat <- bind_rows(dat_I, dat_II, dat_III) %>% #, dat_IV dat_V)
  pivot_wider(names_from = Region, values_from = Proportion)
write.csv(dat, path("results", "EMECO checks", "summary out", "Proportion Wintering Species Achieving GES.csv"),
          row.names = FALSE)

Plot yearly trend of percentage of species exceeding relative abundance threshold for OSPAR I region

# OSPAR I breeding -------------------------------------------------------------

# create empty vector to store yearly percentage of breeding species exceeding threshold
OSPARI.prop.pass.breeding = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARI)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species 
  ### with six year relative abundance mean scores for the year
  OSPARI.total = length(which(!is.na(unsm_dat_OSPARI[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg
  OSPARI.pass = length(which(unsm_dat_OSPARI[, i] > 0.7 & unsm_dat_OSPARI$Eggs > 1)) +
    length(which(unsm_dat_OSPARI[, i] > 0.8 & unsm_dat_OSPARI$Eggs == 1))
  
  # add percentage of breeding species that exceed the threshold for that year
  OSPARI.prop.pass.breeding[i] = OSPARI.pass / OSPARI.total * 100
}

dat_I_breeding <- data.frame(Year = seq(startYearBreeding$I + 5, endYearBreeding$I, 1),
                    Group = "breeding",
                    Proportion = na.omit(OSPARI.prop.pass.breeding))

# OSPAR I breeding surface feeders ---------------------------------------------
unsm_dat_OSPARI_surface <- unsm_dat_OSPARI[unsm_dat_OSPARI$Fun_Group == "Surface Feeder", ]

# create empty vector to store yearly percentage of breeding surface feeder species exceeding threshold
OSPARI.prop.pass.surface = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARI_surface)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding surface feeder species 
  ### with six year relative abundance mean scores for the year
  OSPARI.total.surface = length(which(!is.na(unsm_dat_OSPARI_surface[, i])))
  
  # count the number of breeding surface feeder species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding surface feeder species exceeding 0.8 threshold for species laying 1 egg
  OSPARI.pass.surface = length(which(unsm_dat_OSPARI_surface[, i] > 0.7 & unsm_dat_OSPARI_surface$Eggs > 1)) +
    length(which(unsm_dat_OSPARI_surface[, i] > 0.8 & unsm_dat_OSPARI_surface$Eggs == 1))
  
  # add percentage of breeding surface feeder species that exceed the threshold for that year
  OSPARI.prop.pass.surface[i] = OSPARI.pass.surface / OSPARI.total.surface * 100
}

dat_I_surface <- data.frame(Year = seq(startYearBreeding$I + 5, endYearBreeding$I, 1),
                    Group = "breeding_surface_feeders",
                    Proportion = na.omit(OSPARI.prop.pass.surface))

# OSPAR I breeding water column feeders ----------------------------------------
unsm_dat_OSPARI_water_column <- unsm_dat_OSPARI[unsm_dat_OSPARI$Fun_Group == "Water Column Feeder", ]

# create empty vector to store yearly percentage of breeding water column species exceeding threshold
OSPARI.prop.pass.water.column = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARI_water_column)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding water column species 
  ### with six year relative abundance mean scores for the year
  OSPARI.total.water.column = length(which(!is.na(unsm_dat_OSPARI_water_column[, i])))
  
  # count the number of breeding water column species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding water column species exceeding 0.8 threshold for species laying 1 egg
  OSPARI.pass.water.column = length(which(unsm_dat_OSPARI_water_column[, i] > 0.7 & unsm_dat_OSPARI_water_column$Eggs > 1)) +
    length(which(unsm_dat_OSPARI_water_column[, i] > 0.8 & unsm_dat_OSPARI_water_column$Eggs == 1))
  
  # add percentage of breeding water column feeder species that exceed the threshold for that year
  OSPARI.prop.pass.water.column[i] = OSPARI.pass.water.column / OSPARI.total.water.column * 100
}

dat_I_water_column <- data.frame(Year = seq(startYearBreeding$I + 5, endYearBreeding$I, 1),
                    Group = "breeding_water_column_feeders",
                    Proportion = na.omit(OSPARI.prop.pass.water.column))

# OSPAR I wintering ------------------------------------------------------------

# create empty vector to store yearly percentage of wintering species exceeding threshold
OSPARI.prop.pass.wintering = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARI.dat)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARI.total.wintering = length(which(!is.na(OSPARI.dat[, i])))
  
  # count the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARI.pass.wintering =  length(which(OSPARI.dat[, i] > 0.7 & OSPARI.dat$Eggs > 1)) +
    length(which(OSPARI.dat[, i] > 0.8 & OSPARI.dat$Eggs == 1))
  
  # add percentage of wintering species that exceed the threshold for that year
  OSPARI.prop.pass.wintering[i] = OSPARI.pass.wintering / OSPARI.total.wintering * 100
}

dat_I_wintering <- data.frame(Year = seq(startYearWintering$I + 5, endYearWintering$I, 1),
                    Group = "non-breeding",
                    Proportion = na.omit(OSPARI.prop.pass.wintering))

# OSPAR I wintering wading feeder ----------------------------------------------
OSPARI.dat.wader <- OSPARI.dat[OSPARI.dat$Fun_Group == "Wading Feeder", ]

# create empty vector to store yearly percentage of wintering wading feeder species exceeding threshold
OSPARI.prop.pass.wader = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARI.dat.wader)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering wader feeding species 
  ### with six year relative abundance mean scores for the year
  OSPARI.total.wader = length(which(!is.na(OSPARI.dat.wader[, i])))
  
  # count the number of wintering wading feeder species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering wading feeder species exceeding 0.8 threshold for species laying 1 egg 
  OSPARI.pass.wader =  length(which(OSPARI.dat.wader[, i] > 0.7 & OSPARI.dat.wader$Eggs > 1)) +
    length(which(OSPARI.dat.wader[, i] > 0.8 & OSPARI.dat.wader$Eggs == 1))
  
  # add percentage of wintering wading feeder species that exceed the threshold for that year
  OSPARI.prop.pass.wader[i] = OSPARI.pass.wader / OSPARI.total.wader * 100
}

Plot

# set the x axis year range expanded to 5 year intervals
xAxisMinYear <- plyr::round_any(reduce(append(startYearBreeding$I, startYearWintering$I), min), 5)
xAxisMaxYear <- plyr::round_any(reduce(append(endYearBreeding$I, endYearWintering$I), max), 5, f = ceiling)

# create plot
plot <- ggplot() +
  theme(legend.position = c(0.25, 0.13),
        legend.key.width = unit(20, "mm")) +
  coord_cartesian(xlim = c(xAxisMinYear, xAxisMaxYear),
                  ylim = c(40, 100)) +
  scale_x_continuous(breaks = seq(xAxisMinYear, xAxisMaxYear, 5),
                     name = NULL) +
  scale_y_continuous(breaks = seq(40, 100, 10),
                     name = "Percentage of species with relative abundance above threshold value") +
  scale_colour_manual(name = NULL,
                      breaks = c("breeding", "surface", "water_column", "wintering"),
                      values = c("#3b528b", "#3b528b", "#3b528b", "#5ec962"),
                      labels = c("all breeding", "breeding surface feeders", "breeding water column feeders", "all non-breeding")) +
  scale_linetype_manual(name = NULL,
                      breaks = c("breeding", "surface", "water_column", "wintering"),
                      values = c("solid", "dashed", "dotted", "solid"),
                      labels = c("all breeding", "breeding surface feeders", "breeding water column feeders", "all non-breeding")) +
  # 75% cut off line
  geom_hline(yintercept = 75, linetype = "solid", colour = "black", size = 1.2) +
  geom_line(data = dat_I_breeding,
            mapping = aes(Year, Proportion, linetype = "breeding", colour = "breeding"),  size = 1.5) +
  geom_line(data = dat_I_surface,
            mapping = aes(Year, Proportion, linetype = "surface", colour = "surface"), size = 1.5) +
  geom_line(data = dat_I_water_column,
            mapping = aes(Year, Proportion, linetype = "water_column", colour = "water_column"), size = 1.5) +
  geom_line(data = dat_I_wintering,
            mapping = aes(Year, Proportion, linetype = "wintering", colour = "wintering"), size = 1.5) +
  annotate(geom = "text", x = xAxisMinYear, y = 73, 
           label = "Threshold value", color = "black", size = 3, hjust = 0.1)

# save plot
ggsave(plot = plot,
       filename = path("results", "summary out", "Trends_OSPARI.png"),
       width = 2100, height = 2100, units = "px")

Export results

### output data to check with EMECO
dat <- bind_rows(dat_I_breeding, dat_I_surface, dat_I_water_column, dat_I_wintering) %>%
  pivot_wider(names_from = Group, values_from = Proportion)
write.csv(dat, path("results", "EMECO checks", "summary out", "Proportion OSPARI Species Achieving GES.csv"),
          row.names = FALSE)

Plot yearly trend of percentage of species exceeding relative abundance threshold for OSPAR II region

# OSPAR II breeding -------------------------------------------------------------

# create empty vector to store yearly percentage of breeding species exceeding threshold
OSPARII.prop.pass.breeding = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARII)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species 
  ### with six year relative abundance mean scores for the year
  OSPARII.total = length(which(!is.na(unsm_dat_OSPARII[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg
  OSPARII.pass = length(which(unsm_dat_OSPARII[, i] > 0.7 & unsm_dat_OSPARII$Eggs > 1)) +
    length(which(unsm_dat_OSPARII[, i] > 0.8 & unsm_dat_OSPARII$Eggs == 1))
  
  # add percentage of breeding species that exceed the threshold for that year
  OSPARII.prop.pass.breeding[i] = OSPARII.pass / OSPARII.total * 100
}

dat_II_breeding <- data.frame(Year = seq(startYearBreeding$II + 5, endYearBreeding$II, 1),
                    Group = "breeding",
                    Proportion = na.omit(OSPARII.prop.pass.breeding))

# OSPAR II breeding surface feeders ---------------------------------------------
unsm_dat_OSPARII_surface <- unsm_dat_OSPARII[unsm_dat_OSPARII$Fun_Group == "Surface Feeder", ]

# create empty vector to store yearly percentage of breeding surface feeder species exceeding threshold
OSPARII.prop.pass.surface = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARII_surface)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding surface feeder species 
  ### with six year relative abundance mean scores for the year
  OSPARII.total.surface = length(which(!is.na(unsm_dat_OSPARII_surface[, i])))
  
  # count the number of breeding surface feeder species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding surface feeder species exceeding 0.8 threshold for species laying 1 egg
  OSPARII.pass.surface = length(which(unsm_dat_OSPARII_surface[, i] > 0.7 & unsm_dat_OSPARII_surface$Eggs > 1)) +
    length(which(unsm_dat_OSPARII_surface[, i] > 0.8 & unsm_dat_OSPARII_surface$Eggs == 1))
  
  # add percentage of breeding surface feeder species that exceed the threshold for that year
  OSPARII.prop.pass.surface[i] = OSPARII.pass.surface / OSPARII.total.surface * 100
}

dat_II_surface <- data.frame(Year = seq(startYearBreeding$II + 5, endYearBreeding$II, 1),
                    Group = "breeding_surface_feeders",
                    Proportion = na.omit(OSPARII.prop.pass.surface))

# OSPAR II breeding water column feeders ----------------------------------------
unsm_dat_OSPARII_water_column <- unsm_dat_OSPARII[unsm_dat_OSPARII$Fun_Group == "Water Column Feeder", ]

# create empty vector to store yearly percentage of breeding water column species exceeding threshold
OSPARII.prop.pass.water.column = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARII_water_column)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding water column species 
  ### with six year relative abundance mean scores for the year
  OSPARII.total.water.column = length(which(!is.na(unsm_dat_OSPARII_water_column[, i])))
  
  # count the number of breeding water column species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding water column species exceeding 0.8 threshold for species laying 1 egg
  OSPARII.pass.water.column = length(which(unsm_dat_OSPARII_water_column[, i] > 0.7 & unsm_dat_OSPARII_water_column$Eggs > 1)) +
    length(which(unsm_dat_OSPARII_water_column[, i] > 0.8 & unsm_dat_OSPARII_water_column$Eggs == 1))
  
  # add percentage of breeding water column feeder species that exceed the threshold for that year
  OSPARII.prop.pass.water.column[i] = OSPARII.pass.water.column / OSPARII.total.water.column * 100
}

dat_II_water_column <- data.frame(Year = seq(startYearBreeding$II + 5, endYearBreeding$II, 1),
                    Group = "breeding_water_column_feeders",
                    Proportion = na.omit(OSPARII.prop.pass.water.column))

# OSPAR II wintering ------------------------------------------------------------

# create empty vector to store yearly percentage of wintering species exceeding threshold
OSPARII.prop.pass.wintering = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARII.dat)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARII.total.wintering = length(which(!is.na(OSPARII.dat[, i])))
  
  # count the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARII.pass.wintering =  length(which(OSPARII.dat[, i] > 0.7 & OSPARII.dat$Eggs > 1)) +
    length(which(OSPARII.dat[, i] > 0.8 & OSPARII.dat$Eggs == 1))
  
  # add percentage of wintering species that exceed the threshold for that year
  OSPARII.prop.pass.wintering[i] = OSPARII.pass.wintering / OSPARII.total.wintering * 100
}

dat_II_wintering <- data.frame(Year = seq(startYearWintering$II + 5, endYearWintering$II, 1),
                    Group = "non-breeding",
                    Proportion = na.omit(OSPARII.prop.pass.wintering))

# OSPAR II wintering wading feeder ----------------------------------------------
OSPARII.dat.wader <- OSPARII.dat[OSPARII.dat$Fun_Group == "Wading Feeder", ]

# create empty vector to store yearly percentage of wintering wading feeder species exceeding threshold
OSPARII.prop.pass.wader = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARII.dat.wader)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering wader feeding species 
  ### with six year relative abundance mean scores for the year
  OSPARII.total.wader = length(which(!is.na(OSPARII.dat.wader[, i])))
  
  # count the number of wintering wading feeder species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering wading feeder species exceeding 0.8 threshold for species laying 1 egg 
  OSPARII.pass.wader =  length(which(OSPARII.dat.wader[, i] > 0.7 & OSPARII.dat.wader$Eggs > 1)) +
    length(which(OSPARII.dat.wader[, i] > 0.8 & OSPARII.dat.wader$Eggs == 1))
  
  # add percentage of wintering wading feeder species that exceed the threshold for that year
  OSPARII.prop.pass.wader[i] = OSPARII.pass.wader / OSPARII.total.wader * 100
}

dat_II_wader <- data.frame(Year = seq(startYearWintering$II + 5, endYearWintering$II, 1),
                    Group = "non-breeding_wading_feeders",
                    Proportion = na.omit(OSPARII.prop.pass.wader))

Plot

# set the x axis year range expanded to 5 year intervals
xAxisMinYear <- plyr::round_any(reduce(append(startYearBreeding$II, startYearWintering$II), min), 5)
xAxisMaxYear <- plyr::round_any(reduce(append(endYearBreeding$II, endYearWintering$II), max), 5, f = ceiling)

# create plot
plot <- ggplot() +
  theme(legend.position = c(0.25, 0.13),
        legend.key.width = unit(20, "mm")) +
  coord_cartesian(xlim = c(xAxisMinYear, xAxisMaxYear),
                  ylim = c(20, 100)) +
  scale_x_continuous(breaks = seq(xAxisMinYear, xAxisMaxYear, 5),
                     name = NULL) +
  scale_y_continuous(breaks = seq(20, 100, 10),
                     name = "Percentage of species with relative abundance above threshold value") +
  scale_colour_manual(name = NULL,
                      breaks = c("breeding", "surface", "water_column", "wintering", "wader"),
                      values = c("#3b528b", "#3b528b", "#3b528b", "#5ec962", "#5ec962"), 
                      labels = c("all breeding", "breeding surface feeders", "breeding water column feeders", "all non-breeding", "non-breeding wading feeders")) +
  scale_linetype_manual(name = NULL,
                      breaks = c("breeding", "surface", "water_column", "wintering", "wader"),
                      values = c("solid", "dashed", "dotted", "solid", "dotdash"), 
                      labels = c("all breeding", "breeding surface feeders", "breeding water column feeders", "all non-breeding", "non-breeding wading feeders")) +
  # 75% cut off line
  geom_hline(yintercept = 75, linetype = "solid", colour = "black", size = 1.2) +
  geom_line(data = dat_II_breeding,
            mapping = aes(Year, Proportion, linetype = "breeding", colour = "breeding"),  size = 1.5) +
  geom_line(data = dat_II_surface,
            mapping = aes(Year, Proportion, linetype = "surface", colour = "surface"), size = 1.5) +
  geom_line(data = dat_II_water_column,
            mapping = aes(Year, Proportion, linetype = "water_column", colour = "water_column"), size = 1.5) +
  geom_line(data = dat_II_wintering,
            mapping = aes(Year, Proportion, linetype = "wintering", colour = "wintering"), size = 1.5) +
  geom_line(data = dat_II_wader,
            mapping = aes(Year, Proportion, linetype = "wader", colour = "wader"), size = 1.5) +
  # add percentage of breeding and wintering species exceeding threshold line +
  annotate(geom = "text", x = xAxisMinYear, y = 74, 
           label = "Threshold value", color = "black", size = 3, hjust = 0.1)

# save plot
ggsave(plot = plot,
       filename = path("results", "summary out", "Trends_OSPARII.png"),
       width = 2100, height = 2100, units = "px")

Export results

### output data to check with EMECO
dat <- bind_rows(dat_II_breeding, dat_II_surface, dat_II_water_column, dat_II_wintering, dat_II_wader) %>%
  pivot_wider(names_from = Group, values_from = Proportion)
write.csv(dat, path("results", "EMECO checks", "summary out", "Proportion OSPARII Species Achieving GES.csv"),
          row.names = FALSE)

Plot yearly trend of percentage of species exceeding relative abundance threshold for OSPAR III region

# OSPAR III breeding -------------------------------------------------------------

# create empty vector to store yearly percentage of breeding species exceeding threshold
OSPARIII.prop.pass.breeding = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARIII)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species 
  ### with six year relative abundance mean scores for the year
  OSPARIII.total = length(which(!is.na(unsm_dat_OSPARIII[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg
  OSPARIII.pass = length(which(unsm_dat_OSPARIII[, i] > 0.7 & unsm_dat_OSPARIII$Eggs > 1)) +
    length(which(unsm_dat_OSPARIII[, i] > 0.8 & unsm_dat_OSPARIII$Eggs == 1))
  
  # add percentage of breeding species that exceed the threshold for that year
  OSPARIII.prop.pass.breeding[i] = OSPARIII.pass / OSPARIII.total * 100
}

dat_III_breeding <- data.frame(Year = seq(startYearBreeding$III + 5, endYearBreeding$III, 1),
                    Group = "breeding",
                    Proportion = na.omit(OSPARIII.prop.pass.breeding))

# OSPAR III breeding surface feeders ---------------------------------------------
unsm_dat_OSPARIII_surface <- unsm_dat_OSPARIII[unsm_dat_OSPARIII$Fun_Group == "Surface Feeder", ]

# create empty vector to store yearly percentage of breeding surface feeder species exceeding threshold
OSPARIII.prop.pass.surface = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARIII_surface)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding surface feeder species 
  ### with six year relative abundance mean scores for the year
  OSPARIII.total.surface = length(which(!is.na(unsm_dat_OSPARIII_surface[, i])))
  
  # count the number of breeding surface feeder species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding surface feeder species exceeding 0.8 threshold for species laying 1 egg
  OSPARIII.pass.surface = length(which(unsm_dat_OSPARIII_surface[, i] > 0.7 & unsm_dat_OSPARIII_surface$Eggs > 1)) +
    length(which(unsm_dat_OSPARIII_surface[, i] > 0.8 & unsm_dat_OSPARIII_surface$Eggs == 1))
  
  # add percentage of breeding surface feeder species that exceed the threshold for that year
  OSPARIII.prop.pass.surface[i] = OSPARIII.pass.surface / OSPARIII.total.surface * 100
}

dat_III_surface <- data.frame(Year = seq(startYearBreeding$III + 5, endYearBreeding$III, 1),
                    Group = "breeding_surface_feeders",
                    Proportion = na.omit(OSPARIII.prop.pass.surface))

# OSPAR III breeding water column feeders ----------------------------------------
unsm_dat_OSPARIII_water_column <- unsm_dat_OSPARIII[unsm_dat_OSPARIII$Fun_Group == "Water Column Feeder", ]

# create empty vector to store yearly percentage of breeding water column species exceeding threshold
OSPARIII.prop.pass.water.column = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARIII_water_column)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding water column species 
  ### with six year relative abundance mean scores for the year
  OSPARIII.total.water.column = length(which(!is.na(unsm_dat_OSPARIII_water_column[, i])))
  
  # count the number of breeding water column species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding water column species exceeding 0.8 threshold for species laying 1 egg
  OSPARIII.pass.water.column = length(which(unsm_dat_OSPARIII_water_column[, i] > 0.7 & unsm_dat_OSPARIII_water_column$Eggs > 1)) +
    length(which(unsm_dat_OSPARIII_water_column[, i] > 0.8 & unsm_dat_OSPARIII_water_column$Eggs == 1))
  
  # add percentage of breeding water column feeder species that exceed the threshold for that year
  OSPARIII.prop.pass.water.column[i] = OSPARIII.pass.water.column / OSPARIII.total.water.column * 100
}

dat_III_water_column <- data.frame(Year = seq(startYearBreeding$III + 5, endYearBreeding$III, 1),
                    Group = "breeding_water_column_feeders",
                    Proportion = na.omit(OSPARIII.prop.pass.water.column))

# OSPAR III wintering ------------------------------------------------------------

# create empty vector to store yearly percentage of wintering species exceeding threshold
OSPARIII.prop.pass.wintering = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARIII.dat)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering species 
  ### with six year relative abundance mean scores for the year
  OSPARIII.total.wintering = length(which(!is.na(OSPARIII.dat[, i])))
  
  # count the number of wintering species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering species exceeding 0.8 threshold for species laying 1 egg 
  OSPARIII.pass.wintering =  length(which(OSPARIII.dat[, i] > 0.7 & OSPARIII.dat$Eggs > 1)) +
    length(which(OSPARIII.dat[, i] > 0.8 & OSPARIII.dat$Eggs == 1))
  
  # add percentage of wintering species that exceed the threshold for that year
  OSPARIII.prop.pass.wintering[i] = OSPARIII.pass.wintering / OSPARIII.total.wintering * 100
}

dat_III_wintering <- data.frame(Year = seq(startYearWintering$III + 5, endYearWintering$III, 1),
                    Group = "non-breeding",
                    Proportion = na.omit(OSPARIII.prop.pass.wintering))

# OSPAR III wintering wading feeder ----------------------------------------------
OSPARIII.dat.wader <- OSPARIII.dat[OSPARIII.dat$Fun_Group == "Wading Feeder", ]

# create empty vector to store yearly percentage of wintering wading feeder species exceeding threshold
OSPARIII.prop.pass.wader = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(OSPARIII.dat.wader)) {  # skip first 4 columns and first 5 years
  
  # count total number of wintering wader feeding species 
  ### with six year relative abundance mean scores for the year
  OSPARIII.total.wader = length(which(!is.na(OSPARIII.dat.wader[, i])))
  
  # count the number of wintering wading feeder species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of wintering wading feeder species exceeding 0.8 threshold for species laying 1 egg 
  OSPARIII.pass.wader =  length(which(OSPARIII.dat.wader[, i] > 0.7 & OSPARIII.dat.wader$Eggs > 1)) +
    length(which(OSPARIII.dat.wader[, i] > 0.8 & OSPARIII.dat.wader$Eggs == 1))
  
  # add percentage of wintering wading feeder species that exceed the threshold for that year
  OSPARIII.prop.pass.wader[i] = OSPARIII.pass.wader / OSPARIII.total.wader * 100
}

dat_III_wader <- data.frame(Year = seq(startYearWintering$III + 5, endYearWintering$III, 1),
                    Group = "non-breeding_wading_feeders",
                    Proportion = na.omit(OSPARIII.prop.pass.wader))

Plot

# set the x axis year range expanded to 5 year intervals
xAxisMinYear <- plyr::round_any(reduce(append(startYearBreeding$III, startYearWintering$III), min), 5)
xAxisMaxYear <- plyr::round_any(reduce(append(endYearBreeding$III, endYearWintering$III), max), 5, f = ceiling)

# create plot
plot <- ggplot() +
  theme(legend.position = c(0.25, 0.13),
        legend.key.width = unit(20, "mm")) +
  coord_cartesian(xlim = c(xAxisMinYear, xAxisMaxYear),
                  ylim = c(40, 100)) +
  scale_x_continuous(breaks = seq(xAxisMinYear, xAxisMaxYear, 5),
                     name = NULL) +
  scale_y_continuous(breaks = seq(40, 100, 10),
                     name = "Percentage of species with relative abundance above threshold value") +
  scale_colour_manual(name = NULL,
                      breaks = c("breeding", "surface", "water_column", "wintering", "wader"),
                      values = c("#3b528b", "#3b528b", "#3b528b", "#5ec962", "#5ec962"), 
                      labels = c("all breeding", "breeding surface feeders", "breeding water column feeders", "all non-breeding", "non-breeding wading feeders")) +
  scale_linetype_manual(name = NULL,
                      breaks = c("breeding", "surface", "water_column", "wintering", "wader"),
                      values = c("solid", "dashed", "dotted", "solid", "dotdash"), 
                      labels = c("all breeding", "breeding surface feeders", "breeding water column feeders", "all non-breeding", "non-breeding wading feeders")) +
  # 75% cut off line
  geom_hline(yintercept = 75, linetype = "solid", colour = "black", size = 1.2) +
  geom_line(data = dat_III_breeding,
            mapping = aes(Year, Proportion, linetype = "breeding", colour = "breeding"),  size = 1.5) +
  geom_line(data = dat_III_surface,
            mapping = aes(Year, Proportion, linetype = "surface", colour = "surface"), size = 1.5) +
  geom_line(data = dat_III_water_column,
            mapping = aes(Year, Proportion, linetype = "water_column", colour = "water_column"), size = 1.5) +
  geom_line(data = dat_III_wintering,
            mapping = aes(Year, Proportion, linetype = "wintering", colour = "wintering"), size = 1.5) +
  geom_line(data = dat_III_wader,
            mapping = aes(Year, Proportion, linetype = "wader", colour = "wader"), size = 1.5) +
  # add percentage of breeding and wintering species exceeding threshold line +
  annotate(geom = "text", x = xAxisMinYear, y = 74, 
           label = "Threshold value", color = "black", size = 3, hjust = 0.1)

# save plot
ggsave(plot = plot,
       filename = path("results", "summary out", "Trends_OSPARIII.png"),
       width = 2100, height = 2100, units = "px")

Export results

### output data to check with EMECO
dat <- bind_rows(dat_III_breeding, dat_III_surface, dat_III_water_column, dat_III_wintering, dat_III_wader) %>%
  pivot_wider(names_from = Group, values_from = Proportion)
write.csv(dat, path("results", "EMECO checks", "summary out", "Proportion OSPARIII Species Achieving GES.csv"),
          row.names = FALSE)

Plot yearly trend of percentage of species exceeding relative abundance threshold for OSPAR IV region

# OSPAR IV breeding -------------------------------------------------------------

# create empty vector to store yearly percentage of breeding species exceeding threshold
OSPARIV.prop.pass.breeding = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARIV)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding species 
  ### with six year relative abundance mean scores for the year
  OSPARIV.total = length(which(!is.na(unsm_dat_OSPARIV[, i])))
  
  # count the number of breeding species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding species exceeding 0.8 threshold for species laying 1 egg
  OSPARIV.pass = length(which(unsm_dat_OSPARIV[, i] > 0.7 & unsm_dat_OSPARIV$Eggs > 1)) +
    length(which(unsm_dat_OSPARIV[, i] > 0.8 & unsm_dat_OSPARIV$Eggs == 1))
  
  # add percentage of breeding species that exceed the threshold for that year
  OSPARIV.prop.pass.breeding[i] = OSPARIV.pass / OSPARIV.total * 100
}

dat_IV_breeding <- data.frame(Year = seq(startYearBreeding$IV + 5, endYearBreeding$IV, 1),
                    Group = "breeding",
                    Proportion = na.omit(OSPARIV.prop.pass.breeding))

# OSPAR IV breeding surface feeders ---------------------------------------------
unsm_dat_OSPARIV_surface <- unsm_dat_OSPARIV[unsm_dat_OSPARIV$Fun_Group == "Surface Feeder", ]

# create empty vector to store yearly percentage of breeding surface feeder species exceeding threshold
OSPARIV.prop.pass.surface = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARIV_surface)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding surface feeder species 
  ### with six year relative abundance mean scores for the year
  OSPARIV.total.surface = length(which(!is.na(unsm_dat_OSPARIV_surface[, i])))
  
  # count the number of breeding surface feeder species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding surface feeder species exceeding 0.8 threshold for species laying 1 egg
  OSPARIV.pass.surface = length(which(unsm_dat_OSPARIV_surface[, i] > 0.7 & unsm_dat_OSPARIV_surface$Eggs > 1)) +
    length(which(unsm_dat_OSPARIV_surface[, i] > 0.8 & unsm_dat_OSPARIV_surface$Eggs == 1))
  
  # add percentage of breeding surface feeder species that exceed the threshold for that year
  OSPARIV.prop.pass.surface[i] = OSPARIV.pass.surface / OSPARIV.total.surface * 100
}

dat_IV_surface <- data.frame(Year = seq(startYearBreeding$IV + 5, endYearBreeding$IV, 1),
                    Group = "breeding_surface_feeders",
                    Proportion = na.omit(OSPARIV.prop.pass.surface))

# OSPAR IV breeding water column feeders ----------------------------------------
unsm_dat_OSPARIV_water_column <- unsm_dat_OSPARIV[unsm_dat_OSPARIV$Fun_Group == "Water Column Feeder", ]

# create empty vector to store yearly percentage of breeding water column species exceeding threshold
OSPARIV.prop.pass.water.column = vector()

# iterate through each years six year relative abundance mean values
for (i in 10:ncol(unsm_dat_OSPARIV_water_column)) {  # skip first 4 columns and first 5 years
  
  # count total number of breeding water column species 
  ### with six year relative abundance mean scores for the year
  OSPARIV.total.water.column = length(which(!is.na(unsm_dat_OSPARIV_water_column[, i])))
  
  # count the number of breeding water column species exceeding 0.7 threshold for species laying more than 1 egg +
  # the number of breeding water column species exceeding 0.8 threshold for species laying 1 egg
  OSPARIV.pass.water.column = length(which(unsm_dat_OSPARIV_water_column[, i] > 0.7 & unsm_dat_OSPARIV_water_column$Eggs > 1)) +
    length(which(unsm_dat_OSPARIV_water_column[, i] > 0.8 & unsm_dat_OSPARIV_water_column$Eggs == 1))
  
  # add percentage of breeding water column feeder species that exceed the threshold for that year
  OSPARIV.prop.pass.water.column[i] = OSPARIV.pass.water.column / OSPARIV.total.water.column * 100
}

Plot

# set the x axis year range expanded to 5 year intervals
xAxisMinYear <- plyr::round_any(reduce(append(startYearBreeding$IV, startYearWintering$IV), min), 5)
xAxisMaxYear <- plyr::round_any(reduce(append(endYearBreeding$IV, endYearWintering$IV), max), 5, f = ceiling)

# create plot
plot <- ggplot() +
  theme(legend.position = c(0.25, 0.1),
        legend.key.width = unit(20, "mm")) +
  coord_cartesian(xlim = c(xAxisMinYear, xAxisMaxYear),
                  ylim = c(40, 100)) +
  scale_x_continuous(breaks = seq(xAxisMinYear, xAxisMaxYear, 5),
                     name = NULL) +
  scale_y_continuous(breaks = seq(40, 100, 10),
                     name = "Percentage of species with relative abundance above threshold value") +
  scale_colour_manual(name = NULL,
                      breaks = c("breeding", "surface"), 
                      values = c("#3b528b", "#3b528b"),
                      labels = c("all breeding", "breeding surface feeders")) + 
  scale_linetype_manual(name = NULL,
                      breaks = c("breeding", "surface"),
                      values = c("solid", "dashed"),
                      labels = c("all breeding", "breeding surface feeders")) + 
  # 75% cut off line
  geom_hline(yintercept = 75, linetype = "solid", colour = "black", size = 1.2) +
  geom_line(data = dat_IV_breeding,
            mapping = aes(Year, Proportion, linetype = "breeding", colour = "breeding"),  size = 1.5) +
  geom_line(data = dat_IV_surface,
            mapping = aes(Year, Proportion, linetype = "surface", colour = "surface"), size = 1.5) +
  annotate(geom = "text", x = xAxisMinYear, y = 74, 
           label = "Threshold value", color = "black", size = 3, hjust = 0.1)

# save plot
ggsave(plot = plot,
       filename = path("results", "summary out", "Trends_OSPARIV.png"),
       width = 2100, height = 2100, units = "px")

Export results

### output data to check with EMECO
dat <- bind_rows(dat_IV_breeding, dat_IV_surface) %>% # dat_IV_water_column, dat_IV_wintering, dat_IV_wader) %>%
  pivot_wider(names_from = Group, values_from = Proportion)
write.csv(dat, path("results", "EMECO checks", "summary out", "Proportion OSPARIV Species Achieving GES.csv"),
          row.names = FALSE)

Tidy up

# remove objects not required in subsequent analysis
required_data <- c("startYearBreeding", "endYearBreeding", "startYearWintering", "endYearWintering", 
                   "region_labels", "region_labels_breeding", "region_labels_wintering",
                   "functional_groups")

rm(list=ls()[!ls() %in% c("unsm_dat_OSPARI", "unsm_dat_OSPARIm", "unsm_dat_OSPARIn", "unsm_dat_OSPARIo",
                          "unsm_dat_OSPARII", "unsm_dat_OSPARIIa", "unsm_dat_OSPARIIb", "unsm_dat_OSPARIIc",
                          "unsm_dat_OSPARIId", "unsm_dat_OSPARIIe", "unsm_dat_OSPARIIf",
                          "unsm_dat_OSPARIII", 
                          "unsm_dat_OSPARIV",
                          "OSPARI.dat", "OSPARIn.dat", "OSPARIo.dat",
                          "OSPARII.dat", "OSPARIIa.dat", "OSPARIIb.dat", "OSPARIIc.dat",
                          "OSPARIId.dat", "OSPARIIe.dat", "OSPARIIf.dat",
                          "OSPARIII.dat",
                          required_data)])

Percentage of species within each functional group whose latest six year mean relative abundance exceeds threshold

Latest breeding mean relative abundance

get_latest_breeding_relative_abundance_means <- function(unsm_dat, end_year) {
  
  # end year column
  end_column <- str_glue("X{end_year}") 
  
  # restrict to latest relative abundance year
  unsm_dat <- unsm_dat %>%
    select(Species, spp, Eggs, Fun_Group, {{end_column}})
  names(unsm_dat) <- c("Species", "spp", "Eggs", "Fun_Group", "Selected_Year")
  
  # return latest relative abundance year
  return(unsm_dat)
  
}
unsm_dat_OSPARI_latest <- unsm_dat_OSPARI %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$I)
unsm_dat_OSPARIm_latest <- unsm_dat_OSPARIm %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$Im)
unsm_dat_OSPARIn_latest <- unsm_dat_OSPARIn %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$In)
unsm_dat_OSPARIo_latest <- unsm_dat_OSPARIo %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$Io)
unsm_dat_OSPARII_latest <- unsm_dat_OSPARII %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$II)
unsm_dat_OSPARIIa_latest <- unsm_dat_OSPARIIa %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$IIa)
unsm_dat_OSPARIIb_latest <- unsm_dat_OSPARIIb %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$IIb)
unsm_dat_OSPARIIc_latest <- unsm_dat_OSPARIIc %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$IIc)
unsm_dat_OSPARIId_latest <- unsm_dat_OSPARIId %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$IId)
unsm_dat_OSPARIIe_latest <- unsm_dat_OSPARIIe %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$IIe)
unsm_dat_OSPARIIf_latest <- unsm_dat_OSPARIIf %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$IIf)
unsm_dat_OSPARIII_latest <- unsm_dat_OSPARIII %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$III)
unsm_dat_OSPARIV_latest <- unsm_dat_OSPARIV %>%
  get_latest_breeding_relative_abundance_means(endYearBreeding$IV)

Latest non breeding mean relative abundance

get_latest_non_breeding_relative_abundance_means <- function(dat, end_year) {
  
  # end year column
  end_column <- as.character(end_year)
  
  # restrict to latest relative abundance year
  dat <- dat %>%
    select(Species, Common_name, Eggs, Fun_Group, {{end_column}})
  names(dat) <- c("Species", "Common_name", "Eggs", "Fun_Group", "Selected_Year")
  
  # return latest relative abundance year
  return(dat)
}
OSPARI.dat.latest <- OSPARI.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$I)
OSPARIn.dat.latest <- OSPARIn.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$In)
OSPARIo.dat.latest <- OSPARIo.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$Io)
OSPARII.dat.latest <- OSPARII.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$II)
OSPARIIa.dat.latest <- OSPARIIa.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$IIa)
OSPARIIb.dat.latest <- OSPARIIb.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$IIb)
OSPARIIc.dat.latest <- OSPARIIc.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$IIc)
OSPARIId.dat.latest <- OSPARIId.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$IId)
OSPARIIe.dat.latest <- OSPARIIe.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$IIe)
OSPARIIf.dat.latest <- OSPARIIf.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$IIf)
OSPARIII.dat.latest <- OSPARIII.dat %>%
  get_latest_non_breeding_relative_abundance_means(endYearWintering$III)

Proportion of species with latest six year mean above threshold for each functional group and OSPAR region and subdivision

# breeding
fgStats_breeding <- list()
for (i in region_labels_breeding){
  for (j in functional_groups) {
    
    # breeding species
    fgStats_breeding[[str_glue("OSPAR {i} breeding {j}")]] <-
      latest_proportion_breeding_species_above_threshold(
        latest.dat = get(str_glue("unsm_dat_OSPAR{i}_latest")),
        region = i,
        functional_group = j)
  }
}
fgStats_breeding <- bind_rows(fgStats_breeding)

# wintering
fgStats_wintering <- list()
for (i in region_labels_wintering){
  for (j in functional_groups) {
    
    # non-breeding wintering species
    fgStats_wintering[[str_glue("OSPAR {i} wintering {j}")]] <- 
      latest_proportion_wintering_species_above_threshold(
        latest.dat = get(str_glue("OSPAR{i}.dat.latest")),
        region = i,
        functional_group = j)
  }
}
fgStats_wintering <- bind_rows(fgStats_wintering)
fgStats <- bind_rows(fgStats_breeding, fgStats_wintering)
## Export results
write.csv(fgStats, 
          path("results", "EMECO checks", "summary out", "Functional Group Assessment Table.csv"), 
          row.names = FALSE,
          na = "")

Tidy up

rm(list=ls()[!ls() %in% c("fgStats")])