HerringLarvaeExplore.Rmd

---
title: "Explore Herring Larvae Indices"
author: "Sarah Gaichas"
date: "`r Sys.Date()`"
output:
  html_document:
    code_fold: hide
link-citations: yes
csl: "canadian-journal-of-fisheries-and-aquatic-sciences.csl"
bibliography: zoopindex.bib
urlcolor: blue
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      message = FALSE,
                      warning = FALSE)

library(tidyverse)
library(here)
library(DT)
library(ggiraph)
library(patchwork)
library(VAST)
library(FishStatsUtils)
```

# Introduction

Indices of various copeopod groups have been developed: https://noaa-edab.github.io/zooplanktonindex/CopeModResults.html

Now the question is, are the zooplankton available to herring larvae? We will explore data available on herring larvae in the EcoMon (and previous zooplankton) data.

# Methods

Herring larvae data were added to the input dataset in the updated script https://github.com/NOAA-EDAB/zooplanktonindex/blob/main/data/VASTzoopindex_processinputs.R and all stations were re-mapped to OISST data to fill missing temperature values if necessary. 

Where are herring larvae in each of our seasons?

```{r, fig.width=10}

herringfood_stn <- readRDS(here::here("data/herringfood_stn_all_OISST.rds"))

# make SST column that uses surftemp unless missing or 0

herringfood_stn <- herringfood_stn %>%
  dplyr::mutate(sstfill = ifelse((is.na(sfc_temp)|sfc_temp==0), oisst, sfc_temp),
                season_larv = month %in% c(1:2, 9:12))

herringlarvae_stn_fall <- herringfood_stn %>%
  #ungroup() %>%
  filter(season_ng == "FALL", 
         year > 1981) %>%
  mutate(AreaSwept_km2 = 1, #Elizabeth's code
         #declon = -declon already done before neamap merge
         Vessel = 1,
         Dayofyear = lubridate::yday(date) #as.numeric(as.factor(vessel))-1
  ) %>% 
  dplyr::select(Catch_g = cluhar_100m3, #use megabenwt for individuals input in example
                Year = year,
                Month = month,
                Dayofyear,
                Vessel,
                AreaSwept_km2,
                Lat = lat,
                Lon = lon,
                #btm_temp, #this leaves out many stations
                #sfc_temp, #this leaves out many stations
                #oisst,
                sstfill
  ) %>%
  na.omit() %>%
  as.data.frame()

herringlarvae_stn_sepfeb <- herringfood_stn %>%
  #ungroup() %>%
  dplyr::filter(season_larv == TRUE) %>%
  dplyr::mutate(AreaSwept_km2 = 1, #Elizabeth's code
         #declon = -declon already done before neamap merge
         Vessel = 1,
         Dayofyear = lubridate::yday(date),
         yearshift = ifelse(month < 3, year-1, year)#as.numeric(as.factor(vessel))-1
  ) %>% 
  dplyr::filter(yearshift>1981) |>
  dplyr::select(Catch_g = cluhar_100m3, #use megabenwt for individuals input in example
                Year = yearshift,
                Month = month,
                Dayofyear,
                Vessel,
                AreaSwept_km2,
                Lat = lat,
                Lon = lon,
                #btm_temp, #this leaves out many stations
                #sfc_temp, #this leaves out many stations
                #oisst,
                sstfill
  ) %>%
  na.omit() %>%
  as.data.frame()


herringlarvae_stn_spring <- herringfood_stn %>%
  #ungroup() %>%
  filter(season_ng == "SPRING", 
         year > 1981) %>%
  mutate(AreaSwept_km2 = 1, #Elizabeth's code
         #declon = -declon already done before neamap merge
         Vessel = 1,
         Dayofyear = lubridate::yday(date) #as.numeric(as.factor(vessel))-1
  ) %>% 
  dplyr::select(Catch_g = cluhar_100m3, #use megabenwt for individuals input in example
                Year = year,
                Month = month,
                Dayofyear,
                Vessel,
                AreaSwept_km2,
                Lat = lat,
                Lon = lon,
                #btm_temp, #this leaves out many stations
                #sfc_temp, #this leaves out many stations
                #oisst,
                sstfill
  ) %>%
  na.omit() %>%
  as.data.frame()


nonzerofall <- herringlarvae_stn_fall |>
  dplyr::filter(Catch_g > 0) #,
                #Year > 1981)

nonzerosepfeb <- herringlarvae_stn_sepfeb |>
  dplyr::filter(Catch_g > 0)

nonzerospring <- herringlarvae_stn_spring |>
  dplyr::filter(Catch_g > 0) #,
               # Year > 1981)

Fall <- ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat),  colour = "coral4", size=0.05, alpha=0.1) +
  geom_point(data = nonzerofall, aes(x = Lon, y = Lat), colour = "blue", size=0.5, alpha=1) +
  coord_sf(xlim =c(-78.5, -65.5), ylim = c(33, 45)) + #zoomed to Hatteras and N
  xlab("") +
  ylab("") +
  ggtitle("Fall herring larvae 1982-2022")+
  theme(plot.margin = margin(0, 0, 0, 0, "cm"))

SepFeb <- ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat),  colour = "coral4", size=0.05, alpha=0.1) +
  geom_point(data = nonzerosepfeb, aes(x = Lon, y = Lat), colour = "blue", size=0.5, alpha=1) +
  coord_sf(xlim =c(-78.5, -65.5), ylim = c(33, 45)) + #zoomed to Hatteras and N
  xlab("") +
  ylab("") +
  ggtitle("Sept-Feb herring larvae 1982-2022")+
  theme(plot.margin = margin(0, 0, 0, 0, "cm"))

Spring <- ggplot(data = ecodata::coast) +
  geom_sf() + 
  geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat),  colour = "coral4", size=0.05, alpha=0.1) +
  geom_point(data = nonzerospring, aes(x = Lon, y = Lat), colour = "blue", size=0.5, alpha=1) +
  coord_sf(xlim =c(-78.5, -65.5), ylim = c(33, 45)) + #zoomed to Hatteras and N
  xlab("") +
  ylab("") +
  ggtitle("Spring herring larvae 1982-2022")+
  theme(plot.margin = margin(0, 0, 0, 0, "cm"))
 
