Consumption info for multispecies modeling

Multispecies models rely on either estimated or input parameters governing predator consumption rates to estimate predation mortality on consumed species.

Stomach weights are needed as input to the length-based multispecies model hydra. At present each species has a stomach weight for each size bin that is repeated for each year. We have total consumed weight for each predator age class at each timestep in the detailed diet atlantis output, so would need to map age classes to length bins, sum consumed weight by length bin and divide by total numbers in that length bin to get stomach weight.

Approach 1 based on detailedDietCheck.txt

What we have saved so far is a diet comp in proportion that has had bias and observation error added to proportions. However, we need the stomach weight portion of the equation, which has the survey timing/area/efficiency/selectivity applied to true consumed weight. The hydra input is mean daily stomach weight in grams. From the Atlantis manual, “DetailedDietCheck.txt returns the total consumed biomass since the last output given for each cell (box and layer) of each age group of each functional group.” So daily per capita consumed biomass is this output, summed over prey for total consumption, divided by the numbers from the same survey design, divided by the output step omlist_ss$runpar$outputstep.

The function below should return per-capita consumption by age class using the detailedDietCheck.txt output. This can then be converted to lengthbin as necessary.

om_cons<- function(config = configfile,
                   dietfile = file_diet,
                   usersurvey = usersurvey_file,
                   omlist_ss,
                   n_reps = n_reps,
                   save = TRUE){
  
  source(config)
  
  #Load functional groups
  fgs <- atlantisom::load_fgs(dir=d.name,
                              file_fgs = functional.groups.file)
  
  
  # load or read in saved detailed diet
  if(!file.exists(file.path(d.name,
                            paste0(scenario.name, "detaileddiet.rds")))){
    detaileddiet <- load_detailed_diet_comp(dir = d.name, 
                                            file_diet, 
                                            fgs = fgs)
    
    if(save){
      saveRDS(detaileddiet, file.path(d.name, paste0(scenario.name, "detaileddiet.rds")))
    }
    
  } else {
    detaileddiet <- readRDS(file.path(d.name,
                                      paste0(scenario.name, "detaileddiet.rds")))
  }
  
  #one script for dimension parameters to be used in multiple functions
  source("config/omdimensions.R", local = TRUE)
  
  survObsStomWt <- list()
  
  for (s in usersurvey)
  {
    source(s, local = TRUE)
    
    # survtime doesn't match units of time.days in detaileddiet
    survtime <- survey_sample_full*omlist_ss$runpar$outputstep
    
    # apply survey design to detailed diet
    survey_cons <- create_survey_diet(dat = detaileddiet,
                                      time = survtime,
                                      species = survspp,
                                      boxes = survboxes,
                                      effic = surveffic,
                                      selex = survselex)
    
    # get numbers at ageclass for same survey design
    # note different time units!
    survey_N <- atlantisom::create_survey(dat = omlist_ss$truenums_ss,
                                          time = survey_sample_full,
                                          species = survspp,
                                          boxes = survboxes,
                                          effic = surveffic,
                                          selex = survselex.agecl)
    
    # convert survey N times to cons times in days
    survey_N$time <- survey_N$time*omlist_ss$runpar$outputstep
    
    # get rid of polygon
    survey_totN <- survey_N %>%
      group_by(species, agecl, time) %>%
      summarise(totN = sum(atoutput)) %>%
      ungroup()
    
    # sum over prey to get total consumption in t, divide by N and timestep, convert to g
    survey_totcons <- survey_cons %>%
      group_by(species, agecl, time) %>%
      summarise(totcons = sum(atoutput)) %>%
      ungroup() %>%
      left_join(survey_totN) %>%
      mutate(percap_cons = totcons/totN) %>%
      #mutate(daily_percap_g = percap_cons/omlist_ss$runpar$outputstep*1000000)
      mutate(daily_percap_g = percap_cons*1000000) #try assuming cons is snapshot not cumulative since last timestep
    
    #save survey consumption, takes a long time to generate with lots of reps/species
    #if(save){
    saveRDS(survey_cons, file.path(d.name, paste0(scenario.name, "_",
                                                  survey.name, "surveycons.rds")))
    #}
    
    survObsStomWt[[survey.name]] <- survey_cons
    
  }
  
}

