New Biomass Pools Workflow

library(here)
library(stringr)
library(magrittr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(atlantisom)

Intro

Current atlantisom functions are designed primarily for age structured modeled groups, because the package was originally scoped to generate simulated data to test fishery stock assessments. It is now clear that simulated data testing can also benefit other model types, including multispecies models and food web models. Therefore, there is a need to develop functions to develop simulated data for lower trophic level biomass pools modeled in Atlantis and also included in food web models such as Rpath.

Where is biomass information for non-age structured groups output? To be parallel to survey sampling for fish, it would be good to have it by model output timestep toutinc and polygon, but depth layer is not necessary at this time.

First check content of the .nc output file:

dir <- here("atlantisoutput/NOBA_sacc_30")
file_nc <- "nordic_runresults_01.nc"

file.nc <- file.path(dir, file_nc)

 # Load ATLANTIS output!
  at_out <- RNetCDF::open.nc(con = file.nc)
  #on.exit(RNetCDF::close.nc(at_out))

    # Get info from netcdf file! (Filestructure and all variable names)
  var_names_ncdf <- sapply(seq_len(RNetCDF::file.inq.nc(at_out)$nvars - 1),
    function(x) RNetCDF::var.inq.nc(at_out, x)$name)
  n_timesteps <- RNetCDF::dim.inq.nc(at_out, 0)$length
  n_boxes     <- RNetCDF::dim.inq.nc(at_out, 1)$length
  n_layers    <- RNetCDF::dim.inq.nc(at_out, 2)$length

Which species are biomass pools, are they identified in the groups.csv?

No, biomass pools are defined as the combination of NumCohorts == 1 and GroupType not FISH, FISH_INVERT, SHARK, BIRD, MAMMAL, REPTILE. Biomass pools groups are among those below according to the Atlantis manual, ch 6. p 106:

BIOMASS POOL GROUPS
Phytoplankton
7) Small Phytoplankton SM_PHY
8) Large Phytoplankton LG_PHY (used in combination with IsSiliconDep to define diatoms)
9) Trichodesmium and Cyanobacteria TRICHO (N fixers)
Other primary producers
10) Microphtybenthos MICROPHTYBENTHOS
11) Dinoflagellates DINOFLAG
12) Phytoben PHYTOBEN (typically used to represent macroalgae)
13) Seagrass SEAGRASS
14) Turf TURF
Zooplankton
15) Small Zooplankton SM_ZOO
16) Medium Zooplankton MED_ZOO
17) Large Zooplankton LG_ZOO
18) Jellyfish JELLIES (in the past represented as LG_ZOO)
Large pelagic invetebrates
19)Cephalopod CEP
20) Prawns PWN
Infauna
21) Small Infauna SM_INF
22) Large Infauna LG_INF
Epibenthic organisms
23) Sediment epibenthic filter feeders SED_EP_FF (often used for bivalves, sponges)
24) Benthic grazers SED_EP_OTHER
25) Mobile epibenthos MOB_EP_OTHER (often used for crabs, lobster, octopus)
26) Corals CORAL
27)Sponges SPONGE
Bacteria
28) Pelagic Bacteria PL_BACT
29) Sediment Bacteria SED_BACT
Obligatory detritus groups
30) Carrion CARRION
31) Labile detritus LAB_DET
32) Refractory detritus REF_DET
33) Additional ice and land based groups are also available.
Ice dwelling
34) ICE_BACT
35) ICE_MIXOTROPHS
36) ICE_DIATOMS
37) ICE_ZOOBIOTA
Land dwelling (vegetation only for now)
38) MARSH
39) MANGROVE

For NOBA Atlantis, these are the biomass pools:

fgs <- atlantisom::load_fgs(here("atlantisoutput/NOBA_sacc_30"), "nordic_groups_v04.csv") 

pools <- fgs %>%
  dplyr::filter(NumCohorts==1) %>%
  dplyr::select(Code, Name, Long.Name, NumCohorts, InvertType)

knitr::kable(pools)

Code	Name	Long.Name	NumCohorts	InvertType
KCR	King_crab	Red king crab	1	LG_INF
ZG	Gel_zooplankton	Gelatinous zooplankton	1	LG_ZOO
ZL	Large_zooplankton	Large zooplankton	1	LG_ZOO
ZM	Medium_zooplankton	Medium zooplankton	1	MED_ZOO
ZS	Small_zooplankton	Small zooplankton	1	SM_ZOO
DF	Dinoflag	Dinoflagellates	1	DINOFLAG
PS	Small_phytop	Small phytoplankton	1	SM_PHY
PL	Large_phytop	Large phytoplankton	1	LG_PHY
BC	Predatory_ben	Predatory benthos	1	LG_INF
BD	Detrivore_ben	Detrivore benthos	1	LG_INF
BFF	Benthic_filterf	Benthic filter feeders	1	SED_EP_FF
SPO	Sponges	Sponges	1	SED_EP_FF
COR	Corals	Corals	1	SED_EP_FF
PB	Pel_bact	Pelagic bacteria	1	PL_BACT
BB	Ben_bact	Sediment bacteria	1	SED_BACT
DR	Ref_det	Refractory detritus	1	REF_DET
DL	Lab_det	Labile detritus	1	LAB_DET
DC	Carrion	Carrion	1	CARRION

Here are the relevant outputs in the output.nc file for these groups:

# return everything in var_names_ncdf matching Name in pools

var_names_ncdf[str_detect(var_names_ncdf, str_c(pools$Name, collapse ="|"))]

##  [1] "Large_phytop_N"        "Small_phytop_N"        "Dinoflag_N"           
##  [4] "King_crab_N"           "Gel_zooplankton_N"     "Large_zooplankton_N"  
##  [7] "Medium_zooplankton_N"  "Small_zooplankton_N"   "Pel_bact_N"           
## [10] "Ben_bact_N"            "Detrivore_ben_N"       "Predatory_ben_N"      
## [13] "Ref_det_N"             "Lab_det_N"             "Carrion_N"            
## [16] "Benthic_filterf_N"     "Sponges_N"             "Corals_N"             
## [19] "Benthic_filterf_Cover" "Sponges_Cover"         "Corals_Cover"