Spring + Fall + SepFeb
```

Day of year and months herring larvae found (present, not abundance)

```{r}

herringlarvae_stn_all <- dplyr::bind_rows(herringlarvae_stn_spring, herringlarvae_stn_fall)

hist(herringlarvae_stn_all$Dayofyear, xlim=c(0,365), breaks=366)

hist(herringlarvae_stn_all$Month, xlim=c(0,12), breaks=13)

```

Amount of herring larvae found (sum station volume over all years, not an abundance)

```{r}

sumlarvae <- herringlarvae_stn_all |>
  dplyr::group_by(Month) |>
  dplyr::summarise(totlarvae = sum(Catch_g, na.rm = TRUE))

ggplot2::ggplot(sumlarvae, ggplot2::aes(x=Month, y=totlarvae)) +
  ggplot2::geom_bar(stat = "identity")
```

So fall larvae, could be used to weight fall small copepods. 

What years? Also just summing stations, not an abundance estimate

```{r}

sumlarvaeyr <- herringlarvae_stn_all |>
  dplyr::group_by(Year) |>
  dplyr::summarise(totlarvae = sum(Catch_g, na.rm = TRUE))

ggplot2::ggplot(sumlarvaeyr, ggplot2::aes(x=Year, y=totlarvae)) +
  ggplot2::geom_bar(stat = "identity")
```

## Initial herring larvae model results

```{r}
# from each output folder in pyindex, 
outdir <- here::here("pyindex")
moddirs <- list.dirs(outdir) 
moddirs <- moddirs[-1]
# keep folder name
modnames <- list.dirs(outdir, full.names = FALSE)


# function to apply extracting info
getmodinfo <- function(d.name){
  # read settings
  modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
  modname <- modpath[length(modpath)]
  
  settings <- read.table(file.path(d.name, "settings.txt"), comment.char = "",
    fill = TRUE, header = FALSE)
  
  n_x <- as.numeric(as.character(settings[(which(settings[,1]=="$n_x")+1),2]))
  grid_size_km <- as.numeric(as.character(settings[(which(settings[,1]=="$grid_size_km")+1),2]))
  max_cells <- as.numeric(as.character(settings[(which(settings[,1]=="$max_cells")+1),2]))
  use_anisotropy <- as.character(settings[(which(settings[,1]=="$use_anisotropy")+1),2])
  fine_scale <- as.character(settings[(which(settings[,1]=="$fine_scale")+1),2])
  bias.correct <- as.character(settings[(which(settings[,1]=="$bias.correct")+1),2])
  
  #FieldConfig
  if(settings[(which(settings[,1]=="$FieldConfig")+1),1]=="Component_1"){
    omega1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+2),2])
    omega2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),1])
    epsilon1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),2])
    epsilon2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+5),1])
    beta1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+6),2])
    beta2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+7),1])
  }
  
  if(settings[(which(settings[,1]=="$FieldConfig")+1),1]=="Omega1"){
    omega1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),1])
    omega2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),1])
    epsilon1 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+3),2])
    epsilon2 <- as.character(settings[(which(settings[,1]=="$FieldConfig")+4),2])
    beta1 <- "IID"
    beta2 <- "IID"
  }
  
  
  #RhoConfig
  rho_beta1 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+3),1]))
  rho_beta2 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+3),2]))
  rho_epsilon1 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+4),1]))
  rho_epsilon2 <- as.numeric(as.character(settings[(which(settings[,1]=="$RhoConfig")+4),2]))
  
  # read parameter estimates, object is called parameter_Estimates
  if(file.exists(file.path(d.name, "parameter_estimates.RData"))) {
    load(file.path(d.name, "parameter_estimates.RData"))
    
    AIC <- parameter_estimates$AIC[1]  
    converged <- parameter_estimates$Convergence_check[1]
    fixedcoeff <- unname(parameter_estimates$number_of_coefficients[2])
    randomcoeff <- unname(parameter_estimates$number_of_coefficients[3])
    
  }else if(file.exists(file.path(d.name, "parameter_estimates.txt"))){
    
    reptext <- readLines(file.path(d.name, "parameter_estimates.txt"))
    
    AIC <- as.double(reptext[grep(reptext, pattern = "AIC")+1])
    converged <- reptext[grep(reptext, pattern = "Convergence_check")+1]
    fixedcoeff <- as.integer(stringr::str_split(reptext[grep(reptext, pattern = "number_of_coefficients")+2], 
                                     boundary("word"))[[1]][2])
    randomcoeff <- as.integer(stringr::str_split(reptext[grep(reptext, pattern = "number_of_coefficients")+2], 
                                     boundary("word"))[[1]][3])
    
  }else{
    
    AIC <- NA_real_
    converged <- NA_character_
    fixedcoeff <- NA_integer_
    randomcoeff <- NA_integer_
  }
  
  
  #index <- read.csv(file.path(d.name, "Index.csv"))
  
  
  # return model attributes as a dataframe
  out <- data.frame(modname = modname,
                    n_x = n_x,
                    grid_size_km = grid_size_km,
                    max_cells = max_cells,
                    use_anisotropy = use_anisotropy,
                    fine_scale =  fine_scale,
                    bias.correct = bias.correct,
                    omega1 = omega1,
                    omega2 = omega2,
                    epsilon1 = epsilon1,
                    epsilon2 = epsilon2,
                    beta1 = beta1,
                    beta2 = beta2,
                    rho_epsilon1 = rho_epsilon1,
                    rho_epsilon2 = rho_epsilon2,
                    rho_beta1 = rho_beta1,
                    rho_beta2 = rho_beta2,
                    AIC = AIC,
                    converged = converged,
                    fixedcoeff = fixedcoeff,
                    randomcoeff = randomcoeff#,
                    #index = index
  )
  	
	return(out)

}


modcompare <- purrr::map_dfr(moddirs, getmodinfo)