NOBA

The configuration files for NOBA are in posiedon-dev/config and see also ms-keyrun project repo illustrating full simulated multispecies modeling data. The NOBA model output was producing fairly low per capita consumption using the code above, so we will compare that with what we get from the CC model to see if we have an output file problem across models or if this is a characteristic of the NOBA model.

source(here("config/NOBA2Config.R"))

# this hardcodes a d.name in a different folder because I don't want multiple copies

if(!dir.exists(d.name)){
  d.name <- "/Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_march_2020"
}

# apply surveys to true total consumption (snipped from atlantosom om_diet.R and om_comps.R)
detaileddiet <- readRDS(file.path(d.name, paste0(scenario.name, "detaileddiet.rds")))

usersurvey <- c(here("config/mssurvey_spring.R"),
                here("config/mssurvey_fall.R"))

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

source(here("config/omdimensions.R"))
source(usersurvey[2])

    # survtime doesn't match units of time.days in detaileddiet
    survtime <- survey_sample_full*omlist_ss$runpar$outputstep
    
    # apply survey design to detailed diet
    survey_cons <- create_survey_diet(dat = detaileddiet,
                                      time = survtime,
                                      species = survspp,
                                      boxes = survboxes,
                                      effic = surveffic,
                                      selex = survselex)
    
    # get numbers at ageclass for same survey design
    # note different time units!
    survey_N <- atlantisom::create_survey(dat = omlist_ss$truenums_ss,
                                          time = survey_sample_full,
                                          species = survspp,
                                          boxes = survboxes,
                                          effic = surveffic,
                                          selex = survselex.agecl)
    
    # convert survey N times to cons times in days
    survey_N$time <- survey_N$time*omlist_ss$runpar$outputstep
    
    # get rid of polygon
    survey_totN <- survey_N %>%
      group_by(species, agecl, time) %>%
      summarise(totN = sum(atoutput)) %>%
      ungroup()
    
    # sum over prey to get total consumption in t, divide by N and timestep, convert to g
    survey_totcons <- survey_cons %>%
      group_by(species, agecl, time) %>%
      summarise(totcons = sum(atoutput)) %>%
      ungroup() %>%
      left_join(survey_totN) %>%
      mutate(percap_cons = totcons/totN) %>%
      #mutate(daily_percap_g = percap_cons/omlist_ss$runpar$outputstep*1000000)
      mutate(daily_percap_g = percap_cons*1000000) #try assuming cons is snapshot not cumulative since last timestep

For example, NOBA cod daily per capita intake is tens of grams for older age classes, but is 3-4x lower than expected:

DT::datatable(as.data.frame(survey_totcons%>%filter(species=="North_atl_cod")),
              rownames = FALSE,
              options = list(pageLength = 25, 
                         order = list(list(0, 'asc')))
)

CCA

This is a short run of the CCA model implemented in the latest codebase, so it can produce full age structured output. First we need to process the detailed diet file, and also get the output of run_truth for numbers. At present it seems the CC model is also producing very low per capita consumpton. This suggests that interpretation of the DetailedDietCheck.txt output is incorrect? It doesn’t seem like it is tons consumed at by the population of predators at agecl.

source(here("config/CCConfig_constest.R"))

file_diet <- paste0(scenario.name,"DetailedDietCheck.txt")

save <- TRUE