So we want “N” nitrogen output. In theory this can be achieved with atlantisom::load_ncwith select_variable = "N". Looks like it works, output is N for our selected biomass pool species over all polygons, layers, and times. Output of unique(testNpools$species):

select_groups <- pools$Name
nc_out <- file_nc
bps <- load_bps(dir = dir, fgs = "nordic_groups_v04.csv", file_init = "nordic_biol_v23.nc")
  # Get the boundary boxes
  allboxes <- load_box(dir = dir, file_bgm = "Nordic02.bgm")
  boxes <- get_boundary(allboxes)



testNpools <- load_nc(dir = dir,
                  file_nc = nc_out,
                  bps = bps,
                  fgs = fgs,
                  select_groups = select_groups,
                  select_variable = "N",
                  check_acronyms = TRUE,
                  bboxes = boxes)

## [1] "Read /Users/sarah.gaichas/Documents/0_Data/Atlantis/Poseidon/poseidon-dev/atlantisoutput/NOBA_sacc_30/nordic_runresults_01.nc successfully"

unique(testNpools$species)

##  [1] "Benthic_filterf"    "Sponges"            "Corals"            
##  [4] "King_crab"          "Gel_zooplankton"    "Large_zooplankton" 
##  [7] "Medium_zooplankton" "Small_zooplankton"  "Dinoflag"          
## [10] "Small_phytop"       "Large_phytop"       "Predatory_ben"     
## [13] "Detrivore_ben"      "Pel_bact"           "Ben_bact"          
## [16] "Ref_det"            "Lab_det"            "Carrion"

I don’t think we have a function that aggregates over layers and converts N to biomass but there are pieces in other functions.

New atlantisom function: calc_biomass_pool based on calc_biomass_age, but adding the expansion to volume used in calc_pred_cons.

With the following exceptions!

Benthic groups have N expanded to area

LG_PHY type has N expanded to only water column, not sediment layers.

# we need the volume of each layer as input to this
vol <- load_nc_physics(dir = dir,
                         file_nc = nc_out,
                         physic_variables = "volume",
                         aggregate_layers = FALSE,
                         bboxes = boxes)

# we need area for the benthic inverts instead of volume
# THANK YOU JOE for this tip
area <- load_boxarea(dir = dir,
                     file_bgm = "Nordic02.bgm")

# we need a list of which group types to expand to area instead of volume
benthos <- c("SED_EP_FF", "SED_EP_OTHER", "SED_BAC", 
             "MOB_EP_OTHER", "LG_INF", "SM_INF", 
             "MICROPHYTBENTHOS", "SEAGRASS")



calc_biomass_pool <- function(pooln, vol, area, fgs, biolprm){
  datalist <- list(pooln)

  # Conversion factor from mg N to t wet-weight
  # should only use conversion for non vertebrates
  # also need volume of cell info
  bio_conv <- biolprm$redfieldcn * biolprm$kgw2d / 1000000000
  
  # get fgs info
    # use grouptype column to allocate
  colnames(fgs) <- tolower(colnames(fgs))

  # check for GroupType or InvertType
  if (!("grouptype" %in% colnames(fgs) | "inverttype"%in% colnames(fgs))) {
    stop(paste("The columns GroupType or InvertType ars not in your functional\n",
               "groups file."))
  }

  # change inverttype to grouptype, contents should be the same
  if("inverttype" %in% names(fgs)) names(fgs)[names(fgs) == 'inverttype'] <- 'grouptype'

  fgs$grouptype <- tolower(fgs$grouptype)

  
   pooltype <- fgs |>
    dplyr::select(species=name, grouptype) |>
    dplyr::mutate(pooltype = dplyr::case_when((grouptype %in% c("lg_inf",
                                                                "sm_inf",
                                                                "sed_bact",
                                                                "sed_ep_ff",
                                                                "sed_ep_other",
                                                                "mob_ep_other",
                                                                "coral",
                                                                "sponge")) ~ "benthos",
                                              (grouptype %in% c("lg_phy")) ~ "lg_phy",
                                              TRUE ~ "alllayers"))

  data_names <- c("species", "agecl", "polygon", "layer", "time", "atoutput")

  if (all(sapply(datalist, function(x) all(is.element(names(x), data_names))))){

    names(vol)[names(vol) == "atoutput"] <- "vol"

    sedlayer <- max(vol$layer)

    pooln <- merge(pooln, vol,
                   by = c("time", "polygon", "layer"))

    pooln <- merge(pooln, area,
                   by = c("polygon"))

    pooln <- merge(pooln, pooltype,
                   by = c("species"))

    #pooln$atoutput <- with(pooln, vol * atoutput * bio_conv)

    pooln <- pooln |>
      dplyr::mutate(atoutput = dplyr::case_when(pooltype == "alllayers" ~ vol * atoutput * bio_conv,
                                                pooltype == "lg_phy" ~
                                                  ifelse(layer!=sedlayer, vol * atoutput * bio_conv, 0),
                                                pooltype == "benthos" ~
                                                  ifelse(layer==sedlayer, area * atoutput * bio_conv, 0)))

    # Sum over layers 
    biomass_pools <- aggregate(atoutput ~ species + agecl + time + polygon,
      data = pooln, sum)
  } else {
    stop(paste("Dataframe names do not match with", data_names))
  }

  return(biomass_pools)
}

biolprm <- load_biolprm(dir = dir, file_biolprm = "nordic_biol_incl_harv_v_011_1skg.prm")

## Warning in `[<-.data.frame`(`*tmp*`, , 3:(maxmat + 2), value = structure(list(:
## provided 20 variables to replace 10 variables

biomass_pools <- calc_biomass_pool(pooln = testNpools,
                                   vol = vol,
                                   area = area,
                                   fgs = fgs,
                                   biolprm = biolprm)

Do the new biomass pools match the BiomIndx.txt output when aggregated over polygon? If so they can be the input to the survey function for time and spatial subsetting.