modselect <- modcompare |>
  dplyr::mutate(season = dplyr::case_when(stringr::str_detect(modname, "_fall_") ~ "Fall",
                            stringr::str_detect(modname, "spring") ~ "Spring",
                            stringr::str_detect(modname, "_all_") ~ "Annual",
                            TRUE ~ as.character(NA))) |>
  dplyr::mutate(converged2 = dplyr::case_when(stringr::str_detect(converged, "no evidence") ~ "likely",
                                stringr::str_detect(converged, "is likely not") ~ "unlikely",
                                TRUE ~ as.character(NA))) |>
  dplyr::mutate(copegroup = stringr::str_extract(modname, "[^_]+")) |>
  #dplyr::mutate(modname = str_extract(modname, '(?<=allagg_).*')) |>
  dplyr::group_by(copegroup, season) |>
  dplyr::mutate(deltaAIC = AIC-min(AIC)) |>
  dplyr::select(copegroup, modname, season, deltaAIC, fixedcoeff,
         randomcoeff, use_anisotropy, 
         omega1, omega2, epsilon1, epsilon2, 
         beta1, beta2, AIC, converged2) |>
  dplyr::arrange(copegroup, season, AIC)

# DT::datatable(modselect, rownames = FALSE, 
#               options= list(pageLength = 25, scrollX = TRUE),
#               caption = "Comparison of delta AIC values using Restricted Maxiumum Likelihood (REML) for alternative fixed and random effects model structures. See text for model descriptions.")

# flextable::flextable(modselect) %>%
#                        #dplyr::select(-c(use_anisotropy, 
#          #omega1, omega2, epsilon1, epsilon2, 
#          #beta1, beta2))
#          #) %>%
#   flextable::set_header_labels(modname = "Model name",
#                                season = "Season",
#                                #deltaAIC = "dAIC",
#                                fixedcoeff = "N fixed",
#                                randomcoeff = "N random",
#                                converged2 = "Convergence") |>
#   #flextable::set_caption("Comparison of delta AIC (dAIC) values using Restricted Maxiumum Likelihood (REML) for alternative fixed and random effects model structures, with number of fixed (N fixed) and random (N random) coefficients. See text for model descriptions associated with each model name.") %>%
#   flextable::fontsize(size = 9, part = "all") %>%
#   flextable::colformat_double(digits = 2) |>
#   flextable::set_table_properties(layout = "autofit", width = 1)


```

### Stations by season

Fall sampling for herring larvae was completed in most years aside from GLOBEC. In our definition of spring, herring larvae primarily occur in January and February. We now include a shifted season to better match larval herring availability throughout the time series: September - February. Year corresponds to September-December, and the following January and February are aligned with the previous year (hatch year) in these analyses. 

```{r, results = "asis"}
for(d.name in moddirs[str_detect(moddirs, "herring")]){
  
  modpath <- unlist(str_split(d.name, pattern = "/"))
  modname <- modpath[length(modpath)]
  
  cat(modname, "\n")
  cat(paste0("![](",d.name, "/Data_by_year.png)"))
  cat("\n\n")
  
}
```

### Indices by group, season, and region

```{r}
stratlook <- data.frame(Stratum = c("Stratum_1",
                                    "Stratum_2",
                                    "Stratum_3",
                                    "Stratum_4",
                                    "Stratum_5",
                                    "Stratum_6",
                                    "Stratum_7"),
                        Region  = c("AllEPU",
                                    "her_sp",
                                    "her_fa",
                                    "MAB",
                                    "GB",
                                    "GOM",
                                    "SS"))

# function to apply extracting info
getmodindex <- function(d.name){
  # read settings
  modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
  modname <- modpath[length(modpath)]
  
  if(file.exists(file.path(d.name,"Index.csv"))){
    index <- read.csv(file.path(d.name, "Index.csv"))
  }else{
    stopifnot()
  }
  # return model indices as a dataframe
  out <- data.frame(modname = modname,
                    index
  )
  
  return(out)
}

modcompareindex <- purrr::map_dfr(moddirs, purrr::possibly(getmodindex, otherwise = NULL))

splitoutput <- modcompareindex %>%
  dplyr::mutate(Season = modname |> map(str_split, pattern = "_") |> map_chr(c(1,2))) %>%
  dplyr::mutate(Data = modname |> map(str_split, pattern = "_") |> map_chr(c(1,1))) %>%
  dplyr::mutate(Estimate = ifelse(Estimate == 0, NA, Estimate)) |>
  dplyr::left_join(stratlook) #%>%
  #dplyr::filter(Region %in% c(GOM", "GB", "MAB","SS", "AllEPU")) use all regions

zoomax <- max(splitoutput$Estimate, na.rm=T)


zootsmean <- splitoutput %>%
  dplyr::group_by(modname, Region) %>%
  dplyr::mutate(fmean = mean(Estimate, na.rm=T)) 


```

### Seasonal indices

```{r, fig.height=10}

