knitr::opts_chunk$set(echo = TRUE)
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")
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: outputd.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"))
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.
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"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.