aggbiopools <- biomass_pools %>%
  dplyr::group_by(species, agecl, time) %>%
  dplyr::summarise(totpool = sum(atoutput))

## `summarise()` has grouped output by 'species', 'agecl'. You can override using
## the `.groups` argument.

txtbiopools <- atlantisom::load_bioind(dir, "nordic_runresults_01BiomIndx.txt", fgs) %>%
  dplyr::mutate(time = time/73) %>%
  dplyr::filter(species %in% unique(aggbiopools$species)) %>%
  dplyr::select(species, time, atoutput)

ggplot(txtbiopools, aes(x=time, y=atoutput)) +
  geom_line() +
  geom_point(data = aggbiopools, aes(x = time, y = totpool), colour = "blue", alpha=0.1) +
  theme_bw() +
  facet_wrap(~species, scales = "free_y")

Diagnose mismatches between biopools and BiomInd.txt output

Who matches exactly (enough). Patterns in absolute mismatch are interesting, some species are really close and some always way off. Volume wrong, depth layers wrong? Mine are always higher than the txt output when they are off, so maybe too much expansion somewhere.

Joe pointed out that benthic groups should be expanded to area rather than volume. That correction has now been applied, and matches are better, but still diverge over time for some groups. Not enough to worry about I think.

Originally, large phytoplankton was way off, but Joe also pointed out that LG_PHY N can be in the sediment layer but shouldn’t “count” as biomass. Therefore, this group is expanded to all layers but the sediment layer and that creates a much better match.

mismatch <- txtbiopools |>
  dplyr::left_join(aggbiopools) |>
  dplyr::mutate(mismatch = atoutput - totpool,
                bad = ifelse(abs(mismatch) > 0.01*(atoutput), 1, 0))

## Joining, by = c("species", "time")

ggplot(mismatch, aes(x=time, y=mismatch)) +
  geom_line() +
  theme_bw() +
  facet_wrap(~species, scales = 'free_y')

This is mismatch coded as 0 if less than 1% of txt output, 1 if more.

ggplot(mismatch, aes(x=time, y=bad)) +
  geom_line() +
  theme_bw() +
  facet_wrap(~species)

Large phytoplankton formerly had the largest magnitude mismatch across all groups, now fixed by not expanding N in the sediment layer for this group type.

ggplot(txtbiopools |> filter(species %in% c("Large_phytop", "Small_phytop")), 
       aes(x=time, y=atoutput)) +
  geom_line() +
  geom_point(data = aggbiopools |> filter(species %in% c("Large_phytop", "Small_phytop")), 
             aes(x = time, y = totpool), colour = "blue") +
  theme_bw() +
  facet_wrap(~species, scales = "free_y")

These groups seem close enough to me now.

ggplot(txtbiopools |> filter(species %in% c("Predatory_ben", "Detrivore_ben")), 
       aes(x=time, y=atoutput)) +
  geom_line() +
  geom_point(data = aggbiopools |> filter(species %in% c("Predatory_ben", "Detrivore_ben")), 
             aes(x = time, y = totpool), colour = "blue", alpha=0.1) +
  theme_bw() +
  facet_wrap(~species, scales = "free_y")

Try NEUS?

dir <- here("atlantisoutput/NEUSv2.1.0")
file_nc <- "neus_output.nc"

file.nc <- file.path(dir, file_nc)

 # Load ATLANTIS output!
  at_out <- RNetCDF::open.nc(con = file.nc)
  #on.exit(RNetCDF::close.nc(at_out))

    # Get info from netcdf file! (Filestructure and all variable names)
  var_names_ncdf <- sapply(seq_len(RNetCDF::file.inq.nc(at_out)$nvars - 1),
    function(x) RNetCDF::var.inq.nc(at_out, x)$name)
  n_timesteps <- RNetCDF::dim.inq.nc(at_out, 0)$length
  n_boxes     <- RNetCDF::dim.inq.nc(at_out, 1)$length
  n_layers    <- RNetCDF::dim.inq.nc(at_out, 2)$length

For NEUS Atlantis, these are the biomass pools:

fgs <- atlantisom::load_fgs(here("atlantisoutput/NEUSv2.1.0"), "neus_groups.csv") 

pools <- fgs %>%
  dplyr::filter(NumCohorts==1) %>%
  dplyr::select(Code, Name, LongName, NumCohorts, GroupType)

knitr::kable(pools)

Code	Name	LongName	NumCohorts	GroupType
SCA	Scallop	Sea scallop	1	SED_EP_FF
QHG	Quahog	Ocean quahog	1	SED_EP_FF
CLA	Surf_Clam	Atlantic surf clam	1	SED_EP_FF
BFF	Filter_Other	Other benthic filter feeder	1	SED_EP_FF
BG	Benthic_grazer	Benthic grazer	1	SED_EP_OTHER
LOB	Lobster	Lobster	1	MOB_EP_OTHER
RCB	Red_Crab	Red deep-sea crab	1	MOB_EP_OTHER
BMS	Macrobenth_Shallow	Shallow macrozoobenthos	1	MOB_EP_OTHER
ZL	Carniv_Zoo	Carnivorous zooplankton	1	LG_ZOO
BD	Deposit_Feeder	Deposit Feeder	1	LG_INF
MA	Macroalgae	Macroalgae	1	PHYTOBEN
MB	MicroPB	Microphtybenthos	1	MICROPHTYBENTHOS
SG	Seagrass	Seagrass	1	SEAGRASS
BC	Benthic_Carniv	Benthic Carnivore	1	MOB_EP_OTHER
ZG	Gelat_Zoo	Gelatinous zooplankton	1	LG_ZOO
PL	Diatom	Diatom	1	LG_PHY
DF	DinoFlag	Dinoflagellates	1	DINOFLAG
PS	PicoPhytopl	Pico-phytoplankton	1	SM_PHY
ZM	Zoo	Mesozooplankton	1	MED_ZOO
ZS	MicroZoo	Microzooplankton	1	SM_ZOO
PB	Pelag_Bact	Pelagic Bacteria	1	PL_BACT
BB	Sed_Bact	Sediment Bacteria	1	SED_BACT
BO	Meiobenth	Meiobenthos	1	SM_INF
DL	Lab_Det	Labile detritus	1	LAB_DET
DR	Ref_Det	Refractory detritus	1	REF_DET
DC	Carrion	Carrion	1	CARRION