plot_zooindices <- function(splitoutput, plotdata, plotregions, plottitle){
  
  filterEPUs <- plotregions #c("her_sp", "her_fa", "MAB", "GB", "GOM", "SS", "AllEPU")
  
  seasons <- splitoutput |> dplyr::filter(Data==plotdata) |> dplyr::select(Season) |> dplyr::distinct()
  
  ncols <- dim(seasons)[1]
  
  currentMonth <- lubridate::month(Sys.Date())
  currentYear <- lubridate::year(Sys.Date())
  if (currentMonth > 4) {
    endShade <- currentYear
  } else {
    endShade <- currentYear - 1
  }
  shadedRegion <- c(endShade-9,endShade)
  
  shade.alpha <- 0.3
  shade.fill <- "lightgrey"
  lwd <- 1
  pcex <- 2
  trend.alpha <- 0.5
  trend.size <- 2
  hline.size <- 1
  line.size <- 2
  hline.alpha <- 0.35
  hline.lty <- "dashed"
  label.size <- 5
  hjust.label <- 1.5
  letter_size <- 4
  errorbar.width <- 0.25
  x.shade.min <- shadedRegion[1]
  x.shade.max <- shadedRegion[2]
  
  setup <- list(
    shade.alpha = shade.alpha,
    shade.fill =shade.fill,
    lwd = lwd,
    pcex = pcex,
    trend.alpha = trend.alpha,
    trend.size = trend.size,
    line.size = line.size,
    hline.size = hline.size,
    hline.alpha = hline.alpha,
    hline.lty = hline.lty,
    errorbar.width = errorbar.width,
    label.size = label.size,
    hjust.label = hjust.label,
    letter_size = letter_size,
    x.shade.min = x.shade.min,
    x.shade.max = x.shade.max
  )
  
  
  fix<- splitoutput |>
    dplyr::filter(Data %in% plotdata, #c("calfin"),
                  Region %in% filterEPUs) |>
    dplyr::group_by(Region, Season) |>
    dplyr::summarise(max = max(Estimate, na.rm=T))
  
  p <- splitoutput |>
    dplyr::filter(Data %in% plotdata, #c("calfin"),
                  Region %in% filterEPUs) |>
    dplyr::group_by(Region, Season) |>
    dplyr::left_join(fix) |>
    dplyr::mutate(#Value = Value/resca,
      Mean = as.numeric(Estimate),
      SE = Std..Error.for.Estimate,
      Mean = Mean/max,
      SE = SE/max,
      Upper = Mean + SE,
      Lower = Mean - SE) |>
    ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, linetype = modname, group = modname))+
    ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
                      xmin = setup$x.shade.min , xmax = setup$x.shade.max,
                      ymin = -Inf, ymax = Inf) +
    ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Season), alpha = 0.5)+
    ggplot2::geom_point()+
    ggplot2::geom_line()+
    ggplot2::ggtitle(plottitle)+
    ggplot2::ylab(expression("Relative abundance"))+
    ggplot2::xlab(ggplot2::element_blank())+
    ggplot2::facet_wrap(Region~Season, ncol = ncols, 
                        labeller = label_wrap_gen(multi_line=FALSE))+
    ecodata::geom_gls()+
    ecodata::theme_ts()+
    ecodata::theme_facet()+
    ecodata::theme_title() +
    ggplot2::theme(legend.position = "bottom")
  
  return(p)
}


plot_zooindices(splitoutput = splitoutput, 
             plotdata = "herringlarvae", 
             plotregions = c("her_sp", "her_fa", "MAB", "GB", "GOM", "SS", "AllEPU"), 
             plottitle = "Herring larvae") 
```

Relative density by area.

We now see the spike in 2000 that was observed by @richardson_development_2010. 

```{r}

plotdata <- c("herringlarvae")
plottitle  <-  "Herring larvae"

  fix<- splitoutput |>
    dplyr::filter(Data %in% plotdata #, 
                  #Region %in% filterEPUs
                  ) |>
    dplyr::group_by(Season) |> #Region,
    dplyr::summarise(max = max(Estimate, na.rm=T))
  
  p <- splitoutput |>
    dplyr::filter(Data %in% plotdata #, #c("calfin"),
                  #Region %in% filterEPUs
                  ) |>
    dplyr::group_by(Season) |> #Region, 
    dplyr::left_join(fix) |>
    dplyr::mutate(#Value = Value/resca,
      Mean = as.numeric(Estimate),
      SE = Std..Error.for.Estimate,
      Mean = Mean/max,
      SE = SE/max,
      Upper = Mean + SE,
      Lower = Mean - SE) |>
    ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, linetype = Region, group = Region))+
    #ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
    #                  xmin = setup$x.shade.min , xmax = setup$x.shade.max,
    #                  ymin = -Inf, ymax = Inf) +
    ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Region), alpha = 0.5)+
    ggplot2::geom_point()+
    ggplot2::geom_line()+
    ggplot2::ggtitle(plottitle)+
    ggplot2::ylab(expression("Relative abundance"))+
    ggplot2::xlab(ggplot2::element_blank())+
    ggplot2::facet_wrap(~Season, #Region~ ncol = ncols, 
                        labeller = label_wrap_gen(multi_line=FALSE))+
    #ecodata::geom_gls()+
    ecodata::theme_ts()+
    ecodata::theme_facet()+
    ecodata::theme_title() +
    ggplot2::theme(legend.position = "bottom")
  
 p

```


### Density estimates

```{r, results = "asis"}
for(d.name in moddirs[str_detect(moddirs, "herring")]){
  
  modpath <- unlist(str_split(d.name, pattern = "/"))
  modname <- modpath[length(modpath)]
  
  cat(modname, "\n")
  if(file.exists(paste0(d.name, "/ln_density-predicted.png"))){
    cat(paste0("![](",d.name, "/ln_density-predicted.png)"))
  }
  cat("\n\n")
  
}
```

## Extract the herring larvae data by year

We want to define areas of most dense larvae each year and pull our small copepod index from there.

Maybe quantiles of herring larval density by year?

Plot the data (based on https://github.com/James-Thorson-NOAA/VAST/wiki/Plots-using-ggplot):

Two low years and two high years. Is distribution different?

```{r, fig.height=12}

d.name <- moddirs[str_detect(moddirs, "herringlarvae_sepfeb")]
fit <- readRDS(paste0(d.name, "/fit.rds"))
#fit <- VAST::reload_model(fit) #added to try to make work after restart, no previous VAST run
years <- unique(fit$data_frame$t_i)
years <- c(min(years):max(years))

mdl <- FishStatsUtils::make_map_info(Region = fit$settings$Region,
                     spatial_list = fit$spatial_list,
                     Extrapolation_List = fit$extrapolation_list#,
                     #added to try to make work after restart, no previous VAST run
                     #Include = fit$extrapolation_list[["Area_km2_x"]] > 0 &  rowSums(fit$extrapolation_list[["a_el"]]) >  0
                     )

gmap <- ggplot(data = ecodata::coast) +
  geom_sf() + 
  #aes(x = lon, y = lat, group = group) +
  #geom_polygon(fill="black", colour = "white") +
  scale_color_viridis_c(option = "magma") +  # now make this quantiles...
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.spacing.x=unit(0, "lines"),
        panel.spacing.y=unit(0, "lines"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank() ) +
  coord_sf(xlim=mdl$Xlim, ylim=mdl$Ylim)

## Below shows to you get the model estimate of density, D_gct,
## for each grid (g), category (c; not used here single
## univariate); and year (t); and link it spatially to a lat/lon
## extrapolation point.  You can do this for any _gct or _gc
## variable in the Report.
names(fit$Report)[grepl('_gc|_gct', x=names(fit$Report))]

