Testing atlantisom: generate survey sampled length comps for Norwegian Barents Sea Atlantis

Introduction

This page has visualizations for the NOBA model example, CERES Global Sustainability. For full explanation of methods, see the file linked at the beginning of each section.

Simulate a survey part 4: sample for length composition

Full methods are explained here.

To create a survey, the user specifies the timing of the survey, which species are captured, the spatial coverage of the survey, the species-specific survey efficiency (“q”), and the selectivity at age for each species.

# make defaults that return a standard survey, implement in standard_survey
# users need to map their species groups into these general ones
#   large pelagics/reef associated/burrowers/otherwise non-trawlable
#   pelagics
#   demersals
#   selected flatfish
if(initNOBA) funct.groups <- rename(funct.groups, GroupType = InvertType)

survspp <- funct.groups$Name[funct.groups$IsTurnedOn==1 &
                           funct.groups$GroupType %in% c("FISH", "SHARK")]

if(initCCA) { #Sarah's CCA Grouping
  nontrawl <- c("Shark_C","Yelloweye_rockfish","Benthopel_Fish","Pisciv_S_Fish",
                "Pisciv_T_Fish","Shark_D","Shark_P")
  pelagics <- c("Pisciv_V_Fish","Demersal_S_Fish","Pacific_Ocean_Perch","Mesopel_M_Fish",
                "Planktiv_L_Fish","Jack_mackerel","Planktiv_S_Fish","Pacific_sardine",
                "Anchovy","Herring","Pisciv_B_Fish")
  demersals <- c("Demersal_P_Fish","Planktiv_O_Fish","Demersal_D_Fish",
                 "Demersal_DC_Fish","Demersal_O_Fish","Darkblotched_rockfish",
                 "Demersal_F_Fish","Demersal_E_Fish","Bocaccio_rockfish",
                 "Demersal_B_Fish","Shark_R","Mesopel_N_Fish","Shark_B","Spiny_dogfish",
                 "SkateRay")
  selflats <- c("Pisciv_D_Fish", "Arrowtooth_flounder","Petrale_sole")
}

if(initNEUS) { # Sarah's NEUS Grouping
  nontrawl <- c("Pisciv_T_Fish", "Shark_D", "Shark_P", "Reptile", "Mesopel_M_Fish")
  pelagics <- c("Planktiv_L_Fish", "Planktiv_S_Fish", "Benthopel_Fish", "Pisciv_S_Fish")
  demersals <- c("Pisciv_D_Fish", "Demersal_D_Fish","Demersal_E_Fish", 
                 "Demersal_S_Fish","Demersal_B_Fish","Demersal_DC_Fish",
                 "Demersal_O_Fish","Demersal_F_Fish",
                 "Shark_B", "SkateRay")
  selflats <- c("Pisciv_B_Fish")
}

if(initNOBA) { # Sarah's NOBA Grouping
  nontrawl <- c("Sharks_other", "Pelagic_large","Mesop_fish")
  pelagics <- c("Pelagic_small","Redfish_other","Mackerel","Haddock",
                "Saithe","Redfish","Blue_whiting","Norwegian_ssh","Capelin")
  demersals <- c("Demersals_other","Demersal_large","Flatfish_other","Skates_rays",
                 "Green_halibut","North_atl_cod","Polar_cod","Snow_crab")
  selflats <- c("Long_rough_dab")
}

The following settings are for our example standard survey once per year, most areas, with mixed efficiency and selectivity:

# general specifications for bottom trawl survey, with items defined above commented out to avoid wasting time loading already loaded files:
#   once per year at mid year
# generalized timesteps all models
runpar <- load_runprm(d.name, run.prm.file)
noutsteps <- runpar$tstop/runpar$outputstep
stepperyr <- if(runpar$outputstepunit=="days") 365/runpar$toutinc

midptyr <- round(median(seq(1,stepperyr)))

annualmidyear <- seq(midptyr, noutsteps, stepperyr)

#   ~75-80% of boxes (leave off deeper boxes?)
boxpars <- load_box(d.name, box.file)
boxsurv <- c(2:round(0.8*(boxpars$nbox - 1)))

#   define bottom trawl mixed efficiency
ef.nt <- 0.01 # for large pelagics, reef dwellers, others not in trawlable habitat
ef.pl <- 0.1  # for pelagics
ef.dm <- 0.7  # for demersals
ef.fl <- 1.1  # for selected flatfish

# bottom trawl survey efficiency specification by species group
effnontrawl <- data.frame(species=nontrawl, efficiency=rep(ef.nt,length(nontrawl)))
effpelagics <- data.frame(species=pelagics, efficiency=rep(ef.pl,length(pelagics)))
effdemersals <- data.frame(species=demersals, efficiency=rep(ef.dm,length(demersals)))
effselflats <- data.frame(species=selflats, efficiency=rep(ef.fl,length(selflats)))