Here are the relevant outputs in the output.nc file for these groups:

# return everything in var_names_ncdf matching Name in pools

var_names_ncdf[str_detect(var_names_ncdf, str_c(pools$Name, collapse ="|"))]

##  [1] "Carniv_Zoo_N"         "Carrion_N"            "Deposit_Feeder_N"    
##  [4] "Diatom_N"             "Diatom_S"             "DinoFlag_N"          
##  [7] "Gelat_Zoo_N"          "Lab_Det_N"            "Meiobenth_N"         
## [10] "MicroPB_N"            "MicroPB_S"            "MicroZoo_N"          
## [13] "Pelag_Bact_N"         "PicoPhytopl_N"        "Ref_Det_N"           
## [16] "Sed_Bact_N"           "Zoo_N"                "Benthic_Carniv_N"    
## [19] "Benthic_grazer_N"     "Filter_Other_Cover"   "Filter_Other_N"      
## [22] "Lobster_N"            "Macroalgae_Cover"     "Macroalgae_N"        
## [25] "Macrobenth_Shallow_N" "MicroPB_Cover"        "Quahog_Cover"        
## [28] "Quahog_N"             "Scallop_Cover"        "Scallop_N"           
## [31] "Seagrass_Cover"       "Seagrass_N"           "Surf_Clam_Cover"     
## [34] "Surf_Clam_N"

NEUS has some Cover types as well as N outputs.

select_groups <- pools$Name
nc_out <- file_nc
bps <- load_bps(dir = dir, fgs = "neus_groups.csv", file_init = "neus_init.nc")
  # Get the boundary boxes
  allboxes <- load_box(dir = dir, file_bgm = "neus_tmerc_RM2.bgm")
  boxes <- get_boundary(allboxes)



testNpools <- load_nc(dir = dir,
                  file_nc = nc_out,
                  bps = bps,
                  fgs = fgs,
                  select_groups = select_groups,
                  select_variable = "N",
                  check_acronyms = TRUE,
                  bboxes = boxes)

## Warning in load_nc(dir = dir, file_nc = nc_out, bps = bps, fgs = fgs, select_groups = select_groups, : Some selected groups are not active in the model run. Check 'IsTurnedOn' in fgs
##  Red_Crab

## [1] "Read /Users/sarah.gaichas/Documents/0_Data/Atlantis/Poseidon/poseidon-dev/atlantisoutput/NEUSv2.1.0/neus_output.nc successfully"

unique(testNpools$species)

##  [1] "Scallop"            "Quahog"             "Surf_Clam"         
##  [4] "Filter_Other"       "Benthic_grazer"     "Lobster"           
##  [7] "Macrobenth_Shallow" "Macroalgae"         "Seagrass"          
## [10] "Benthic_Carniv"     "Carniv_Zoo"         "Deposit_Feeder"    
## [13] "MicroPB"            "Gelat_Zoo"          "Diatom"            
## [16] "DinoFlag"           "PicoPhytopl"        "Zoo"               
## [19] "MicroZoo"           "Pelag_Bact"         "Sed_Bact"          
## [22] "Meiobenth"          "Lab_Det"            "Ref_Det"           
## [25] "Carrion"

# we need the volume of each layer as input to this
vol <- load_nc_physics(dir = dir,
                         file_nc = nc_out,
                         physic_variables = "volume",
                         aggregate_layers = FALSE,
                         bboxes = boxes)

# we need area for the benthic inverts instead of volume
# THANK YOU JOE for this tip
area <- load_boxarea(dir = dir,
                     file_bgm = "neus_tmerc_RM2.bgm")

biolprm <- load_biolprm(dir = dir, file_biolprm = "at_biology.prm")

## Warning in `[<-.data.frame`(`*tmp*`, , 3:(maxmat + 2), value = structure(list(:
## provided 20 variables to replace 10 variables

## Warning in load_biolprm(dir = dir, file_biolprm = "at_biology.prm"): NAs
## introduced by coercion

biomass_pools <- calc_biomass_pool(pooln = testNpools,
                                   vol = vol,
                                   area = area,
                                   fgs = fgs,
                                   biolprm = biolprm)

Do the new biomass pools match the BiomIndx.txt output when aggregated over polygon? If so they can be the input to the survey function for time and spatial subsetting.

aggbiopools <- biomass_pools %>%
  dplyr::group_by(species, agecl, time) %>%
  dplyr::summarise(totpool = sum(atoutput))

## `summarise()` has grouped output by 'species', 'agecl'. You can override using
## the `.groups` argument.

txtbiopools <- atlantisom::load_bioind(dir, "neus_outputBiomIndx.txt", fgs) %>%
  dplyr::mutate(time = time/73) %>%
  dplyr::filter(species %in% unique(aggbiopools$species)) %>%
  dplyr::select(species, time, atoutput)

ggplot(txtbiopools, aes(x=time, y=atoutput)) +
  geom_line() +
  geom_point(data = aggbiopools, aes(x = time, y = totpool), colour = "blue", alpha=0.1) +
  theme_bw() +
  facet_wrap(~species, scales = "free_y")

mismatch

mismatch <- txtbiopools |>
  dplyr::left_join(aggbiopools) |>
  dplyr::mutate(mismatch = atoutput - totpool,
                bad = ifelse(abs(mismatch) > 0.01*(atoutput), 1, 0))

## Joining, by = c("species", "time")

ggplot(mismatch, aes(x=time, y=mismatch)) +
  geom_line() +
  theme_bw() +
  facet_wrap(~species, scales = 'free_y')

This is mismatch coded as 0 if less than 1% of txt output, 1 if more.

ggplot(mismatch, aes(x=time, y=bad)) +
  geom_line() +
  theme_bw() +
  facet_wrap(~species)

Learned that we should not hardcode the sediment layer! Added a sedlayer variable to the function.