D_gt <- fit$Report$D_gct[,1,] # drop the category
dimnames(D_gt) <- list(cell=1:nrow(D_gt), year=years)
## tidy way of doing this, reshape2::melt() does
## it cleanly but is deprecated
D_gt <- D_gt %>% as.data.frame() %>%
    tibble::rownames_to_column(var = "cell") %>%
    pivot_longer(-cell, names_to = "Year", values_to='D')
D <- merge(D_gt, mdl$PlotDF, by.x='cell', by.y='x2i')

saveRDS(D, "Dherr_sepfeb.rds")

g <- gmap +
  geom_point(data=D, aes(Lon, Lat, color=log(as.vector(D)), group=NULL),
             ## These settings are necessary to avoid
             ## overlplotting which is a problem here. May need
             ## to be tweaked further.
             size=.3, stroke=0,shape=16) + facet_wrap('Year')
#g


highyears <- c(1992, 2000)
lowyears <- c(1983, 2019)

Dsub <- D |> dplyr::filter(Year %in% c(lowyears, highyears))

g <- gmap +
  geom_point(data=Dsub, aes(Lon, Lat, color=log(as.vector(D)), group=NULL)
             #,
             ## These settings are necessary to avoid
             ## overlplotting which is a problem here. May need
             ## to be tweaked further.
             #size=.3, stroke=0,shape=16
             ) + facet_wrap('Year')

g

```

Calculate quantiles of distribution across the time series. First sum all cells over time, then these are the quantiles of summed log density:

```{r, fig.width=10}

# one option, sum D_gct over all years then take quantiles in space

D <- readRDS("Dherr_sepfeb.rds")

Dtot <- D |> 
  dplyr::group_by(cell) |>
  dplyr::mutate(Dsum = sum(D, na.rm=TRUE),
                logD = log(as.vector(Dsum))) |>
  dplyr::select(!c(D, Year)) |>
  dplyr::distinct()

Dvec <- terra::vect(Dtot, geom=c("Lon", "Lat"))

q <- quantile(log(as.vector(Dtot$Dsum)), probs=seq(0,1,0.1))
q2 <- quantile(Dtot$logD, probs=seq(0,1,0.1))

qvec <- terra::quantile(Dvec, probs=c(0.6, 0.65, 0.7, 0.75, 0.8))

quants <- classInt::classIntervals(log(as.vector(Dtot$Dsum)), 
                                   style = "quantile", n = 10)

#quants$brks


D60pct <- Dtot |>
  dplyr::filter(log(as.vector(Dsum))>qvec["60%","logD"]) |>
  dplyr::rename("60th" = "Include")

D65pct <- Dtot |>
  dplyr::filter(log(as.vector(Dsum))>qvec["65%","logD"]) |>
  dplyr::rename("65th" = "Include")

D70pct <- Dtot |>
  dplyr::filter(log(as.vector(Dsum))>qvec["70%","logD"]) |>
  dplyr::rename("70th" = "Include")

D75pct <- Dtot |>
  dplyr::filter(log(as.vector(Dsum))>qvec["75%","logD"]) |>
  dplyr::rename("75th" = "Include")

D80pct <- Dtot |>
  dplyr::filter(log(as.vector(Dsum))>qvec["80%","logD"]) |>
  dplyr::rename("80th" = "Include")

qvec

```


Mapping the summed density below.

```{r}
g <- gmap +
  geom_point(data=Dtot, aes(Lon, Lat, color=log(as.vector(Dsum)), group=NULL)
             ,
             ## These settings are necessary to avoid
             ## overlplotting which is a problem here. May need
             ## to be tweaked further.
             size=1.5, stroke=0,shape=16
             ) 

g
```


Mapping the quantiles from 60% (white) to 80% (light blue).

Quantiles below 70% include densities south of Long Island and off Cape Hatteras. The WG thinks the larvae found to the south may be lost to the population so we likely don't want to use a quantile below the 70th percentile.

The 80th percentile starts to show a gap on Georges Bank and along the southwest coast of Maine. This may be slicing things too finely to cover general herring larval habitat. 

I think this leaves us with using either the 70th (light green) or 75th (dark green) percentile of summed density over all years to define herring larvae relevant habitat for the small copepods index.

```{r}
 
gq <- gmap +
 geom_point(data=D60pct, 
             aes(x=Lon, y=Lat, color = "60th"), 
            # color="white"
            #)
            # ,
             ## These settings are necessary to avoid
             ## overlplotting which is a problem here. May need
             ## to be tweaked further.
             size=1.5, stroke=0,shape=16
             )  +
  
  geom_point(data=D65pct, 
             aes(x=Lon, y=Lat, color = "65th"), 
             #color="yellow"
             #)
             #,
             ## These settings are necessary to avoid
             ## overlplotting which is a problem here. May need
             ## to be tweaked further.
             size=1.5, stroke=0,shape=16
             )  +
  
  geom_point(data=D70pct, 
             aes(x=Lon, y=Lat, color = "70th"), 
             #color="green"
             #)
             #,
             ## These settings are necessary to avoid
             ## overlplotting which is a problem here. May need
             ## to be tweaked further.
             size=1.5, stroke=0,shape=16
             )  +
  
  geom_point(data=D75pct, 
             aes(x=Lon, y=Lat, color = "75th"), 
             #color="darkgreen"
              #)
             #,
             ## These settings are necessary to avoid
             ## overlplotting which is a problem here. May need
             ## to be tweaked further.
             size=1.5, stroke=0,shape=16
             )  +
  
   geom_point(data=D80pct, 
             aes(x=Lon, y=Lat, color = "80th"), 
             #color="lightblue"
             #)
             #,
             ## These settings are necessary to avoid
             ## overlplotting which is a problem here. May need
             ## to be tweaked further.
             size=1.5, stroke=0,shape=16
             )  +
  
  ggplot2::scale_color_manual(name = "Quantile",
                        breaks = c("60th", "65th", "70th", "75th", "80th"),
                        values = c("60th" = "white", "65th" = "yellow", 
                                   "70th" = "green", "75th" = "darkgreen", 
                                   "80th" = "lightblue") )
  
