In r4atlantis/atlantisom: Generates Data Sets From Atlantis Scenarios

knitr::opts_chunk$set(echo = TRUE)

Setup

First, you will want to load packages and initialize the correct file and working directory names.

library(tidyr)
library(dplyr)
library(ggplot2)
library(data.table)
library(here)

If you want to load the local of the atlantisom R package, use devtools:

#library(devtools)
#package.dir <- "C:/Users/chris/Documents/Github/atlantisom"
#package.dir <- here()
#devtools::load_all(package.dir)
library(atlantisom)

Or you can install directly from the Github repository.

devtools::install_github("r4atlantis\atlantisom")

Initializing input files and directories

You will first need to tell atlantisom to know where to look for the output and input files from your atlantis model run. Here, we will give the directory where the Atlantis inputs and outputs are stored d.name, the location of the functional groups file functional.group.file (.csv), the biomass pools file biomass.pools.file (.nc), the box locations file box.file (.bgm), an initial.conditions.file (.nc), the biology .prm file biol.prm.file (.prm), and the run .prm file run.prm.file (.prm). You will also need to specify a scenario name, which will be used to define the output files (i.e. output is stored in a number of netCDF files of the format: output.nc). All of these files should be stored in d.name.

#ss.directory <- here("inst/extdata/Sardine_SS_files")

#Need a way to manage data that doesn't rely on specific paths
d.name <- here("non_synced_data/CalCurrentSummitScenario1")
source(here("config/CCConfig.R"))

Getting the "true" operating model values

There are a number of functions in the package that begin with the prefix load that load various files. See documentation if you'd only like to load one file. The atlantisom::run_truth() function uses the above file definitions and calls a number of the load functions to read in all of the atlantis output. Note: this call reads in a number of large .nc files, so it will take a few minutes to return.

#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

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

Now the R object results with the comprehensive results from the Atlantis model has been read in.

Simulating "sampling" from an Atlantis OM

There are a number of functions to "sample" biological data from this comprehensive Atlantis output. These include create_survey, sample_ages, sample_diet, sample_fish, and sample_survey_biomass. Below, we use the sample_survey_biomass function to extract time series of a biomass index for a survey and fisheries CPUE for Pacific sardine.

species=c("Pacific_sardine")

####### Boxes = spatial area, use all for now
# effic - ?
# sel - age-based selectivity curve
boxes <- unique(results$nums$polygon)
effic <- data.frame(species=species, efficiency=0.5)
sel<- data.frame(species=rep(species,10),
               agecl=c(1:10),
              selex=rep(1,10))

#Extract survey numbers
tmp <- create_survey(dat=results$nums, time=timeall, species=species, boxes=boxes, effic=effic, selex=sel)

#wtAtAge comes from empirical data from assessment, in kg
wtAtAge <- data.frame(species=rep("Pacific_sardine",10),
                agecl=1:10,
                wtAtAge=c( 0.0542, 0.0837,  0.1103, 0.1323, 0.1497, 0.163,  0.1729, 0.1801, 0.1854, 0.1941))

#Set survey CV at 0.2
cv <- data.frame(species="Pacific_sardine", cv=c(0.2))

#Sample the survey biomass
survObsBiom <- sample_survey_biomass(dat=tmp,cv=cv,wtAtAge)

#Define selectivity curve for the fishery
selCatch<- data.frame(species=rep(species,10),
                 agecl=c(1:10),
                 selex=(1/(1+exp(-0.5*(1:10-2)))))

#How do we distinguish survey sampling from fisheries sampling? Maybe define by historical sampling locations?
fishery <- create_survey(dat=results$nums, time=seq(1,99,1), species=species, boxes=boxes, effic=effic, selex=selCatch)

# Make CV of fishery CPUE slightly larger
cv <- data.frame(species="Pacific_sardine", cv=c(0.4))

# Get catch biomass
catchObsBiom <- sample_survey_biomass(dat=fishery,cv=cv,wtAtAge)

We now have the full time series of catch and survey biomass from the atlantis output for Pacific sardine. Next, we read in the dummy SS data (.dat) file and replace the abundance indices with the sampled time series.

#Read in dummy data
file <- system.file("extdata","Sardine_SS_files",file="sardEM_3_3.dat",package="atlantisom")
sardine.dat <- r4ss::SS_readdat_3.30(file)

#Replace CPUE observations for first time series with output from Atlantis - last 25 years
# SS units are metric tons, so we divide by 1000
sardine.dat$CPUE[1:75,"obs"] <- survObsBiom$atoutput[25:99]/1000

#Replace CPUE observations for second time series with output from Atlantis - last 25 years
# SS units are metric tons, so we divide by 1000
sardine.dat$CPUE[76:150,"obs"] <- catchObsBiom$atoutput[25:99]/1000
r4ss::SS_writedat_3.30(datlist = sardine.dat,
                       outfile = file.path(assessment_dir, "sardEM_3_3.dat"))