With corrections for LG_PHY things now line up in both models.

ggplot(txtbiopools |> filter(species %in% c("Diatom", "DinoFlag")), 
       aes(x=time, y=atoutput)) +
  geom_line() +
  geom_point(data = aggbiopools |> filter(species %in% c("Diatom", "DinoFlag")), 
             aes(x = time, y = totpool), colour = "blue") +
  theme_bw() +
  facet_wrap(~species, scales = "free_y")

Remaining changes to workflow, once corrected

Add calc_biomass_pool to run_truth

With an if statement? Or let the species selection take care of that.

Decision: add it and if no biomass pools species selected it should return nothing

Added to run_truth:

• area using existing load_boxarea
• pooln using existing load_nc with select_variable = “N”
• biomass_pools using new calc_biomass_pool

Returns biomass_pools as part of truth list

Testing–have to ask for N data for only biomass pool groups. If there are none we should trap this.

dir <- here("atlantisoutput","NOBA_sacc_30")
file_fgs <- "nordic_groups_v04.csv"
biomass.pools.file <- "nordic_biol_v23.nc"
file_biolprm <- "nordic_biol_incl_harv_v_011_1skg.prm"
file_bgm <- "Nordic02.bgm"
file_init <- "nordic_biol_v23.nc"
file_runprm <- "nordic_run_v04.xml"
scenario <- "nordic_runresults_01"
verbose <- TRUE
annage <- TRUE
file_fish <- "NoBAFisheries.csv"
  

  # Read in information
  # Read in the functional groups csv since that is used by many functions
  fgs <- load_fgs(dir = dir, file_fgs = file_fgs)
  
  
    #Get just the names of active functional groups
  funct.group.names <- fgs %>%
    filter(IsTurnedOn == 1) %>%
    select(Name) %>%
    .$Name
  
  
  select_groups <- funct.group.names
  
  
  # 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_catch <- paste0(scenario, 'CATCH.nc')
  dietcheck <- paste0(scenario, 'DietCheck.txt')
  nc_out <- paste0(scenario, ".nc")
  nc_prod <- paste0(scenario, "PROD.nc")
  file_catchfish <- file.path(dir,
    paste0(scenario, "CatchPerFishery.txt"))
  file_catch <- paste0(scenario, "Catch.txt")

  
  if(annage){
    if(!file.exists(paste0(file.path(dir,paste0(scenario, 'ANNAGEBIO.nc'))))){
      stop("ANNAGEBIO.nc file not found")
    }
    if(!file.exists(paste0(file.path(dir,paste0(scenario, 'ANNAGECATCH.nc'))))){
      stop("ANNAGECATCH.nc file not found")
    }
    nc_annagebio <- paste0(scenario, 'ANNAGEBIO.nc')
    nc_annagecatch <- paste0(scenario, 'ANNAGECATCH.nc')
  }

  # Get the boundary boxes
  allboxes <- load_box(dir = dir, file_bgm = file_bgm)
  boxes <- get_boundary(allboxes)

  # Get box area for benthic biomass pool calc
  area <- load_boxarea(dir = dir,
                       file_bgm = file_bgm)
    
  #Extract from NetCDF files
  # Need: dir, file_nc, bps, fgs, select_groups, select_variable,
  # check_acronyms, bboxes

  nums <- load_nc(dir = dir,
                  file_nc = nc_out,
                  bps = bps,
                  fgs = fgs,
                  select_groups = select_groups,
                  select_variable = "Nums",
                  check_acronyms = TRUE,
                  bboxes = boxes)
  if(verbose) message("Numbers read in.")

  resn <- load_nc(dir = dir,
                  file_nc = nc_out,
                  bps = bps,
                  fgs = fgs,
                  select_groups = select_groups,
                  select_variable = "ResN",
                  check_acronyms = TRUE,
                  bboxes = boxes)
  if(verbose) message("Reserve nitrogen read in.")

  structn <- load_nc(dir = dir,
                  file_nc = nc_out,
                  bps = bps,
                  fgs = fgs,
                  select_groups = select_groups,
                  select_variable = "StructN",
                  check_acronyms = TRUE,
                  bboxes = boxes)
  if(verbose) message("Structural nitrogen read in.")

  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)
  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)
  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.")
  
  # biomass pools added 2023 for Rpath testing
  pools <- fgs %>%
    dplyr::filter(NumCohorts == 1,
                  IsTurnedOn == 1) %>%
    dplyr::select(Name) %>%
    .$Name

  pooln <- load_nc(dir = dir,
                   file_nc = nc_out,
                   bps = bps,
                   fgs = fgs,
                   select_groups = pools,
                   select_variable = "N",
                   check_acronyms = TRUE,
                   bboxes = boxes)
 if(verbose) message("Biomass pools N read in.")

Test wrapper functions with biomass pool groups

Only need species and index here, comps don’t apply.

In om_species need to check for biomass pools output and append that to the truebio_ss object? If one object can use existing om_index on bio pools. However, will that be a problem if trying to combine the truenums and truebio objects? They won’t be the same length. Perhaps leave separate at this level and combine only if needed.

If a single true bio object is created, may have to modify selectivty definition for biomass pools in om_index. If left separate can apply no selectivity to biomass pool groups, only efficiency.

Decision: leave the truebiopool_ss object separate and process separately with a different survey config file and a separate call to om_index. Or just bind them together for biomass based surveys if it doesn’t break anything.

Running om_init from the mskeyrun repo with my local NOBA atlantis outputs produced expected output:

> NOBAom <- om_init(here("data-raw/simulated-data/config/NOBA_sacc38Config.R"))
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01PROD.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01PROD.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01CATCH.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01CATCH.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01CATCH.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01ANNAGEBIO.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01ANNAGECATCH.nc successfully"
[1] "Read /Users/sarah.gaichas/Documents/0_Data/ms-keyrun/simulated-data/atlantisoutput/NOBA_sacc_38/nordic_runresults_01ANNAGECATCH.nc successfully"
> NOBAom$truth$biomass_pools
               species agecl time polygon     atoutput