gq


# another option, pull from fit$report Omega 1 and 2?

```


## Make a new extrapolation grid using herring larvae density 70th percentile

These methods are similar to those used in the bluefish RTA for the forage index nearshore strata 3 miles from shore.

First make the 70th percentile+ points into an `sf` object, then intersect that object with the built in `FishStatsUtils::northwest_atlantic_grid`:

```{r, eval=FALSE}

# methods from https://stackoverflow.com/questions/78335772/find-outer-edge-of-polygon-in-r
# after much trial and error, concave_hull is what we want, only available in newer sf

# dataframe to sf object
D70pct_sf <- sf::st_as_sf(D70pct, coords = c("Lon", "Lat")) 
  
# concave hull in newest sf only works with GEOS>3.11
D70pct_ls <- D70pct_sf |>
  sf::st_union() |>
  sf::st_concave_hull(ratio=0.1) |>
  sf::st_cast(to ="LINESTRING") |>
  sf::st_cast(to ="POLYGON") |>
  sf::st_set_crs(sf::st_crs(ecodata::coast))

# just in case 75th too
D75pct_sf <- sf::st_as_sf(D75pct, coords = c("Lon", "Lat")) 

D75pct_ls <- D75pct_sf |>
  sf::st_union() |>
  sf::st_concave_hull(ratio=0.1) |>
  sf::st_cast(to ="LINESTRING") |>
  sf::st_cast(to ="POLYGON") |>
  sf::st_set_crs(sf::st_crs(ecodata::coast))

# Dont need this? set crs from ecodata::coast
# # set bounding boxes
# neus.xmin=-77
# neus.xmax=-65
# neus.ymin=35
# neus.ymax=45
# 
# neus.bbox1 <- sf::st_set_crs(sf::st_as_sf(as(raster::extent(neus.xmin, neus.xmax, neus.ymin, neus.ymax), "SpatialPolygons")), sf::st_crs(ecodata::coast))
# 
# neus.bbox2 <- sf::st_set_crs(sf::st_as_sf(as(raster::extent(-78, -74, 42, 45), "SpatialPolygons")), sf::st_crs(ecodata::coast)) # smaller bounding box to get rid of extra lines on the map 
# 
# neuscoast <- ecodata::coast |> 
#   sf::st_intersection(neus.bbox1) |>
#   sf::st_difference(neus.bbox2) # gets rid of extra non coastal line 

# intersect buffer with the current FishStatsUtils::northwest_atlantic_grid
# make northwest atlantic grid into sf object
nwagrid_sf  <-  sf::st_as_sf(FishStatsUtils::northwest_atlantic_grid, coords = c("Lon","Lat")) %>%
  sf::st_set_crs(sf::st_crs(ecodata::coast))

# intersect, rearrange in same format as nwatl grid, and save
D70pct_nwa <- sf::st_intersection(nwagrid_sf,D70pct_ls) %>% #native pipe wont do dots
  dplyr::mutate(Lon = as.numeric(sf::st_coordinates(.)[,1]),
                Lat = as.numeric(sf::st_coordinates(.)[,2])) |>
  sf::st_set_geometry(NULL) |>
  #dplyr::select(-featurecla) |>
  dplyr::select(stratum_number, Lon, Lat, everything())

write_rds(D70pct_nwa, here("spatialdat","D70pct_nwa.rds"))

# intersect, rearrange in same format as nwatl grid, and save
D75pct_nwa <- sf::st_intersection(nwagrid_sf,D75pct_ls) %>% #native pipe wont do dots
  dplyr::mutate(Lon = as.numeric(sf::st_coordinates(.)[,1]),
                Lat = as.numeric(sf::st_coordinates(.)[,2])) |>
  sf::st_set_geometry(NULL) |>
  #dplyr::select(-featurecla) |>
  dplyr::select(stratum_number, Lon, Lat, everything())

write_rds(D75pct_nwa, here("spatialdat","D75pct_nwa.rds"))


```

The portions of `nortwest_atlantic_grid` intersecting with the 70th and 75th percentile of herring larval density (1982-2022) were saved in the spatialdat folder. Next, we define new strata based on that intersection and make a new extrapolation grid. Then we can use this grid and call the new strata as strata.limits when running the small copepods model.

Right now just make a grid for the 70th percentile; we can make one for 75th if needed later.

```{r, eval=FALSE}

D70pct_nwa <- readRDS(here("spatialdat/D70pct_nwa.rds"))

D70pct_nwast <- D70pct_nwa %>%
  dplyr::mutate(strat2 = 1) %>% #herring larvae = 1
  dplyr::right_join(FishStatsUtils::northwest_atlantic_grid) %>%
  dplyr::mutate(strat2 = replace_na(strat2, 2)) %>% #replace NA with 2 for non-larval
  dplyr::mutate(stratum_number2 = as.numeric(paste0(stratum_number, strat2))) %>%
  dplyr::select(-strat2)

saveRDS(D70pct_nwast, here("spatialdat","D70pct_nwa_strat2.rds"))

# new lookups

herring_spring <- c(01010, 01020, 01030, 01040, 01050, 01060, 01070, 01080, 01090, 
                    01100, 01110, 01120, 01130, 01140, 01150, 01160, 01170, 01180, 
                    01190, 01200, 01210, 01220, 01230, 01240, 01250, 01260, 01270, 
                    01280, 01290, 01300, 01360, 01370, 01380, 01390, 01400, 01610, 
                    01620, 01630, 01640, 01650, 01660, 01670, 01680, 01690, 01700, 
                    01710, 01720, 01730, 01740, 01750, 01760)
herring_fall <- c(01050, 01060, 01070, 01080, 01090, 01100, 01110, 01120, 01130, 
                  01140, 01150, 01160, 01170, 01180, 01190, 01200, 01210, 01220, 
                  01230, 01240, 01250, 01260, 01270, 01280, 01290, 01300, 01360, 
                  01370, 01380, 01390, 01400)

MAB <- c(1010:1080, 1100:1120, 1600:1750, 3010:3450, 3470, 3500, 3510)
GB  <- c(1090, 1130:1210, 1230, 1250, 3460, 3480, 3490, 3520:3550)
GOM <- c(1220, 1240, 1260:1290, 1360:1400, 3560:3830)
SS  <- c(1300:1352, 3840:3990)