#Load functional groups
  fgs <- atlantisom::load_fgs(dir=d.name,
                              file_fgs = functional.groups.file)
  
  
  # load or read in saved detailed diet
  if(!file.exists(file.path(d.name,
                            paste0(scenario.name, "detaileddiet.rds")))){
    detaileddiet <- load_detailed_diet_comp(dir = d.name, 
                                            file_diet, 
                                            fgs = fgs)
    
    if(save){
      saveRDS(detaileddiet, file.path(d.name, paste0(scenario.name, "detaileddiet.rds")))
    }
    
  } else {
    detaileddiet <- readRDS(file.path(d.name,
                                      paste0(scenario.name, "detaileddiet.rds")))
  }
  
# Age structured single species subset
  sppsubset <- fgs %>%
    filter(Name %in% c("Yelloweye_rockfish",
                       "Darkblotched_rockfish",
                       "Bocaccio_rockfish",
                       "Pacific_Ocean_Perch",
                       "Arrowtooth_flounder",
                       "Petrale_sole",
                       "Jack_mackerel",
                       "Pacific_sardine",
                       "Anchovy",
                       "Herring",
                       "Spiny_dogfish")
    )
  

# Get true data and species subset
  
CCAtestom <- om_init(here("config/CCConfig_constest.R"))

CCAtestom_ms <- om_species(sppsubset$Name, CCAtestom)

omlist_ss <- CCAtestom_ms 

# fall census survey
source(here("config/omdimensions.R"))
source(here("config/mssurvey_CCtest.R"))
       
      
       
    # survtime doesn't match units of time.days in detaileddiet
    survtime <- survey_sample_full*omlist_ss$runpar$outputstep
    
    # apply survey design to detailed diet
    survey_cons <- create_survey_diet(dat = detaileddiet,
                                      time = survtime,
                                      species = survspp,
                                      boxes = survboxes,
                                      effic = surveffic,
                                      selex = survselex)
    
    # get numbers at ageclass for same survey design
    # note different time units!
    survey_N <- atlantisom::create_survey(dat = omlist_ss$truenums_ss,
                                          time = survey_sample_full,
                                          species = survspp,
                                          boxes = survboxes,
                                          effic = surveffic,
                                          selex = survselex.agecl)
    
    # convert survey N times to cons times in days
    survey_N$time <- survey_N$time*omlist_ss$runpar$outputstep
    
    # get rid of polygon
    survey_totN <- survey_N %>%
      group_by(species, agecl, time) %>%
      summarise(totN = sum(atoutput)) %>%
      ungroup()
    
    # sum over prey to get total consumption in t, divide by N and timestep, convert to g
    survey_totcons <- survey_cons %>%
      group_by(species, agecl, time) %>%
      summarise(totcons = sum(atoutput)) %>%
      ungroup() %>%
      left_join(survey_totN) %>%
      mutate(percap_cons = totcons/totN) %>%
      #mutate(daily_percap_g = percap_cons/omlist_ss$runpar$outputstep*1000000)
      mutate(daily_percap_g = percap_cons*1000000) #try assuming cons is snapshot not cumulative since last timestep

For example, CCA arrowtooth flounder should have a daily per capita intake in at least tens of grams for older age classes:

knitr::kable(as.data.frame(survey_totcons%>%filter(species=="Arrowtooth_flounder")))