1             Ben_bact     1    0       1 1.376317e+05
2        Detrivore_ben     1    0       1 1.088091e+05
3             Dinoflag     1    0       1 1.312067e+04
4      Gel_zooplankton     1    0       1 1.569070e+05
5              Lab_det     1    0       1 2.581300e+07
6         Large_phytop     1    0       1 3.268895e+05
7    Large_zooplankton     1    0       1 4.772586e+05
8   Medium_zooplankton     1    0       1 1.118751e+06
9             Pel_bact     1    0       1 3.268895e+04
10       Predatory_ben     1    0       1 1.088091e+05
11             Ref_det     1    0       1 8.904920e+06
12        Small_phytop     1    0       1 3.280167e+05
13   Small_zooplankton     1    0       1 9.152906e+05
14             Sponges     1    0       1 4.621541e+05
15            Ben_bact     1    1       1 8.875378e+05
16     Benthic_filterf     1    1       1 1.163897e-04
17             Carrion     1    1       1 1.700432e-10
18              Corals     1    1       1 1.148013e-04
19       Detrivore_ben     1    1       1 1.096034e+05
20            Dinoflag     1    1       1 4.509383e+01
21     Gel_zooplankton     1    1       1 1.563221e+05
22           King_crab     1    1       1 0.000000e+00
23             Lab_det     1    1       1 8.511963e+07
24        Large_phytop     1    1       1 8.999719e+03
25   Large_zooplankton     1    1       1 5.237310e+05
26  Medium_zooplankton     1    1       1 5.390765e+06
27            Pel_bact     1    1       1 1.230456e+06
28       Predatory_ben     1    1       1 1.089640e+05
29             Ref_det     1    1       1 8.142682e+06
30        Small_phytop     1    1       1 1.177500e+03
31   Small_zooplankton     1    1       1 7.373987e+07
32             Sponges     1    1       1 4.647222e+05

Added lines to om_species to include the truebiopool outputs:

#'Species-specific atlantisom outputs
#'@description A wrapper function to reduce a set of \code{atlantisom} objects from atlantis
#'output to a single species or subset of species for creating assessment input datasets.
#'Optionally includes diet data (set diet=TRUE in function call, default is FALSE)
#'@param species The species to sample in the survey and fishery (a vector)
#'@param omlist output of \code{om_init}
#'@return Returns an omlist_ss list object containing dataframes and lists:
#' \itemize{
#'  \item{species_ss, species name(s)}
#'  \item{code_ss, species code from funct.groups}
#'  \item{truetotbio_ss, dataframe output of \code{load_bioind}, species_ss only}
#'  \item{truecatchbio_ss, dataframe output of \code{load_catch}, species_ss only}
#'  \item{YOY_ss, dataframe young of year output, code_ss only}
#'  \item{truenums_ss, numbers at age output of \code{run_truth}, species_ss only}
#'  \item{truebio_ss, biomass at age output of \code{run_truth}, species_ss only}
#'  \item{truebiopool_ss, biomass pool groups output of \code{run_truth}, species_ss only}
#'  \item{trueresn_ss, reserve nitrogen output of \code{run_truth}, species_ss only}
#'  \item{truestructn_ss, structural nitrogen output of \code{run_truth}, species_ss only}
#'  \item{truecatchnum_ss, fishery catch at age output of \code{run_truth}, species_ss only}
#'  \item{truecons_ss, total consumption output of \code{run_truth}, species_ss only}
#'  \item{truecatchtons_ss, fishery catch by fleet output of \code{run_truth}, code_ss only}
#'  \item{truedisctons_ss, fishery discard by fleet output of \code{run_truth}, code_ss only}
#'  \item{truenumsage_ss, true numbers at annual age output of \code{run_truth}, species_ss only}
#'  \item{truecatchage_ss, fishery catch at annual age output of \code{run_truth}, species_ss only}
#'  \item{truediscage_ss, fishery discard at annual age output of \code{run_truth}, species_ss only}
#'  \item{funct.groups_ss, dataframe of species characteristics, species_ss only}
#'  \item{biol, list of biological parameters passed from input omlist}
#'  \item{boxpars, dataframe of spatial parameters}
#'  \item{runpar, list of run parameters passed from input omlist}
#' },
#'
#'@export
#'
#'@family wrapper functions
#'@author Sarah Gaichas
#'
#'@examples
#'\dontrun{
#' # assuming CC3om is output of om_init(here("config/CC3config.r"))
#'CC3om_sardine <- om_species(c("Pacific_sardine"), CC3om)
#'
#'CC3om_2spp <- om_species(c("Pacific_sardine", "Mesopel_M_Fish"), CC3om)
#'}
#'
om_species <- function(species = spp, omlist, save = TRUE,
                       removefullom = TRUE){
  # spp format c("speciesname1", "speciesname2")
  if(!all(species %in% omlist$funct.group.names)) stop("species name not found")
  species_ss <- species

  #subset species true bio
  truetotbio_ss <- omlist$truetotbio[omlist$truetotbio$species %in% species_ss,]

  #subset species true catch
  truecatchbio_ss <- omlist$truecatchbio[omlist$truecatchbio$species %in% species_ss,]

  #subset species YOY
  # get code matching species name to split YOY file
  code_ss <- omlist$funct.groups$Code[which(omlist$funct.groups$Name %in% species_ss)]
  # cut to a single species in YOY file
  YOY_ss <- omlist$YOY %>%
    select(any_of(c("Time", paste0(code_ss, ".0"))))
  # reformat to be like all the other objects

  # numbers at agecl at full resolution (all polygons and layers)
  truenums_ss <- omlist$truth$nums[omlist$truth$nums$species %in% species_ss,]

  # biomass at agecl at full resolution (all polygons and layers)
  truebio_ss <- omlist$truth$biomass_ages[omlist$truth$biomass_ages$species %in% species_ss,]

  # biomass pool groups at full resolution (all polygons and layers)
  truebiopool_ss <- omlist$truth$biomass_pools[omlist$truth$biomass_pools$species %in% species_ss,]

  # reserve nitrogen at agecl at full resolution
  trueresn_ss <- omlist$truth$resn[omlist$truth$resn$species %in% species_ss,]

  # structural nitrogen at agecl at full resolution
  truestructn_ss <- omlist$truth$structn[omlist$truth$structn$species %in% species_ss,]

  # catch in numbers at agecl at full resoluation (all polygons, no layer in output)
  truecatchnum_ss <- omlist$truth$catch[omlist$truth$catch$species %in% species_ss,]

  # catch in tons at full resolution by fleet (all polygons, no layer or agecl)
  truecatchtons_ss <- omlist$truth$catchtons[omlist$truth$catchtons$species %in% code_ss,]

  # discard in tons at full resolution by fleet (all polygons, no layer or agecl)
  truedisctons_ss <- omlist$truth$disctons[omlist$truth$disctons$species %in% code_ss,]

  # consumption (biomass_eaten) at agecl at full resolution (all polygons, no layer in output)
  # based on atlantisom::calc_pred_cons()
  truecons_ss <- omlist$truth$biomass_eaten[omlist$truth$biomass_eaten$species %in% species_ss,]

  # if annage output exists, add in, otherwise fill with NULL
  if("numsage" %in% names(omlist$truth)){
    truenumsage_ss <- omlist$truth$numsage[omlist$truth$numsage$species %in% species_ss,]
  } else {truenumsage_ss <- NULL}

  if("catchage" %in% names(omlist$truth)){
    truecatchage_ss <- omlist$truth$catchage[omlist$truth$catchage$species %in% species_ss,]
  } else {truecatchage_ss <- NULL}

  if("discage" %in% names(omlist$truth)){
    truediscage_ss <- omlist$truth$discage[omlist$truth$discage$species %in% species_ss,]
  } else {truediscage_ss <- NULL}

  # original biomass_eaten was not correct by polygon, and also appeared inflated
  # a separate set of functions will get true diet by polygon from the
  # DetailedDietCheck.txt file
  # now biomass_eaten is based on atlantisom::calc_pred_cons() and included above by default

  # if(diet){
  #   truebioeaten_ss <- omlist$truth$biomass_eaten[omlist$truth$biomass_eaten$species %in% species_ss,]
  # } else {truebioeaten_ss <- NULL}

  #subset species biol parameters? no, may miss some globals that aren't by species

  #subset species functional group info
  funct.groups_ss <- omlist$funct.groups[omlist$funct.groups$Code %in% code_ss,]

  #keep the runpar and also need boxes for survey selection

  omlist_ss <- list("species_ss" = species_ss,
                    "code_ss" = code_ss,
                    "truetotbio_ss" = truetotbio_ss,
                    "truebiopool_ss" = truebiopool_ss,
                    "truecatchbio_ss" = truecatchbio_ss,
                    "YOY_ss" = YOY_ss,
                    "truenums_ss" = truenums_ss,
                    "truebio_ss" = truebio_ss,
                    "trueresn_ss" = trueresn_ss,
                    "truestructn_ss" = truestructn_ss,
                    "truecatchnum_ss" = truecatchnum_ss,
                    "truecons_ss" = truecons_ss,
                    "truecatchtons_ss" = truecatchtons_ss,
                    "truedisctons_ss" = truedisctons_ss,
                    "truenumsage_ss" = truenumsage_ss,
                    "truecatchage_ss" = truecatchage_ss,
                    "truediscage_ss" = truediscage_ss,
                    "funct.group_ss" = funct.groups_ss,
                    "biol" = omlist$biol,
                    "boxpars" = omlist$boxpars,
                    "runpar" = omlist$runpar)

  if(save){
    saveRDS(omlist_ss, file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))
  }
  if(removefullom) rm(omlist) #not removing passed data object

  return(omlist_ss)
}