# MAB EPU
MAB2 <- D70pct_nwast %>% 
  dplyr::filter(stratum_number %in% MAB) %>%
  dplyr::select(stratum_number2) %>%
  dplyr::distinct()

# MAB herring larvae area
MAB2herr <- MAB2 %>%
  dplyr::filter(stratum_number2 %% 10 == 1) 

# MAB outside larval area
MAB2out <- MAB2 %>%
  dplyr::filter(stratum_number2 %% 10 == 2) 

# Georges Bank EPU
GB2 <- D70pct_nwast %>% 
  dplyr::filter(stratum_number %in% GB) %>%
  dplyr::select(stratum_number2) %>%
  dplyr::distinct()

# GB herring larvae
GB2herr <- GB2 %>%
  dplyr::filter(stratum_number2 %% 10 == 1) 

#GB outside larval area
GB2out <- GB2 %>%
  dplyr::filter(stratum_number2 %% 10 == 2)

# gulf of maine EPU 
GOM2 <- D70pct_nwast %>%
  dplyr::filter(stratum_number %in% GOM) %>%
  dplyr::select(stratum_number2) %>%
  dplyr::distinct()

# GOM herring larvae
GOM2herr <- GOM2 %>%
  dplyr::filter(stratum_number2 %% 10 == 1) 

#GOM outside larval area
GOM2out <- GOM2 %>%
  dplyr::filter(stratum_number2 %% 10 == 2)

# scotian shelf EPU 
SS2 <- D70pct_nwast %>%
  dplyr::filter(stratum_number %in% SS) %>%
  dplyr::select(stratum_number2) %>%
  dplyr::distinct()

# SS herring larvae
SS2herr <- SS2 %>%
  dplyr::filter(stratum_number2 %% 10 == 1) 

#SS outside larval area
SS2out <- SS2 %>%
  dplyr::filter(stratum_number2 %% 10 == 2)

# whole herring larval area
herrlarv <- dplyr::bind_rows(MAB2herr, GB2herr, GOM2herr, SS2herr)

# outside herring larval area
nolarv <- dplyr::bind_rows(MAB2out, GB2out, GOM2out, SS2out)

# spring herring NEFSC BTS
her_spr2 <- D70pct_nwast %>%
  dplyr::filter(stratum_number %in% herring_spring) %>%
  dplyr::select(stratum_number2) %>%
  dplyr::distinct()

# fall herring NEFSC BTS
her_fall2 <- D70pct_nwast %>%
  dplyr::filter(stratum_number %in% herring_fall) %>%
  dplyr::select(stratum_number2) %>%
  dplyr::distinct()

# all epus
allEPU2 <- D70pct_nwast %>%
  dplyr::filter(stratum_number %in% c(MAB, GB, GOM, SS)) %>%
  dplyr::select(stratum_number2) %>%
  dplyr::distinct()

```


Modify the function from `FishStatsUtils::Prepare_NWA_Extrapolation_Data_Fn` to make a new grid with updated strata

```{r newprepextrap, eval=FALSE}

Prepare_NWA_Extrapolation_Data_Fn_skg <- function (strata.limits = NULL, 
                                               epu_to_use = c("All", "Georges_Bank", "Mid_Atlantic_Bight", "Scotian_Shelf", "Gulf_of_Maine", "Other")[1],
                                               projargs = NA, zone = NA, flip_around_dateline = FALSE, ...) 
{
    if (is.null(strata.limits)) {
        strata.limits = list(All_areas = 1:1e+05)
    }
    message("Using strata ", strata.limits)
    if (any(tolower(epu_to_use) %in% "all")) {
        epu_to_use <- c("Georges_Bank", "Mid_Atlantic_Bight", 
            "Scotian_Shelf", "Gulf_of_Maine", "Other")
    }
    utils::data(northwest_atlantic_grid, package = "FishStatsUtils")
    Data_Extrap <- D70pct_nwast
    Tmp = cbind(BEST_DEPTH_M = 0, BEST_LAT_DD = Data_Extrap[, 
        "Lat"], BEST_LON_DD = Data_Extrap[, "Lon"])
    if (length(strata.limits) == 1 && strata.limits[1] == "EPU") {
        Data_Extrap <- Data_Extrap[Data_Extrap$EPU %in% epu_to_use, 
            ]
        Data_Extrap$EPU <- droplevels(Data_Extrap$EPU)
        a_el = matrix(NA, nrow = nrow(Data_Extrap), ncol = length(epu_to_use), 
            dimnames = list(NULL, epu_to_use))
        Area_km2_x = Data_Extrap[, "Area_in_survey_km2"]
        for (l in 1:ncol(a_el)) {
            a_el[, l] = ifelse(Data_Extrap[, "EPU"] %in% epu_to_use[[l]], 
                Area_km2_x, 0)
        }
    }
    else {
        a_el = as.data.frame(matrix(NA, nrow = nrow(Data_Extrap), 
            ncol = length(strata.limits), dimnames = list(NULL, 
                names(strata.limits))))
        Area_km2_x = Data_Extrap[, "Area_in_survey_km2"]
        for (l in 1:ncol(a_el)) {
            a_el[, l] = ifelse(Data_Extrap[, "stratum_number2"] %in% 
                strata.limits[[l]], Area_km2_x, 0)
        }
    }
    tmpUTM = project_coordinates(X = Data_Extrap[, "Lon"], Y = Data_Extrap[, 
        "Lat"], projargs = projargs, zone = zone, flip_around_dateline = flip_around_dateline)
    Data_Extrap = cbind(Data_Extrap, Include = 1)
    Data_Extrap[, c("E_km", "N_km")] = tmpUTM[, c("X", "Y")]
    Return = list(a_el = a_el, Data_Extrap = Data_Extrap, zone = attr(tmpUTM, 
        "zone"), projargs = attr(tmpUTM, "projargs"), flip_around_dateline = flip_around_dateline, 
        Area_km2_x = Area_km2_x)
    return(Return)
}


```

Now define new `strata.limits`. *Needed? We do this at runtime*

```{r stratalimits, eval=FALSE}

