atlantisom demo notebook

require(tidyr)
require(dplyr)
library(ggplot2)
library(data.table)
library(here)
library(ggforce)
library(ggthemes)
devtools::install_github("r4atlantis/atlantisom")

These packages have more recent versions available.
Which would you like to update?

 1:   All
 2:   CRAN packages only
 3:   None
 4:   chron     (2.3-53  -> 2.3-54) [CRAN]
 5:   cli       (1.1.0   -> 2.0.0 ) [CRAN]
 6:   digest    (0.6.20  -> 0.6.23) [CRAN]
 7:   ellipsis  (0.2.0.1 -> 0.3.0 ) [CRAN]
 8:   pkgconfig (2.0.2   -> 2.0.3 ) [CRAN]
 9:   purrr     (0.3.2   -> 0.3.3 ) [CRAN]
10:   R6        (2.4.0   -> 2.4.1 ) [CRAN]
11:   Rcpp      (1.0.1   -> 1.0.3 ) [CRAN]
12:   rlang     (0.4.0   -> 0.4.2 ) [CRAN]
13:   RNetCDF   (1.9-1   -> 2.1-1 ) [CRAN]
14:   xml2      (1.2.0   -> 1.2.2 ) [CRAN]
Enter one or more numbers separated by spaces, or an empty line to cancel

trying URL 'https://cran.cnr.berkeley.edu/bin/windows/contrib/3.5/tidyr_1.0.0.zip'
Content type 'application/zip' length 1292969 bytes (1.2 MB)
downloaded 1.2 MB

package ‘tidyr’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\chris\AppData\Local\Temp\RtmpSw1oMF\downloaded_packages
  
  
  
√  checking for file 'C:\Users\chris\AppData\Local\Temp\RtmpSw1oMF\remotes27e422e0570b\r4atlantis-atlantisom-b43c34e/DESCRIPTION' (380ms)

  
  
  
-  preparing 'atlantisom': (1.1s)
   checking DESCRIPTION meta-information ...
  
√  checking DESCRIPTION meta-information

  
  
  
-  checking for LF line-endings in source and make files and shell scripts (501ms)

  
-  checking for empty or unneeded directories

  
     NB: this package now depends on R (>= 3.5.0)
     WARNING: Added dependency on R >= 3.5.0 because serialized objects in  serialize/load version 3 cannot be read in older versions of R.  File(s) containing such objects:  'atlantisom/nemow_demo/CC_2063_OA_OFF_22/CCV3catchN_census_sard_hake.rds'  WARNING: Added dependency on R >= 3.5.0 because serialized objects in  serialize/load version 3 cannot be read in older versions of R.  File(s) containing such objects:  'atlantisom/nemow_demo/CC_2063_OA_OFF_22/CCV3catchage_census_sard_hake.rds'  WARNING: Added dependency on R >= 3.5.0 because serialized objects in  serialize/load version 3 cannot be read in older versions of R.  File(s) containing such objects:  'atlantisom/nemow_demo/CC_2063_OA_OFF_22/CCV3catchlengthwt_census_sard_hake.rds'  WARNING: Added dependency on R >= 3.5.0 because serialized objects in  serialize/load version 3 cannot be read in older versions of R.  File(s) containing such objects:  'atlantisom/nemow_demo/CC_2063_OA_OFF_22/CCV3lengthwt_census_sard_hake.rds'  WARNING: Added dependency on R >= 3.5.0 because serialized objects in  serialize/load version 3 cannot be read in older versions of R.  File(s) containing such objects:  'atlantisom/nemow_demo/CC_2063_OA_OFF_22/CCV3natage_census_sard_hake.rds'  WARNING: Added dependency on R >= 3.5.0 because serialized objects in  serialize/load version 3 cannot be read in older versions of R.  File(s) containing such objects:  'atlantisom/nemow_demo/CC_2063_OA_OFF_22/CCV3surveyBcensus.rds'  WARNING: Added dependency on R >= 3.5.0 because serialized objects in  serialize/load version 3 cannot be read in older versions of R.  File(s) containing such objects:  'atlantisom/nemow_demo/CC_2063_OA_OFF_22/CCV3surveyNcensus.rds'