Now test modified om_species in mskeyrun repo:

> NOBAom_fw <- om_species(fwspp$Name, NOBAom, save = FALSE) #save by hand, don't overwrite 11 species omlist_ss
> NOBAom_fw$truebiopool_ss
               species agecl time polygon     atoutput
1             Ben_bact     1    0       1 1.376317e+05
2        Detrivore_ben     1    0       1 1.088091e+05
3             Dinoflag     1    0       1 1.312067e+04
4      Gel_zooplankton     1    0       1 1.569070e+05
5              Lab_det     1    0       1 2.581300e+07
6         Large_phytop     1    0       1 3.268895e+05
7    Large_zooplankton     1    0       1 4.772586e+05
8   Medium_zooplankton     1    0       1 1.118751e+06
9             Pel_bact     1    0       1 3.268895e+04
10       Predatory_ben     1    0       1 1.088091e+05
11             Ref_det     1    0       1 8.904920e+06
12        Small_phytop     1    0       1 3.280167e+05
13   Small_zooplankton     1    0       1 9.152906e+05
14             Sponges     1    0       1 4.621541e+05
15            Ben_bact     1    1       1 8.875378e+05
16     Benthic_filterf     1    1       1 1.163897e-04
17             Carrion     1    1       1 1.700432e-10
18              Corals     1    1       1 1.148013e-04
19       Detrivore_ben     1    1       1 1.096034e+05
20            Dinoflag     1    1       1 4.509383e+01
21     Gel_zooplankton     1    1       1 1.563221e+05
22           King_crab     1    1       1 0.000000e+00
23             Lab_det     1    1       1 8.511963e+07
24        Large_phytop     1    1       1 8.999719e+03
25   Large_zooplankton     1    1       1 5.237310e+05
26  Medium_zooplankton     1    1       1 5.390765e+06
27            Pel_bact     1    1       1 1.230456e+06
28       Predatory_ben     1    1       1 1.089640e+05
29             Ref_det     1    1       1 8.142682e+06
30        Small_phytop     1    1       1 1.177500e+03
31   Small_zooplankton     1    1       1 7.373987e+07
32             Sponges     1    1       1 4.647222e+05
33            Ben_bact     1    2       1 3.080604e+06
34     Benthic_filterf     1    2       1 1.199844e-04
35             Carrion     1    2       1 1.698761e-10
36              Corals     1    2       1 1.168009e-04