species	agecl	time	totcons	totN	percap_cons	daily_percap_g
Arrowtooth_flounder	2	365	0.0000006	26203309.8	0e+00	0.0000000
Arrowtooth_flounder	2	730	0.0000020	26312266.5	0e+00	0.0000001
Arrowtooth_flounder	2	1095	0.0000019	24872346.5	0e+00	0.0000001
Arrowtooth_flounder	2	1460	0.0000012	15855594.8	0e+00	0.0000001
Arrowtooth_flounder	3	365	0.2361153	15708545.3	0e+00	0.0150310
Arrowtooth_flounder	3	730	0.5687898	16401427.4	0e+00	0.0346793
Arrowtooth_flounder	3	1095	0.5574132	16718155.0	0e+00	0.0333418
Arrowtooth_flounder	3	1460	0.6081520	18053962.5	0e+00	0.0336852
Arrowtooth_flounder	4	365	0.1738271	9417156.4	0e+00	0.0184586
Arrowtooth_flounder	4	730	0.4186986	9832692.3	0e+00	0.0425823
Arrowtooth_flounder	4	1095	0.4103810	10022737.4	0e+00	0.0409450
Arrowtooth_flounder	4	1460	0.4477898	10823708.4	0e+00	0.0413712
Arrowtooth_flounder	5	365	0.1167518	5645442.6	0e+00	0.0206807
Arrowtooth_flounder	5	730	0.2812394	5894574.2	0e+00	0.0477116
Arrowtooth_flounder	5	1095	0.2756387	6008552.9	0e+00	0.0458744
Arrowtooth_flounder	5	1460	0.3008103	6488936.0	0e+00	0.0463574
Arrowtooth_flounder	6	365	0.0747134	3384354.4	0e+00	0.0220761
Arrowtooth_flounder	6	730	0.1799696	3533708.3	1e-07	0.0509294
Arrowtooth_flounder	6	1095	0.1763902	3602074.4	0e+00	0.0489691
Arrowtooth_flounder	6	1460	0.1925054	3890156.2	0e+00	0.0494853
Arrowtooth_flounder	7	365	0.0465334	2028868.1	0e+00	0.0229356
Arrowtooth_flounder	7	730	0.1120871	2118401.4	1e-07	0.0529112
Arrowtooth_flounder	7	1095	0.1098614	2159432.4	1e-07	0.0508751
Arrowtooth_flounder	7	1460	0.1198996	2332138.0	1e-07	0.0514119
Arrowtooth_flounder	8	365	0.0285327	1216278.2	0e+00	0.0234590
Arrowtooth_flounder	8	730	0.0687274	1269951.3	1e-07	0.0541181
Arrowtooth_flounder	8	1095	0.0673626	1294548.2	1e-07	0.0520356
Arrowtooth_flounder	8	1460	0.0735204	1398093.4	1e-07	0.0525862
Arrowtooth_flounder	9	365	0.0173360	729137.8	0e+00	0.0237760
Arrowtooth_flounder	9	730	0.0417577	761316.8	1e-07	0.0548494
Arrowtooth_flounder	9	1095	0.0409296	776064.4	1e-07	0.0527399
Arrowtooth_flounder	9	1460	0.0446700	838137.6	1e-07	0.0532967
Arrowtooth_flounder	10	365	0.0104766	437107.4	0e+00	0.0239680
Arrowtooth_flounder	10	730	0.0252349	456395.7	1e-07	0.0552917
Arrowtooth_flounder	10	1095	0.0247339	465234.7	1e-07	0.0531644
Arrowtooth_flounder	10	1460	0.0269944	502449.0	1e-07	0.0537256

Approach 2 based on PROD.nc

An existing atlantisom function calc_pred_diet() was tested and found to produce consumption estimates orders of magnitude higher than those derived from detailedDietCheck.txt.

Here we develop new functions that don’t apply global diet composition to consumption but rather just get total consumption from PROD.nc. The challenge is figuring out what the units are of “Eat” and “Grazing” and how to appropriately expand them to the population level (water volume or numbers in a polygon).

First we need to get “Eat” and “Grazing” from PROD.nc and “Vol” from the .nc file. Here is our model setup from the config file. Lets look at NOBA_sacc_38, which is being used for mskeyrun simulated dataset generation. THIS OUTPUT IS NOT COMPARABLE TO NOBA_March_2020 above.

# WARNING!
# hardcoding because making copies of this in each repo is too cumbersome
# this currently lives in mskeyrun to make the simulated datasets with no climate
# change as appropriate!

dir <- "/Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38"