-  building 'atlantisom_0.1.tar.gz'

  
   

* installing *source* package 'atlantisom' ...
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
  converting help for package 'atlantisom'
    finding HTML links ... done
    SS_write_biol                           html  
    SS_write_comps                          html  
    SS_write_ts                             html  
    aggregateData                           html  
    aggregateDensityData                    html  
    aggregate_comps                         html  
    calc_Z                                  html  
    calc_age2length                         html  
    calc_biomass_age                        html  
    calc_pred_diet                          html  
    calc_stage2age                          html  
    calc_timestep2time                      html  
    create_env                              html  
    create_fishery_subset                   html  
    finding level-2 HTML links ... done

Rd warning: C:/Users/chris/AppData/Local/Temp/RtmpCOxRni/R.INSTALL660541e4c54/atlantisom/man/create_fishery_subset.Rd:20: missing link 'create_fishery_data'
    create_survey                           html  
Rd warning: C:/Users/chris/AppData/Local/Temp/RtmpCOxRni/R.INSTALL660541e4c54/atlantisom/man/create_survey.Rd:20: missing link 'create_fishery_data'
    get_boundary                            html  
    load_biolprm                            html  
    load_box                                html  
    load_boxarea                            html  
    load_bps                                html  
    load_catch                              html  
    load_diet_comp                          html  
    load_fgs                                html  
    load_meta                               html  
    load_nc                                 html  
    load_nc_physics                         html  
    load_runprm                             html  
    load_yoy                                html  
    reformat_compositions                   html  
Rd warning: C:/Users/chris/AppData/Local/Temp/RtmpCOxRni/R.INSTALL660541e4c54/atlantisom/man/reformat_compositions.Rd:22: missing link 'create_fishery_data'
    run_stocksynthesis                      html  
    run_truth                               html  
    sample_ages                             html  
Rd warning: C:/Users/chris/AppData/Local/Temp/RtmpCOxRni/R.INSTALL660541e4c54/atlantisom/man/sample_ages.Rd:20: missing link 'create_fishery_data'
    sample_diet                             html  
    sample_fish                             html  
Rd warning: C:/Users/chris/AppData/Local/Temp/RtmpCOxRni/R.INSTALL660541e4c54/atlantisom/man/sample_fish.Rd:20: missing link 'create_fishery_data'
    sample_survey_biomass                   html  
Rd warning: C:/Users/chris/AppData/Local/Temp/RtmpCOxRni/R.INSTALL660541e4c54/atlantisom/man/sample_survey_biomass.Rd:20: missing link 'create_fishery_data'
    sample_survey_numbers                   html  
Rd warning: C:/Users/chris/AppData/Local/Temp/RtmpCOxRni/R.INSTALL660541e4c54/atlantisom/man/sample_survey_numbers.Rd:20: missing link 'create_fishery_data'
    vbgf_func                               html  
    write_meta                              html  
** building package indices
** installing vignettes
** testing if installed package can be loaded
*** arch - i386
*** arch - x64
* DONE (atlantisom)
In R CMD INSTALL

library(atlantisom)

Introduction

This demonstrates the Atlantis to atlantisom to SS3 workflow with the California Current run that includes climate forcing, recruitment variability, and the fishing scenario outlined for our project in Norway: “output_CC_2063_OA_OFF_22”.

We will use Atlantis model output to generate time series and composition data for a species “sardine”, that doesn’t require any age class splitting.

species_ss <- c("Pacific_sardine")
age_classes <- 1:10
source(here("config/NEMOWConfig.R"))
needed.files <- ls()
print(unlist(lapply(needed.files, FUN=get)))

 [1] "1"                                                                      
 [2] "2"                                                                      
 [3] "3"                                                                      
 [4] "4"                                                                      
 [5] "5"                                                                      
 [6] "6"                                                                      
 [7] "7"                                                                      
 [8] "8"                                                                      
 [9] "9"                                                                      