efficmix <- bind_rows(effnontrawl, effpelagics, effdemersals, effselflats)

#   mixed selectivity (using 10 agecl for all species)
#     flat=1 for large pelagics, reef dwellers, others not in trawlable habitat
#     sigmoid 0 to 1 with 0.5 inflection at agecl 3 for pelagics, reaching 1 at agecl 5, flat top
#     sigmoid 0 to 1 with 0.5 inflection at agecl 5 for most demersals and flatfish, reaching 1 at agecl 7, flat top
#     dome shaped 0 to 1 at agecl 6&7 for selected demersals, falling off to 0.7 by agecl 10

sigmoid <- function(a,b,x) {
  1 / (1 + exp(-a-b*x))
}

# survey selectivity specification by species group
selnontrawl <- data.frame(species=rep(nontrawl, each=10),
                          agecl=rep(c(1:10),length(nontrawl)),
                          selex=rep(1.0,length(nontrawl)*10))
selpelagics <- data.frame(species=rep(pelagics, each=10),
                          agecl=rep(c(1:10),length(pelagics)),
                          selex=sigmoid(5,1,seq(-10,10,length.out=10)))
seldemersals <- data.frame(species=rep(demersals, each=10),
                          agecl=rep(c(1:10),length(demersals)),
                          selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selselflats <- data.frame(species=rep(selflats, each=10),
                          agecl=rep(c(1:10),length(selflats)),
                          selex=sigmoid(1,1,seq(-10,10,length.out=10)))

selexmix <- bind_rows(selnontrawl, selpelagics, seldemersals, selselflats)

# use this constant 0 cv for testing
surv_cv_0 <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))

#   define bottom trawl survey cv by group
cv.nt <- 1.0 # for large pelagics, reef dwellers, others not in trawlable habitat
cv.pl <- 0.5  # for pelagics
cv.dm <- 0.3  # for demersals
cv.fl <- 0.3  # for selected flatfish

# specify cv by species groups
surv_cv_nontrawl <- data.frame(species=nontrawl, cv=rep(cv.nt,length(nontrawl)))
surv_cv_pelagics <- data.frame(species=pelagics, cv=rep(cv.pl,length(pelagics)))
surv_cv_demersals <- data.frame(species=demersals, cv=rep(cv.dm,length(demersals)))
surv_cv_selflats <- data.frame(species=selflats, cv=rep(cv.fl,length(selflats)))

surv_cv_mix <- bind_rows(surv_cv_nontrawl, surv_cv_pelagics, surv_cv_demersals, surv_cv_selflats)

Here we use create_survey on the numbers output of run_truth to create the survey census of age composition.

survey_testNstd <- create_survey(dat = truth$nums,
                                 time = annualmidyear,
                                 species = survspp,
                                 boxes = boxsurv,
                                 effic = efficmix,
                                 selex = selexmix)

# consider saving this interim step if it takes a long time go generate

Next we apply the aggregateDensityData function can to resn and structn for survey times, species, and boxes.

# aggregate true resn per survey design
survey_aggresnstd <- aggregateDensityData(dat = truth$resn,
                                 time = annualmidyear,
                                 species = survspp,
                                 boxes = boxsurv)

# aggregate true structn per survey design
survey_aggstructnstd <- aggregateDensityData(dat = truth$structn,
                                 time = annualmidyear,
                                 species = survspp,
                                 boxes = boxsurv)

Now we should have inputs to sample_fish on the same scale, and they need to be aggregated across boxes into a single biological sample for the whole survey. The sample_fish applies the median for aggregation and does not apply multinomial sampling if sample=FALSE in the function call.

#   define n fish for biological sampling by group
#   this could easily be a vector or time series, constant here
ns.nt <- 25 # for large pelagics, reef dwellers, others not in trawlable habitat
ns.pl <- 1000  # for pelagics
ns.dm <- 1000  # for demersals
ns.fl <- 1000  # for selected flatfish

effNnontrawl <- data.frame(species=nontrawl, effN=rep(ns.nt,length(nontrawl)))
effNpelagics <- data.frame(species=pelagics, effN=rep(ns.pl,length(pelagics)))
effNdemersals <- data.frame(species=demersals, effN=rep(ns.dm,length(demersals)))
effNselflats <- data.frame(species=selflats, effN=rep(ns.fl,length(selflats)))

effNmix <- bind_rows(effNnontrawl, effNpelagics, effNdemersals, effNselflats)