file_fgs <- "nordic_groups_v04.csv"
file_biolprm <- "nordic_biol_incl_harv_v_011_1skg.prm"
file_bgm <- "Nordic02.bgm"
file_init <- "nordic_biol_v23.nc"
file_runprm <- "nordic_run_v01.xml"
scenario <- "nordic_runresults_01"

select_groups <- mskeyrun::simFocalSpecies$Name # 11 mskeyrun fish species

verbose <- FALSE

This code is part of atlantisom::run_truth(). With our selected species from mskeyrun, only eat will be created. grazing applies to non-age structured groups.

If this works for fish, test for inverts too

  # Read in the functional groups csv since that is used by many functions
  fgs <- load_fgs(dir = dir, file_fgs = file_fgs)
  # Read in the biomass pools
  bps <- load_bps(dir = dir, fgs = file_fgs, file_init = file_init)
  # Read in the biological parameters
  biol <- load_biolprm(dir = dir, file_biolprm = file_biolprm)
  # Read in the run parameters
  runprm <- load_runprm(dir = dir, file_runprm = file_runprm)

  nc_out <- paste0(scenario, ".nc")
  nc_prod <- paste0(scenario, "PROD.nc")
  
  # Get the boundary boxes
  allboxes <- load_box(dir = dir, file_bgm = file_bgm)
  boxes <- get_boundary(allboxes)



  eat <- load_nc(dir = dir,
                     file_nc = nc_prod,
                     bps = bps,
                     fgs = fgs,
                     select_groups = select_groups,
                     select_variable = "Eat",
                     check_acronyms = TRUE,
                     bboxes = boxes)

## [1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01PROD.nc successfully"

  if(verbose) message("Eaten read in.")

  grazing <- load_nc(dir = dir,
                 file_nc = nc_prod,
                 bps = bps,
                 fgs = fgs,
                 select_groups = select_groups,
                 select_variable = "Grazing",
                 check_acronyms = TRUE,
                 bboxes = boxes)

## [1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01PROD.nc successfully"

  if(verbose) message("Grazing read in.")

  vol <- load_nc_physics(dir = dir,
                         file_nc = nc_out,
                         physic_variables = "volume",
                         aggregate_layers = FALSE,
                         bboxes = boxes)
  if(verbose) message("Volume read in.")

Now try the new function atlantisom::calc_pred_cons() which uses the volume of all water column layers to expand the atoutput (assumed to be mg N per cubic m at that timestep) to bio_eaten in tons.

This is the function code