[10] "10"                                                                     
[11] "outputCCV3BiomIndx.txt"                                                 
[12] "CalCurrentV3_Biol.prm"                                                  
[13] "DIVCalCurrentV3_Biol.nc"                                                
[14] "CalCurrentV3_utm.bgm"                                                   
[15] "outputCCV3Catch.txt"                                                    
[16] "C:/Users/chris/Documents/GitHub/atlantisom/nemow_demo/CC_2063_OA_OFF_22"
[17] "CalCurrentV3Groups.csv"                                                 
[18] "DIVCalCurrentV3_Biol.nc"                                                
[19] "TRUE"                                                                   
[20] "CalCurrentV3_run.xml"                                                   
[21] "CCV3"                                                                   
[22] "Pacific_sardine"                                                        

The below example requires a folder in your project directory named “atlantisom” containing the above files from an Atlantis output.

#Check for needed files
#Load functional groups
funct.groups <- load_fgs(dir=d.name,
                         file_fgs = functional.groups.file)
#Get just the names of active functional groups
funct.group.names <- funct.groups %>% 
  filter(IsTurnedOn == 1) %>%
  select(Name) %>%
  .$Name

This can be run with species_ss in select_groups to get results only for that species. WARNING: running with all species takes a long time! (50 minutes for this output)

  #Store all loaded results into an R object
sardine_truth <- if(!file.exists(file.path(d.name, 
                          "sardine_truth.RData"))){
  #Store all loaded results into an R object
  sardine_truth <- run_truth(scenario = scenario.name,
                     dir = d.name,
                     file_fgs = functional.groups.file,
                     file_bgm = box.file,
                     select_groups = "Pacific_sardine",
                     file_init = initial.conditions.file,
                     file_biolprm = biol.prm.file,
                     file_runprm = run.prm.file,
                     verbose = TRUE
  )
} else {
  sardine_truth <- get(load(file.path(d.name,
                              "sardine_truth.RData")))
}

This truth object shows the true population dynamics for sardine from Atlantis which we need for reference to compare with estimated population dynamics from the stock assessment estimation models.

Then we will have a function that implements the “survey” that produces the stock assessment inputs and sends those to the SS file writing functions. Here the survey functions are still separate.

Derive “Data” from Atlantis to give to SS3

Specify our survey sampling. This could come from other information such as the overlap of actual survey stations with OM polygons, experiments evaluating survey selectivity and efficiency, actual sample-based survey cv, etc.

Here are specifications which approximate the true sardine acoustic survey

# 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(0,stepperyr)))
# a survey that takes place once per year mid year
annualmidyear <- seq(midptyr, noutsteps, stepperyr)
# learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(runpar$outputstepunit=="days") 365/runpar$toutfinc
# should return all model areas
boxpars <- load_box(d.name, box.file)
boxall <- c(0:(boxpars$nbox - 1))
source(here("config/sardine_survey.R"))

Get input survey biomass for SS (these have q, selectivity, cv)

#Sample survey - 3rd timestep simulates a summer survey
survey_out <- create_survey(dat=sardine_truth$nums, 
                            time=survey_sample_full, 
                            species=species, 
                            boxes=boxall, 
                            effic=effic, 
                            selex=sel)
#Try a biomass based survey for comparison
survey_outB <- create_survey(dat=sardine_truth$biomass_ages,
                            time=survey_sample_full,
                            species=species,
                            boxes=boxall,
                            effic=effic,
                            selex=sel)
#Set effective sample size for age compositions
effN <- surveyEffN
highEffN <- data.frame(species=species, effN=rep(effN, length(species)))
#Sample fish for age composition
age_comp_data <- sample_fish(survey_out, highEffN)
# aggregate true resn per survey design
survey_aggresnstd <- aggregateDensityData(dat = sardine_truth$resn,
                                          time = survey_sample_full,
                                          species = species,
                                          boxes = boxall)