strata.limits <- as.list(c("AllEPU" = allEPU2,
                           "her_sp" = her_spr2,
                           "her_fa" = her_fall2,
                           "her_larv" = herrlarv,
                           "no_larv" = nolarv,
                           "MAB" = MAB2,
                           "GB" = GB2,
                           "GOM" = GOM2,
                           "SS" = SS2
))
```

Make the new extrapolation list:

```{r extraplist, eval=FALSE}

Extrapolation_List  <-  Prepare_NWA_Extrapolation_Data_Fn_skg( strata.limits=strata.limits)

saveRDS(Extrapolation_List, file = here("spatialdat/CustomExtrapolationList.rds"))

```

Did it work? Plot

```{r}

newstrat <- readRDS(here("spatialdat/D70pct_nwa_strat2.rds"))

herrlarv_area <- newstrat |>
  dplyr::filter(stratum_number2 %% 10 == 1)

ggplot2::ggplot(data = ecodata::coast) +
  ggplot2::geom_sf() + 
  ggplot2::geom_point(data = FishStatsUtils::northwest_atlantic_grid, ggplot2::aes(x = Lon, y = Lat), size=0.05, colour = "brown", alpha=0.1) +
  ggplot2::geom_point(data = herrlarv_area, ggplot2::aes(x = Lon, y = Lat), size=0.05, colour = "green",  alpha=0.5) +
  ggplot2::coord_sf(xlim = c(-78, -65.5), ylim = c(35, 45))+
  ggplot2::ggtitle("Herring larvae area: northwest_atlantic_grid")

```

# Results: Small copepods in the herring larval area

Using the new grid, run the smallcopeALL_sepfeb dataset in VAST to get an index in the herring larvae area as well as in the others. 

Script for doing this is https://github.com/NOAA-EDAB/zooplanktonindex/blob/main/VASTscripts/VASTunivariate_zoopindex_smcopeALL_herrlarvarea.R

The model converged, diagnostics look fine. 

Visualize these small copepod model indices:

```{r, fig.height=10, fig.asp=2}

# add the new herring larval stratum

stratlook2 <- data.frame(Stratum = c("Stratum_1",
                                    "Stratum_2",
                                    "Stratum_3",
                                    "Stratum_4",
                                    "Stratum_5",
                                    "Stratum_6",
                                    "Stratum_7",
                                    "Stratum_8",
                                    "Stratum_9"),
                        Region  = c("AllEPU",
                                    "her_sp",
                                    "her_fa",
                                    "her_larv",
                                    "no_larv",
                                    "MAB",
                                    "GB",
                                    "GOM",
                                    "SS"))

# only larvarea models have this strata set
splitoutput2 <- modcompareindex %>%
  dplyr::filter(str_detect(modname, "larvarea")) |>
  dplyr::mutate(Season = modname |> map(str_split, pattern = "_") |> map_chr(c(1,2))) %>%
  dplyr::mutate(Data = modname |> map(str_split, pattern = "_") |> map_chr(c(1,1))) %>%
  dplyr::mutate(Estimate = ifelse(Estimate == 0, NA, Estimate)) |>
  dplyr::left_join(stratlook2)

# tack on what should be the same model in the original extrapolation grid to make sure nothing went wrong
# the same strata should have the same trends
smcopesepfeb <- splitoutput |>
  dplyr::filter(modname == "smallcopeALL_sepfeb_yrshift_500_biascorrect")

splitoutput2 <- dplyr::bind_rows(splitoutput2, smcopesepfeb)

plot_zooindices(splitoutput = splitoutput2, 
             plotdata = "smallcopeALL", 
             plotregions = c("her_larv", "no_larv", "her_sp", "her_fa", "MAB", "GB", "GOM", "SS", "AllEPU"), 
             plottitle = "Small copepods Sept-Feb") 


```

We get the same time series for the comparable areas using the original and new extrapolation grid, as expected. (All time series overlap.)

Compare small copepods inside and outside herring larval area:

```{r}
plot_zooindices(splitoutput = splitoutput2, 
             plotdata = "smallcopeALL", 
             plotregions = c("her_larv", "no_larv"), 
             plottitle = "Small copepods Sept-Feb") 
```

On the same plot to compare scale:

```{r}
plotdata <- c("smallcopeALL")
plottitle  <-  "Small copepods Sept-Feb inside (herr_larv) and outside (no_larv) herring larval area"

  fix<- splitoutput2 |>
    dplyr::filter(Data %in% plotdata, 
                  Region %in% c("her_larv", "no_larv")
                  ) |>
    dplyr::group_by(Season) |> #Region,
    dplyr::summarise(max = max(Estimate, na.rm=T))
  
  p <- splitoutput2 |>
    dplyr::filter(Data %in% plotdata, #c("calfin"),
                  Region %in% c("her_larv", "no_larv")
                  ) |>
    dplyr::group_by(Season) |> #Region, 
    dplyr::left_join(fix) |>
    dplyr::mutate(#Value = Value/resca,
      Mean = as.numeric(Estimate),
      SE = Std..Error.for.Estimate,
      Mean = Mean/max,
      SE = SE/max,
      Upper = Mean + SE,
      Lower = Mean - SE) |>
    ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, linetype = Region, group = Region))+
    #ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
    #                  xmin = setup$x.shade.min , xmax = setup$x.shade.max,
    #                  ymin = -Inf, ymax = Inf) +
    ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Region), alpha = 0.5)+
    ggplot2::geom_point()+
    ggplot2::geom_line()+
    ggplot2::ggtitle(plottitle)+
    ggplot2::ylab(expression("Relative abundance"))+
    ggplot2::xlab(ggplot2::element_blank())+
    ggplot2::facet_wrap(~Season, #Region~ ncol = ncols, 
                        labeller = label_wrap_gen(multi_line=FALSE))+
    #ecodata::geom_gls()+
    ecodata::theme_ts()+
    ecodata::theme_facet()+
    ecodata::theme_title() +
    ggplot2::theme(legend.position = "bottom")
  
 p
```

The small copepods inside the herring larval area have been about the same or more than those outside recently. 

We'll see if any of this helps as a covariate on recruitment...

# References