Now try a new om_index that rbinds the truebio and truebiopool objects together for a biomass based survey

#'Generate index data from atlantisom
#'
#'#'@description A wrapper function to create survey and fishery index data for assessment input.
#'Takes the output of \code{om_species}. Wrapper can generate replicates. Saves output as .rds
#'Results for more than one survey are generated with multiple survey config files and
#'saved as separate .rds files.
#'@param usersurvey survey config file in format of /config/usersurvey.R
#'@param userfishery fishery config file in format of /config/fisherycensus.R
#'@param omlist_ss output of \code{om_species}
#'@param n_reps number of replicate indices to be generated
#'@template save
#'@return Returns list objects containing dataframes of survey biomass index and total catch:
#' \itemize{
#'  \item{survObsBiomB, list of replicate dataframes of observed survey biomass (tons)}
#'  \item{fishObsCatchB, list of replicate dataframes of observed fishery catch (tons)}
#' },
#'
#'@export
#'
#'@family wrapper functions
#'@author Sarah Gaichas
#'
#'@examples
#'\dontrun{
#' # assuming CC3om is output of om_init(here("config/CC3config.r"))
#' # and CC3om_sardine <- om_species(c("Pacific_sardine"), CC3om)
#'
#'CC3om_sard_ind <- om_index(usersurvey = here("config/usersurvey.R"),
#'     userfishery = here("config/fisherycensus.R"),
#'     omlist_ss = CC3om_sardine,
#'     n_reps = 5,
#'     save = TRUE)
#'}
#'
om_index <- function(usersurvey = usersurvey_file,
                     userfishery = userfishery_file,
                     omlist_ss,
                     n_reps = n_reps,
                     save = TRUE){

  #one script for dimension parameters to be used in multiple functions
  source("config/omdimensions.R", local = TRUE)

  # user options for survey--default is a census with mid-year sample
  # allows muliple surveys
  survObsBiomBs <- list()

  for (s in usersurvey)
  {
    source(s, local = TRUE)

    # add in biomass pools if they exist
    if("truebiopool_ss" %in% names(omlist_ss)){
      truebio_ss <- rbind(omlist_ss$truebio_ss, omlist_ss$truebiopool_ss)
    }else{
      truebio_ss <- omlist_ss$truebio_ss
    }

    #biomass based fishery independent survey index
    # this uses result$biomass_ages to sample biomass directly, no need for wt@age est
    survey_B <- atlantisom::create_survey(dat = truebio_ss,
                                          time = survtime,
                                          species = survspp,
                                          boxes = survboxes,
                                          effic = surveffic,
                                          selex = survselex)

    # call sample_survey_biomass with a bunch of 1000s for weight at age
    # in the code it multiplies atoutput by wtatage/1000 so this allows us to use
    # biomass directly
    wtage <- data.frame(species=rep(names(age_classes), n_age_classes),
                        agecl=unlist(sapply(n_age_classes,seq)),
                        wtAtAge=rep(1000.0,sum(n_age_classes)))

    # this is the step to repeat n_reps time if we want different realizations
    # of the same survey design specified above; only observation error differs
    # using the census cv of 0 will produce identical reps!
    survObsBiomB <- list()
    for(i in 1:n_reps){
      survObsBiomB[[i]] <- atlantisom::sample_survey_biomass(survey_B, surv_cv, wtage)
    }

    #save survey indices, takes a long time to generate with lots of reps/species
    if(save){
      saveRDS(survObsBiomB, file.path(d.name, paste0(scenario.name, "_",
                                                     survey.name, "surveyB.rds")))
    }

    survObsBiomBs[[survey.name]] <- survObsBiomB
  }

  #configure the fishery, a default is in config/fisherycensus.R
  #fishery configuration can specify only area and time of observation
  #fishery species inherited from omlist_ss
  #this is total catch not by fleet, so only one "fishery"

  # 2023 update, can now get fleet specific catch by polygon
  # user options for fishery--default is a census in all areas for all fleets
  # allows multiple fisheries
  fishObsCatchBs <- list()

  for(f in userfishery){

    source(f, local = TRUE)

    # 2023 update: we can now subset fishery catch from CATCH.nc using fleet output
    # create_fishery_subset as currently written aggregates across fleets and polygons
    # change so that fleets to be aggregated together are specified in the userfishery file

    fishery_C <- atlantisom::create_fishery_subset(dat = omlist_ss$truecatchtons_ss,
                                                   time = fishtime,
                                                   fleets = fishfleets,
                                                   species = fishspp,
                                                   boxes = fishboxes)

    fishObsCatchB <- list()
    for(i in 1:n_reps){
      fishObsCatchB[[i]] <- atlantisom::sample_fishery_totcatch(fishery_C, fish_cv)
    }

    if(save){
      saveRDS(fishObsCatchB, file.path(d.name, paste0(scenario.name, "_",
                                                      fishery.name, "fishCatch.rds")))
    }

    fishObsCatchBs[[fishery.name]] <- fishObsCatchB

  }

  indices <- list("survObsBiomB" = survObsBiomBs,
                  "fishObsCatchB" = fishObsCatchBs)

  return(indices)
}

It works as tested from the mskeyrun SimData.Rmd (both age structured and biopool groups are in the output):

> NOBAom_fw_ind$survObsBiomB$BTS_fall_allbox_effic1_fw[[1]] |> dplyr::group_by(species) |> dplyr::summarise(survBtot = sum(atoutput))
# A tibble: 40 × 2
   species         survBtot
   <chr>              <dbl>
 1 Beard_seal      8.14e+ 4
 2 Ben_bact        9.36e+12
 3 Benthic_filterf 7.86e+10
 4 Carrion         7.46e- 7
 5 Corals          3.19e+11
 6 Demersal_large  6.47e+ 6
 7 Demersals_other 1.81e+ 8
 8 Detrivore_ben   2.59e+10
 9 Dinoflag        1.31e+ 6
10 Fin_whale       3.87e+ 7
# … with 30 more rows
# ℹ Use `print(n = ...)` to see more rows

Plots of the new mskeyrun::simSurveyIndexFW object are here.