# aggregate true structn per survey design
survey_aggstructnstd <- aggregateDensityData(dat =sardine_truth$structn,
                                             time = survey_sample_full,
                                             species = species,
                                             boxes = boxall)
ss_structnstd <- sample_fish(survey_aggstructnstd,
                             effN,
                             sample=FALSE)
ss_resnstd <- sample_fish(survey_aggresnstd,
                          effN,
                          sample=FALSE)
#Extract length composition data
ss_length_stdsurv <- calc_age2length(structn = ss_structnstd,
                                     resn = ss_resnstd,
                                     nums = age_comp_data,
                                     biolprm = sardine_truth$biolprm, fgs = sardine_truth$fgs,
                                     CVlenage = CVs$lenage, remove.zeroes=TRUE)
#Need to replace with interp function
wtAtAge <- ss_length_stdsurv$muweight %>%
  select(species, agecl, time, wtAtAge = atoutput) %>%
  mutate(wtAtAge = wtAtAge/1000)
# CV function
cv <- data.frame(species=species_ss, cv=CVs$survey)
#Sample survey biomass
survObsBiom <- sample_survey_biomass(dat=survey_out,cv=cv,wtAtAge)
# check against survey with truth$biomass_ages output
wtage1 <- data.frame(species=rep(species_ss, each=max(age_classes)),
                    agecl=rep(c(age_classes),length(species_ss)),
                    wtAtAge=rep(1000.0,length(species_ss)*max(age_classes)))
survObsBiomB <- sample_survey_biomass(dat=survey_outB,cv=cv,wtage1)
# survey numbers, not sure which SS needs?
survObsNum <- sample_survey_numbers(dat=survey_out,cv=cv)

Get composition inputs for SS (survey and fishery catch at age, survey and fishery lengths, survey and fishey weight at age).

Because catch composition information goes in to the assessment as a proportion, we can use fishery catch at age from this legacy codebase even with absolute catch numbers likely half what they should be overall.

# Survey length comps and wtage done above to get survey ts using nums*wtage approach
#We end up using CAAL for the survey below, so let's generate fishery age comps instead
effN <- fisheryEffN/fstepperyr
effN <- data.frame(species=species, effN=effN)
#catch at age each timestep summed over polygons
# catch at age by area and timestep
catch_numbers <-  create_fishery_subset(dat = sardine_truth$catch,
                                         time = fish_times,
                                         species = species,
                                         boxes = boxall)
catch_numsss_samp <- sample_fish(catch_numbers, effN)
rm(catch_numbers)
gc()

           used  (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells  1182750  63.2   20832527 1112.6  40688531 2173.1
Vcells 46855901 357.5  128263088  978.6 250506075 1911.3

#Get weights
# aggregate true resn per fishery subset design
catch_aggresnss <- aggregateDensityData(dat = sardine_truth$resn,
                                 time = fish_times,
                                 species = species,
                                 boxes = boxall)
# aggregate true structn fishery subsetdesign
catch_aggstructnss <- aggregateDensityData(dat = sardine_truth$structn,
                                 time = fish_times,
                                 species = species,
                                 boxes = boxall)
#dont sample these, just aggregate them using median
catch_structnss_samp <- sample_fish(catch_aggstructnss, effN, sample = FALSE)
catch_resnss_samp <-  sample_fish(catch_aggresnss, effN, sample = FALSE)
# these fishery lengths and weight at age are each output timestep
catch_lengthwt_samp <- calc_age2length(structn = catch_structnss_samp,
                                 resn = catch_resnss_samp,
                                 nums = catch_numsss_samp,
                                 biolprm = sardine_truth$biolprm, 
                                 fgs = sardine_truth$fgs,
                                 maxbin = maxbin,
                                 CVlenage = CVs$lenage, remove.zeroes=TRUE)