# apply default sample fish as before to get numbers
numsstd <- sample_fish(survey_testNstd, effNmix)

#dont sample these, just aggregate them using median (effNmix does nothing)
structnstd <- sample_fish(survey_aggstructnstd, effNmix, sample = FALSE)

resnstd <-  sample_fish(survey_aggresnstd, effNmix, sample = FALSE)

# now cut these down to a single species for testing

ss_numsstd <- numsstd[numsstd$species == "Green_halibut",]
ss_structnstd <- structnstd[structnstd$species == "Green_halibut",]
ss_resnstd <- resnstd[resnstd$species == "Green_halibut",]

ss_length_stdsurv <- calc_age2length(structn = ss_structnstd,
                                 resn = ss_resnstd,
                                 nums = ss_numsstd,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

Plot samples one species:

Sampling looks really different for greenland halibut than census did. Check this!

I’m using a shiny app for checking, which I’m not sure I want to launch from the html, so this code block is not evaluated. Can do interactively instead. It takes a few minutes to load the info:

# lets visualize which boxes we sampled and where greenland halibut are in NOBA
# use https://github.com/Atlantis-Ecosystem-Model/shinyrAtlantis
# I've already installed it, see readme for instructions
library(shinyrAtlantis)
obj <- make.sh.prm.object(file.path(d.name, box.file), 
                          file.path(d.name, functional.groups.file), 
                          file.path(d.name, biol.prm.file))
sh.prm(obj)

Our survey uses boxes 2:47 as parameterized here (somewhat randomly for testing). This is the NOBA map (copied from the shiny app):

Greenland halibut adults have a fairly high fraction in box 1, which we are missing with our survey, and we miss 48:57 where there are high (but not the highest) fraction of juveniles (screenshot of shiny app):

Greenland halibut have 2 years per output age class. Our selectivity acts on the age class. We have them in the demersals group, so the selectivity curve is:

So we won’t see any fish until they are almost 10, when they are pretty big, and we fully sample the larger fish. This combined with missing some juvenile areas may explain the strange shapes of the sampled length frequencies.

Here is the same survey missing the areas as described, but with full selectivity for all greenland halibut age classes:

ss.name <- funct.group.names[funct.group.names == "Green_halibut"]

selex1 <- data.frame(species=rep(funct.group.names, each=10),
                     agecl=rep(c(1:10),length(funct.group.names)),
                     selex=rep(1.0,length(funct.group.names)*10))

# get survey nums with full (no) selectivity
ss_survey_testNstd_nosel <- create_survey(dat = truth$nums,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv,
                                 effic = efficmix,
                                 selex = selex1)

# now sample fish nums from this
ss_numsstd_nosel <- sample_fish(ss_survey_testNstd_nosel, effNmix)

# aggregate resn and struct n haven't changed because selectivity doesn't apply
# use the ones made above in the calc call with new nums

ss_length_stdsurv_nosel <- calc_age2length(structn = ss_structnstd,
                                 resn = ss_resnstd,
                                 nums = ss_numsstd_nosel,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

This should give the difference between missing key survey strata and survey selectivity, here set to capture all age classes perfectly:

We may want to extend the length bins for something like greenland halibut, which get larger than the hardcoded largest bin in calc_age2length, 150 cm. Both census and these sampled length comps are chopped off at 150 cm for this species.

I will likely change the calc_age2length code to pass the upper bin of 150 as a default but allow the user to change it in the function call.

Samples for all species and SAVE the file

Now run for all species in the survey.

length_stdsurv <- calc_age2length(structn = structnstd,
                                 resn = resnstd,
                                 nums = numsstd,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

#save for later use, takes a long time to generate
saveRDS(length_stdsurv, file.path(d.name, paste0(scenario.name, "length_stdsurv.rds")))

Runtime for census lengths was ~2.5 hours for NOBA once we decided once per year was plenty. For standard survey lengths runtime was not reduced because the dimensions (species, timesteps) were identical.

Compare sampled lf with census from here:

length_censussurvsamp <- readRDS(file.path(d.name, paste0(scenario.name, "length_censussurvsamp.rds")))

length_stdsurv <- readRDS(file.path(d.name, paste0(scenario.name, "length_stdsurv.rds")))

Comparison plots for other species. Blue bars are the census length sample and red outlines are the survey sample (standardized for comparison; census total n lengths is 1e+8 while sample total n lengths is 1000):

Here as with cohort sampling, the impact of survey selectivity is much larger than sampling error with effN of 1000. We could try effN that is lower, but 1000 lengths may be reasonable from a survey. A sample_lengths function could add further measuremet error on top of this, if desired (though measurement error at the 1 cm level may be negligible).