#' Calculate eaten total biomass in tonnes for each functional group.
#'
#' Function to calculate the eaten total biomass in tonnes:
#' "consumption" for each functional group
#' per time, polygon, and ageclass.
#'
#' @param eat A \code{data.frame} containing the consumed biomass in mg N for
#'   each functional group, time, ageclass, and polygon.
#'   The \code{data.frame} must originate from \code{\link{load_nc}}
#'   using \code{select_variable = "Eat"}.
#' @param grazing A \code{data.frame} containing the consumed biomass in mg N for
#'   each functional group, time, ageclass, and polygon.
#'   The \code{data.frame} must originate from \code{\link{load_nc}}
#'   using \code{select_variable = "Grazing"}.
#' @param vol A \code{data.frame} containing the volume for
#'   each time, polygon, and layer.
#'   The \code{data.frame} must originate from \code{\link{load_nc_physics}}
#'   using \code{physic_variables = "volume"}.
#' @template biolprm
#' @return A code{data.frame} in tidy format with the following columns:
#'   species, agecl, time, polygon, atoutput (bio_eaten)
#'
#' @family calc functions
#' @importFrom magrittr %>%
#' @export
#' @author Alexander Keth, Sarah Gaichas
#'
#' @examples
#' d <- system.file("extdata", "SETAS_Example", package = "atlantisom")
#' fgs <- load_fgs(d, "Functional_groups.csv")
#' bps <- load_bps(dir = d, fgs = "Functional_groups.csv",
#'   file_init = "Initial_condition.nc")
#' runprm <- load_runprm(d, "Run_settings.xml")
#' boxes <- get_boundary(load_box(dir = d, file_bgm = "Geography.bgm"))
#' groups <- load_fgs(dir = d, "Functional_groups.csv")
#' groups <- groups[groups$IsTurnedOn > 0, "Name"]
#' eat <- load_nc(dir = d, file_nc = "outputsPROD.nc",
#'   bps = bps, fgs = fgs, select_groups = groups,
#'   select_variable = "Eat", check_acronyms = TRUE,
#'   bboxes = boxes)
#' grazing <- load_nc(dir = d, file_nc = "outputsPROD.nc",
#'   bps = bps, fgs = fgs, select_groups = groups,
#'   select_variable = "Grazing", check_acronyms = TRUE,
#'   bboxes = boxes)
#' vol <- load_nc_physics(dir = d, file_nc = "outputs.nc",
#'   physic_variables = "volume", aggregate_layers = FALSE,
#'   bboxes = boxes)
#' biolprm <- load_biolprm(dir = d,
#'   file_biolprm = "Biology.prm")
#' runprm <- load_runprm(dir = d, file_runprm = "Run_settings.xml")
#' calcs <- calc_pred_cons(eat = eat, grazing = grazing,
#'   vol = vol, biolprm = biolprm, runprm = runprm)
#' rm(calcs)
#'
calc_pred_cons <- function(eat, grazing, vol, biolprm, runprm){

  # Conversion factor from mg N to tonnes wet-weight
  bio_conv <- biolprm$redfieldcn * biolprm$kgw2d / 1000000000

  # Eat and grazing are per species, agecl, polygon, and time.
  # Therefore, we need to aggregate the vol over layers.

  # SKG: The main question is the UNITS of eat and grazing.
  # manual says this is either mg/m3/day or mg/day/individ
  # calculations here assume mg/m3/day so multiplying by box vol
  # is this m3 assuming the entire box volume or only where groups feed?
  # do we need to match species to layers to get sppropriate volume?
  # VERTnight_XXX and VERTday_XXX in biol.prm

  # all layers? even sediment? try filtering out layer > 7 to remove sediment volume
  # WARNING: ensure that layer 7 (processed via atlantisom) isn't unique to NOBA
  vol <- dplyr::filter(vol, layer<7)

  vol <- aggregate(atoutput ~ polygon + time, data = vol, sum)
  colnames(vol)[colnames(vol) == "atoutput"] <- "vol"

  # Combine eat and grazing! Calculate eaten biomass
  biomass_eaten <- rbind(eat, grazing)
  biomass_eaten <- merge(biomass_eaten, vol,
    by = c("time", "polygon"))

  # todo: previously AK noted that the bio_eaten values were too high
  # the new values need to be checked to see if they give realistic
  # results
  # removed layer 7 in volume which should expand mg N to only water?
  # total consumption and per capita from this look ok
  # snapshot so should be daily value?

  biomass_eaten$bio_eaten <- with(biomass_eaten,
    atoutput * vol * bio_conv)

  return(biomass_eaten)
}

This is the current result

biomass_eaten <- calc_pred_cons(eat = eat, grazing = grazing, vol = vol, biolprm = biol, runprm = runprm)

codtotcons <- biomass_eaten %>% 
  dplyr::filter(species=="North_atl_cod") %>%
  dplyr::group_by(time, agecl) %>% 
  dplyr::summarise(totcons = sum(bio_eaten))

Plotting a cod example total consumption by age class

ggplot2::ggplot(codtotcons, aes(x=time, y=totcons)) + geom_line() + facet_wrap(~agecl, scales="free")

Per capita consumption

Sarah Gaichas

24 June, 2022

Consumption info for multispecies modeling

Approach 1 based on detailedDietCheck.txt

NOBA

CCA

Approach 2 based on PROD.nc