Atlantis to SS for California Current sardines and hake: full fishing and climate scenario
Sarah Gaichas and Christine Stawitz
01 October, 2019
library(tidyr)
require(dplyr)
library(ggplot2)
library(data.table)
library(here)
library(ggforce)
library(ggthemes)
library(atlantisom)Introduction
Now we test the atlantisom to SS workflow with a new 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 go from atlantis model output to SS input in this page for a single species, “sardine”, that doesn’t require any age class splitting.
As of 31 May 2019, the run_truth function has been modified to do the catch numbers fix for legacy models.
initCCA <- TRUE
initNEUS <- FALSE
initNOBA <- FALSE
species_ss <- c("Pacific_sardine", "Mesopel_M_Fish")
#function to make a config file? need one for each atlantis run
if(initCCA) source(here("config/CC3Config.R"))
if(initNEUS) source(here("config/NEUSConfig.R"))
if(initNOBA) source(here("config/NOBAConfig.R"))#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) %>%
.$NameThis can be run with species_ss in select_groups to get results only for that species. Here I’m running with all species; took about 50 minutes to return.
# default run_truth setup will save the file, so check for that first
if(!file.exists(file.path(d.name,
paste0("output", scenario.name, "run_truth.RData")))){
#Store all loaded results into an R object
truth <- 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,
verbose = TRUE
)
} else {
truth <- get(load(file.path(d.name,
paste0("output", scenario.name, "run_truth.RData"))))
}In the big picture, we will have the true population dynamics from Atlantis which we need for reference to compare with estimated population dynamics from the stock assessment estimation models.
Ideally we will have a function that gets the whole truth (and nothing but the truth) from Atlantis and stores it in a sensible way for later comparison and calculation of performance metrics.
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.
Truth from OM
The following settings should achieve a survey that samples all Atlantis model output timesteps, all species, and all model polygons, with perfect efficiency and full selectivity for all ages:
# make a function for this
source(here("config/census_spec.R"))Get true total biomass, total numbers. This code gets results for all species and saves the output true B and N.
# this uses result$biomass_ages to sample biomass directly, no need for wt@age est
survey_testBall <- create_survey(dat = truth$biomass_ages,
time = timeall,
species = survspp,
boxes = boxall,
effic = effic1,
selex = selex1)
# make up a constant 0 cv for testing
surv_cv <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))
# 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(survspp, each=10),
agecl=rep(c(1:10),length(survspp)),
wtAtAge=rep(1000.0,length(survspp)*10))
surveyB_frombio <- sample_survey_biomass(survey_testBall, surv_cv, wtage)
#save for later use, takes a long time to generate
saveRDS(surveyB_frombio, file.path(d.name, paste0(scenario.name, "surveyBcensus.rds")))
# this uses result$nums and a new function to get survey index in numbers (abundance)
survey_testNall <- create_survey(dat = truth$nums,
time = timeall,
species = survspp,
boxes = boxall,
effic = effic1,
selex = selex1)
# as above, make up a constant 0 cv for testing
surv_cv <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))
surveyN <- sample_survey_numbers(survey_testNall, surv_cv)
#save for later use, takes a long time to generate
saveRDS(surveyN, file.path(d.name, paste0(scenario.name, "surveyNcensus.rds")))Plot true biomass and numbers for both “hake” and “sardines.”
trueB <- readRDS(file.path(d.name, paste0(scenario.name,"surveyBcensus.rds")))
trueB_ss <- trueB %>%
filter(species==species_ss)
# read Atlantis output files
atBtxt2 <- read.table(file.path(d.name, paste0("output", scenario.name, "BiomIndx.txt")), header=T)
groupslookup <- funct.groups %>%
filter(IsTurnedOn > 0)
atBtxt2tidy <- atBtxt2 %>%
select(Time:DIN) %>%
#select(Time, FPL:DIN) %>%
rename_(.dots=with(groupslookup, setNames(as.list(as.character(Code)), Name))) %>%
gather(species, biomass, -Time) %>%
filter(species %in% c("Mesopel_M_Fish", "Pacific_sardine"))
plotB <-ggplot() +
geom_line(data=trueB_ss, aes(x=time/stepperyr,y=atoutput, color="survey census B"),
alpha = 10/10) +
geom_point(data=atBtxt2tidy, aes(x=Time/365,y=biomass, color="txt output true B"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
plotB +
facet_wrap(~species, scales="free") 
trueN <- readRDS(file.path(d.name, paste0(scenario.name,"surveyNcensus.rds")))
trueN_ss <- trueN %>%
filter(species==species_ss)
plotN <-ggplot() +
geom_line(data=trueN_ss, aes(x=time/stepperyr,y=atoutput, color="survey census N"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
plotN +
facet_wrap(~species, scales="free") 
Get true population age comp, length comp, and weight at age. Now we go to single species to speed things up.
# wrap these individual calls into a function
survey_N_ss <- create_survey(dat = truth$nums,
time = timeall,
species = species_ss,
boxes = boxall,
effic = effic1,
selex = selex1)
# apply default sample fish to get numbers, effNhigh samples all, see census_spec
# NOPE! effNhigh is not big enough for census of sardines, bigger breaks function
#numssshigh <- sample_fish(survey_N_ss, effNhigh)
numsss <- aggregate(survey_N_ss$atoutput,list(survey_N_ss$species,survey_N_ss$agecl,survey_N_ss$time),sum)
names(numsss) <- c("species","agecl","time","atoutput")
numssshigh <- data.frame(species = numsss$species,
agecl = numsss$agecl,
polygon = NA,
layer = NA,
time = numsss$time,
atoutput = numsss$atoutput)
# save age comps
saveRDS(numssshigh, file.path(d.name, paste0(scenario.name, "natage_census_sard_hake.rds")))
# aggregate true resn per survey design
aggresnss <- aggregateDensityData(dat = truth$resn,
time = timeall,
species = species_ss,
boxes = boxall)
# aggregate true structn per survey design
aggstructnss <- aggregateDensityData(dat = truth$structn,
time = timeall,
species = species_ss,
boxes = boxall)
#dont sample these, just aggregate them using median
structnss <- sample_fish(aggstructnss, effNhigh, sample = FALSE)
resnss <- sample_fish(aggresnss, effNhigh, sample = FALSE)
# this function does the length sampling and also weight at age
length_census_ss <- calc_age2length(structn = structnss,
resn = resnss,
nums = numssshigh,
biolprm = truth$biolprm, fgs = truth$fgs,
CVlenage = 0.1, remove.zeroes=TRUE)
#save for later use, takes a long time to generate
saveRDS(length_census_ss, file.path(d.name, paste0(scenario.name, "lengthwt_census_sard_hake.rds")))Plot some of the age, length comps, and wtatage for selected timesteps (after burnin)
Natage <- readRDS(file.path(d.name, paste0(scenario.name, "natage_census_sard_hake.rds")))
Natagesub <- Natage %>%
filter(time %in% c(155:195))
for(sp in unique(Natage$species)){
Natagesubsp <- Natagesub %>%
filter(species==sp)
Natageplot <- ggplot(Natagesubsp, aes(x=agecl, y=atoutput)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste(scenario.name,
sp))
Natageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 1, scales="free")
Natageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 2, scales="free")
Natageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 3, scales="free")
Natageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 4, scales="free")
}lengthwt_census_ss <- readRDS(file.path(d.name, paste0(scenario.name, "lengthwt_census_sard_hake.rds")))
len <- lengthwt_census_ss$natlength
lensub <- len %>%
filter(time %in% c(155:195))
for(sp in unique(len$species)){
lensubsp <- lensub %>%
filter(species==sp)
lfplot <- ggplot(len, aes(upper.bins)) +
geom_bar(aes(weight = atoutput)) +
theme_tufte() +
labs(subtitle = paste(scenario.name,
len$species))
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 1, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 2, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 3, scales="free_y")
}wtage_ann <- lengthwt_census_ss$muweight %>%
filter(time %in% annualmidyear)
# reverse to show agecl time series of wt
wageplot <- ggplot(wtage_ann, aes(time, atoutput)) +
geom_line(aes(colour = factor(agecl))) +
theme_tufte() +
theme(legend.position = "bottom") +
xlab("model timestep (5 per year)") +
ylab("average individual weight (g)") +
labs(subtitle = paste0(scenario.name, " annual mid year sample"))
wageplot + facet_wrap(c("species"), scales="free_y")
Get true population parameters (M in progress, S-R, growth)
Get true catch in tons (may not be different from input for our project)
truecatchbio <- load_catch(d.name, file_catch = catch.file, fgs = funct.groups)
truecatchbio_ss <- truecatchbio[truecatchbio$species == species_ss,]
# note: time output of this file is in days
# what model timestep is this output matching?
# rescale catch in biomass time from days to years
truecatchbio_ss$time <- as.integer(truecatchbio_ss$time/365)
plotcatch <- ggplot() +
geom_line(data=truecatchbio_ss, aes(x=time,y=atoutput, color="true catch bio"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
plotcatch +
facet_wrap(~species, scales="free") 
Get true catch at age and other comps (may not be different from input for our project). First see if catch output is at same timescale as population by checking toutfinc: 5
# get survey nums with full (no) selectivity
catchN_ss <- create_fishery_subset(dat = truth$catch,
time = timeall,
species = species_ss,
boxes = boxall)
# try summing this to get true catch in nums by timestep with surv_cv = 0
totcatchN <- sample_survey_numbers(catchN_ss, surv_cv)
# save true total catch in numbers by timestep, need to sum to year
saveRDS(totcatchN, file.path(d.name, paste0(scenario.name, "catchN_census_sard_hake.rds")))
# apply default sample fish as before to get numbers
# no, dont apply sample, effN isnt high enough and can't be made high enough
#catch_numssshigh <- sample_fish(catchN_ss, effNhigh)
# just aggregte true numbers over layer and polygon
catch_numsss <- aggregate(catchN_ss$atoutput,list(catchN_ss$species,catchN_ss$agecl,catchN_ss$time),sum)
names(catch_numsss) <- c("species","agecl","time","atoutput")
catch_numssshigh <- data.frame(species = catch_numsss$species,
agecl = catch_numsss$agecl,
polygon = NA,
layer = NA,
time = catch_numsss$time,
atoutput = catch_numsss$atoutput)
# save true catch at age
saveRDS(catch_numssshigh, file.path(d.name, paste0(scenario.name, "catchage_census_sard_hake.rds")))
# aggregate true resn per survey or fishery subset design
catch_aggresnss <- aggregateDensityData(dat = truth$resn,
time = timeall,
species = species_ss,
boxes = boxall)
# aggregate true structn per survey or fishery subsetdesign
catch_aggstructnss <- aggregateDensityData(dat = truth$structn,
time = timeall,
species = species_ss,
boxes = boxall)
#dont sample these, just aggregate them using median
catch_structnss <- sample_fish(catch_aggstructnss, effNhigh, sample = FALSE)
catch_resnss <- sample_fish(catch_aggresnss, effNhigh, sample = FALSE)
# get catch lengths and wtage at every timestep
catch_length_census_ss <- calc_age2length(structn = catch_structnss,
resn = catch_resnss,
nums = catch_numssshigh,
biolprm = truth$biolprm, fgs = truth$fgs,
maxbin = 150,
CVlenage = 0.1, remove.zeroes=TRUE)
#save for later use, takes a long time to generate
saveRDS(catch_length_census_ss, file.path(d.name, paste0(scenario.name, "catchlengthwt_census_sard_hake.rds")))This is true catch at age for sardine:
catch_numssshigh <- readRDS(file.path(d.name, paste0(scenario.name, "catchage_census_sard_hake.rds")))
catchage <- catch_numssshigh %>%
filter(species == "Pacific_sardine")
catageplot <- ggplot(catchage, aes(x=agecl, y=atoutput)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste(scenario.name,
catchage$species))
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 1, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 2, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 3, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 4, scales="free")
Do true catch numbers times weight at age match catch in biomass? It was about 16 times too high relative to catch.txt output, and shifted one year back (which is due to the floor function, and an easy fix).
June 6: just corrected the catch nums correction in run_truth and I hope it matches now.
totcatchN <- readRDS(file.path(d.name, paste0(scenario.name, "catchN_census_sard_hake.rds")))
totcatchann <- totcatchN %>%
mutate(yr = floor(time/fstepperyr)+1) %>%
group_by(species, yr) %>%
summarise(anncatchage = sum(atoutput))
plottotC <- ggplot(totcatchann, aes(x=yr, y=anncatchage)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste("Catch in N ", scenario.name,
totcatchann$species))
plottotC
catch_numssshigh <- readRDS(file.path(d.name, paste0(scenario.name, "catchage_census_sard_hake.rds")))
catch_lengthwt_census_ss <- readRDS(file.path(d.name, paste0(scenario.name, "catchlengthwt_census_sard_hake.rds")))
CwtAtAge <- catch_lengthwt_census_ss$muweight %>%
select(species, agecl, time, wtAtAge = atoutput) %>%
mutate(wtAtAge = wtAtAge/1000)
catchbioN_ss <- merge(catch_numssshigh,CwtAtAge,by=c("species","agecl", "time"),all.x=T) %>%
mutate(biomass = atoutput*wtAtAge/1000) %>%
mutate(yr = floor(time/fstepperyr)+1) %>%
group_by(species, yr) %>%
summarise(anncatchB = sum(biomass))
plotcatch <- ggplot() +
geom_line(data=truecatchbio_ss, aes(x=time,y=atoutput, color="true catch bio from txt"),
alpha = 10/10) +
geom_point(data=catchbioN_ss, aes(x=yr,y=anncatchB, color="true catch bio from N"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
plotcatch +
facet_wrap(~species, scales="free") 
Sadly, still a problem with the catch correction if it should match catch.txt when converted to biomass. We are now off by about a factor of 2 (catch nums * wtAtAge is half true catch bio from the txt file), but this differs slightly by species:
sardcatchB <- truecatchbio_ss #%>%
# now done above
# mutate(yr = as.integer(time/365))
sardcatchBcomp <- merge(sardcatchB, catchbioN_ss)
sardcatchBcomp <- mutate(sardcatchBcomp, mult=atoutput/anncatchB)
#mean(sardcatchBcomp$mult)
#summary(sardcatchBcomp$mult)
Filter(Negate(is.null), tapply(sardcatchBcomp$mult, sardcatchBcomp$species, summary))## $Mesopel_M_Fish
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 2.305 5.475 3.217 631.051
##
## $Pacific_sardine
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.184 1.544 16.645 2.806 5382.851“Data” from OM for EMs
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.
Now we are going single species. First we generate the inputs for sardine.
Here we are using Christine’s specifications, which approximate the acoustic survey index used in the sardine assessment:
# from atlantisom/config/sardine_config.R, to be generalized here
#Set species name
species <- "Pacific_sardine" #change to species_ss if we want both
#Directory with SS files
model_dir <- "Sardine_SS_files/"
#Name of SS data file
datfile_name <- "sardEM_3_3.dat"
#Name of atlantis model output - this is an RData object from the runtruth() function
#atlantis_output_name <- "outputCCV3run_truth.RData"
#Efficiency
eff_case <- 0.5
#Age classes of your species
age_classes <- 1:10
#Selectivity function - a vector equivalent to length (age classes)
sel_by_age <- rep(1,length(age_classes))
#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5
#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May,
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)
survey_sample_time <- midptyr #3 this is more late summer-fall, changed to 2
#The last timestep to sample
total_sample <- noutsteps-1 #495
#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
total_sample, by=timestep)
#Effective sample size for composition data (annual total samples)
surveyEffN <- 1000
fisheryEffN <- 1000
#CVs for length at age, catch, and survey
CVs <- list("lenage"=0.1, "fishery"=0.01, "survey"=0.1)
#Number of years of data to pull
nyears <- 50
#Atlantis initialization period in years
burnin <- 30
#Maximum size of a fish in cm
maxbin <- 150
#Years to fish, assuming fishing output is annual and years index 1,2,3 etc
#fish_years <- burnin:(burnin+nyears-1)
# fishery output timestep (fstepperyr in census_spec.R) for this run is 5 per year
# here we are summing the catch in numbers for input in to SS.
# To use the total catch biomass output from catch.txt the above would work
# because the output is annual but in units of DAYS so needs time/365
#The last timestep to sample
total_sample <- noutsteps-1 #495
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample)
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[5], max(fish_times), by=fstepperyr) #last timestep
#fish_years <- unique(floor(fish_times/fstepperyr)+1) # my original
fish_years <- unique(floor(fish_times/fstepperyr)) #from Christine's new sardine_config.R
# potential problem: while they are from the same model year,
# fish_years dont match survey_years timesteps (survey mid year, fishery end)
# so when matching use the index, not the value itself
#Years to survey, assuming survey is once per year
#survey_years <- survey_sample_full[(burnin+1):(burnin+nyears)] # my original
survey_years <- survey_sample_full[burnin:(burnin+nyears-1)] #from Christine's new sardine_config.R
#Month of survey/fishing
survey_month <- 7
fishing_month <- 1
#Efficiency of the survey
effic <- data.frame(species=species, efficiency=eff_case)
# Assume uniform selectivity, same as selex1 here
sel<- data.frame(species=rep(species,length(age_classes)),
agecl=age_classes,
selex=sel_by_age)Get input survey biomass for SS (these have q, selectivity, cv)
#Sample survey - 3rd timestep simulates a summer survey
survey_out <- create_survey(dat=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=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 = truth$resn,
time = survey_sample_full,
species = species,
boxes = boxall)
# aggregate true structn per survey design
survey_aggstructnstd <- aggregateDensityData(dat =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 = truth$biolprm, fgs = 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=10),
agecl=rep(c(1:10),length(species_ss)),
wtAtAge=rep(1000.0,length(species_ss)*10))
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)Do some comparisons of biomass surveys with truth.
plotB <-ggplot() +
geom_point(data=atBtxt2tidy, aes(x=Time/365,y=biomass, color="txt output true B"),
alpha = 10/10) +
geom_line(data=trueB_ss, aes(x=time/stepperyr,y=atoutput, color="survey census B"),
alpha = 10/10) +
geom_point(data=survObsBiomB, aes(x=time/stepperyr,y=atoutput, color="survey sample B from biomass_ages"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
plotB +
facet_wrap(~species, scales="free") 
plotB <-ggplot() +
geom_point(data=atBtxt2tidy, aes(x=Time/365,y=biomass, color="txt output true B"),
alpha = 10/10) +
geom_line(data=trueB_ss, aes(x=time/stepperyr,y=atoutput, color="survey census B"),
alpha = 10/10) +
geom_point(data=survObsBiom, aes(x=time/stepperyr,y=atoutput, color="survey sample B from nums * wtAtAge"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
plotB +
facet_wrap(~species, scales="free") 
Get input total catch for SS (here using numbers). SS takes a fishery cv which can be applied here to get a catch time series, or we can input true catch.
PLEASE DON’T USE THIS FOR CCA. We have determined above that catch in numbers from legacy atlantis code, even corrected, is still off by a factor of 2. We leave the info here because this code can be used for atlantis codebases after December 2015.
# catch at age by area and timestep
catch_numbers <- create_fishery_subset(dat = truth$catch,
time = fish_times,
species = species,
boxes = boxall)
if (fstepperyr>1){
#need to take all output time steps and sum catch over each year
catch_numbers_area <- catch_numbers %>%
#assign a temporary year column, aggregate nums over this year
mutate(yr = floor(time/fstepperyr)+1) %>%
group_by(species, agecl, polygon, yr) %>%
summarise(anncatch_area=sum(atoutput))
# catch in numbers summed over polygon
# and agecl after time has been aggregated into 1 yr
catch_total <- catch_numbers_area %>%
group_by(yr) %>%
summarise(catch=sum(anncatch_area))
}
# catch in numbers summed over polygon and agecl if output is annual
if (fstepperyr==1){
catch_total <- catch_numbers %>%
group_by(time) %>%
summarise(catch=sum(atoutput))
}
# if we want to apply the CV to this we need to write sample_total_catch
#cv <- data.frame(species="Pacific_sardine", cv=CVs$fishery)
#catchObsNum <- sample_total_catch(dat=catch_numbers, cv=cv)Check this against true total catch numbers above. Should be the same (and it is, but both are incorrect in the absolute scale due to the CATCH.nc output problem in older codebases).
compareC <- ggplot() +
geom_line(data=totcatchann, aes(x=yr, y=anncatchage, color="census catch N"),
alpha = 10/10) +
geom_point(data=catch_total, aes(x=yr,y=catch, color="sample catch N"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
compareC +
facet_wrap(~species, scales="free") 
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_numsss_samp <- sample_fish(catch_numbers, effN)
#Get weights
# aggregate true resn per fishery subset design
catch_aggresnss <- aggregateDensityData(dat = truth$resn,
time = fish_times,
species = species,
boxes = boxall)
# aggregate true structn fishery subsetdesign
catch_aggstructnss <- aggregateDensityData(dat = 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 = truth$biolprm, fgs = truth$fgs,
maxbin = maxbin,
CVlenage = CVs$lenage, remove.zeroes=TRUE)Write data for Stock Synthesis
To read in the Stock Synthesis files, first ensure you have the latest version of r4ss. You will need to ensure the relevant Stock Synthesis files are stored in the inst/extdata/ folder, with each species example in its own folder.
#devtools::install_github("r4ss/r4ss", dependencies = FALSE)
require(r4ss)## Loading required package: r4ss#source("./config/sardine_config.R") this is above for illustration
stocksynthesis.data <- r4ss::SS_readdat_3.30(paste0("./inst/extdata/",
model_dir,
datfile_name))## Running SS_readdat_3.30## SS_readdat_3.30 - read version = 3.30## use_meanbodywt (0/1): 0## N_lbinspop:## use_lencomp (0/1): 1## N_lbins: 48## N_agebins: 10## use_MeanSize_at_Age_obs (0/1): 0## N_environ_variables: 0## Read of data file complete. Final value = 999Now changing catch units to tons!
Write time series biomass (tons) and catch (tons) data for SS using SS_write_ts:
#Test writing CPUE in biomass
stocksynthesis.data <- SS_write_ts(ss_data_list = stocksynthesis.data,
ts_data = list(survObsBiom$atoutput[fish_years],
truecatchbio_ss$atoutput[truecatchbio_ss$time %in% fish_years]),
CVs = c(CVs$survey,
CVs$fishery),
data_years = list((survObsBiom$time[fish_years]-survey_sample_time)/timestep+1, fish_years),
sampling_month = list(rep(survey_month,nyears),
rep(fishing_month,nyears)),
units = c("biomass","biomass"),
data_type=c("CPUE","catch"),
fleets = c(2,1))
stocksynthesis.data$CPUE## year seas index obs se_log
## 1 30 7 2 1202512 0.1
## 2 31 7 2 1509004 0.1
## 3 32 7 2 1131980 0.1
## 4 33 7 2 714126 0.1
## 5 34 7 2 424091 0.1
## 6 35 7 2 337944 0.1
## 7 36 7 2 201896 0.1
## 8 37 7 2 163246 0.1
## 9 38 7 2 123381 0.1
## 10 39 7 2 104026 0.1
## 11 40 7 2 146043 0.1
## 12 41 7 2 183381 0.1
## 13 42 7 2 228366 0.1
## 14 43 7 2 417989 0.1
## 15 44 7 2 465290 0.1
## 16 45 7 2 423393 0.1
## 17 46 7 2 513879 0.1
## 18 47 7 2 415702 0.1
## 19 48 7 2 481611 0.1
## 20 49 7 2 464315 0.1
## 21 50 7 2 362996 0.1
## 22 51 7 2 431447 0.1
## 23 52 7 2 339960 0.1
## 24 53 7 2 252302 0.1
## 25 54 7 2 244934 0.1
## 26 55 7 2 308850 0.1
## 27 56 7 2 423603 0.1
## 28 57 7 2 691614 0.1
## 29 58 7 2 975129 0.1
## 30 59 7 2 1067609 0.1
## 31 60 7 2 1375942 0.1
## 32 61 7 2 1354634 0.1
## 33 62 7 2 1225909 0.1
## 34 63 7 2 1651157 0.1
## 35 64 7 2 1558874 0.1
## 36 65 7 2 1021866 0.1
## 37 66 7 2 1049372 0.1
## 38 67 7 2 876197 0.1
## 39 68 7 2 1239095 0.1
## 40 69 7 2 1398204 0.1
## 41 70 7 2 1164691 0.1
## 42 71 7 2 1114676 0.1
## 43 72 7 2 1009397 0.1
## 44 73 7 2 1196984 0.1
## 45 74 7 2 1113603 0.1
## 46 75 7 2 1251009 0.1
## 47 76 7 2 1138017 0.1
## 48 77 7 2 935018 0.1
## 49 78 7 2 903728 0.1
## 50 79 7 2 1094709 0.1stocksynthesis.data$catch## year seas fleet catch catch_se
## 1 30 1 1 0 0.01
## 2 31 1 1 319865 0.01
## 3 32 1 1 249922 0.01
## 4 33 1 1 215679 0.01
## 5 34 1 1 220781 0.01
## 6 35 1 1 216810 0.01
## 7 36 1 1 203003 0.01
## 8 37 1 1 193817 0.01
## 9 38 1 1 185898 0.01
## 10 39 1 1 180176 0.01
## 11 40 1 1 176713 0.01
## 12 41 1 1 175377 0.01
## 13 42 1 1 182145 0.01
## 14 43 1 1 64694 0.01
## 15 44 1 1 78720 0.01
## 16 45 1 1 91349 0.01
## 17 46 1 1 103954 0.01
## 18 47 1 1 116472 0.01
## 19 48 1 1 127937 0.01
## 20 49 1 1 137312 0.01
## 21 50 1 1 140571 0.01
## 22 51 1 1 138208 0.01
## 23 52 1 1 136455 0.01
## 24 53 1 1 135849 0.01
## 25 54 1 1 140378 0.01
## 26 55 1 1 1018985 0.01
## 27 56 1 1 424892 0.01
## 28 57 1 1 157952 0.01
## 29 58 1 1 101970 0.01
## 30 59 1 1 73612 0.01
## 31 60 1 1 141543 0.01
## 32 61 1 1 185664 0.01
## 33 62 1 1 315029 0.01
## 34 63 1 1 324810 0.01
## 35 64 1 1 278613 0.01
## 36 65 1 1 245273 0.01
## 37 66 1 1 174959 0.01
## 38 67 1 1 178687 0.01
## 39 68 1 1 156851 0.01
## 40 69 1 1 249209 0.01
## 41 70 1 1 328456 0.01
## 42 71 1 1 366142 0.01
## 43 72 1 1 264745 0.01
## 44 73 1 1 240040 0.01
## 45 74 1 1 270350 0.01
## 46 75 1 1 268981 0.01
## 47 76 1 1 245319 0.01
## 48 77 1 1 309845 0.01
## 49 78 1 1 202432 0.01
## 50 79 1 1 254846 0.01Writing age composition data
#age_comp_data
#Get the age bins
age_bin_names <- names(stocksynthesis.data$agecomp)[10:length(names(stocksynthesis.data$agecomp))]
age_bins <- sub("a","",age_bin_names)
age_comp_flat <- reformat_compositions(age_comp_data,
round.places = 4,
comp_type = "agecomp")
## Write age composition data for survey
stocksynthesis.data <- SS_write_comps(ss_data_list = stocksynthesis.data,
comp_matrix = list(age_comp_flat[burnin:(burnin+nyears-1),]),
data_rows = list(stocksynthesis.data$styr:(stocksynthesis.data$styr+nyears-1)),
sampling_month = list(rep(survey_month,nyears)),
data_type = c("agecomp"),
fleet_number = c(1),
bins = list(age_bins),
caal_bool = c(FALSE))
stocksynthesis.data$agecomp## Yr Seas FltSvy Gender Part Ageerr Lbin_lo Lbin_hi Nsamp a1 a2
## 1 1 7 1 0 0 1 -1 -1 1 0.106 0.015
## 2 2 7 1 0 0 1 -1 -1 1 0.079 0.080
## 3 3 7 1 0 0 1 -1 -1 1 0.062 0.105
## 4 4 7 1 0 0 1 -1 -1 1 0.231 0.055
## 5 5 7 1 0 0 1 -1 -1 1 0.025 0.322
## 6 6 7 1 0 0 1 -1 -1 1 0.117 0.028
## 7 7 7 1 0 0 1 -1 -1 1 0.351 0.119
## 8 8 7 1 0 0 1 -1 -1 1 0.052 0.391
## 9 9 7 1 0 0 1 -1 -1 1 0.019 0.069
## 10 10 7 1 0 0 1 -1 -1 1 0.536 0.006
## 11 11 7 1 0 0 1 -1 -1 1 0.726 0.165
## 12 12 7 1 0 0 1 -1 -1 1 0.166 0.649
## 13 13 7 1 0 0 1 -1 -1 1 0.121 0.163
## 14 14 7 1 0 0 1 -1 -1 1 0.748 0.047
## 15 15 7 1 0 0 1 -1 -1 1 0.060 0.740
## 16 16 7 1 0 0 1 -1 -1 1 0.012 0.040
## 17 17 7 1 0 0 1 -1 -1 1 0.460 0.005
## 18 18 7 1 0 0 1 -1 -1 1 0.219 0.380
## 19 19 7 1 0 0 1 -1 -1 1 0.165 0.222
## 20 20 7 1 0 0 1 -1 -1 1 0.110 0.182
## 21 21 7 1 0 0 1 -1 -1 1 0.362 0.082
## 22 22 7 1 0 0 1 -1 -1 1 0.192 0.335
## 23 23 7 1 0 0 1 -1 -1 1 0.175 0.215
## 24 24 7 1 0 0 1 -1 -1 1 0.199 0.163
## 25 25 7 1 0 0 1 -1 -1 1 0.349 0.165
## 26 26 7 1 0 0 1 -1 -1 1 0.746 0.079
## 27 27 7 1 0 0 1 -1 -1 1 0.131 0.691
## 28 28 7 1 0 0 1 -1 -1 1 0.248 0.100
## 29 29 7 1 0 0 1 -1 -1 1 0.100 0.232
## 30 30 7 1 0 0 1 -1 -1 1 0.080 0.091
## 31 31 7 1 0 0 1 -1 -1 1 0.277 0.067
## 32 32 7 1 0 0 1 -1 -1 1 0.065 0.281
## 33 33 7 1 0 0 1 -1 -1 1 0.161 0.055
## 34 34 7 1 0 0 1 -1 -1 1 0.125 0.152
## 35 35 7 1 0 0 1 -1 -1 1 0.173 0.116
## 36 36 7 1 0 0 1 -1 -1 1 0.048 0.202
## 37 37 7 1 0 0 1 -1 -1 1 0.063 0.054
## 38 38 7 1 0 0 1 -1 -1 1 0.529 0.029
## 39 39 7 1 0 0 1 -1 -1 1 0.092 0.480
## 40 40 7 1 0 0 1 -1 -1 1 0.015 0.138
## 41 41 7 1 0 0 1 -1 -1 1 0.052 0.022
## 42 42 7 1 0 0 1 -1 -1 1 0.032 0.036
## 43 43 7 1 0 0 1 -1 -1 1 0.101 0.019
## 44 44 7 1 0 0 1 -1 -1 1 0.528 0.050
## 45 45 7 1 0 0 1 -1 -1 1 0.034 0.590
## 46 46 7 1 0 0 1 -1 -1 1 0.013 0.029
## 47 47 7 1 0 0 1 -1 -1 1 0.033 0.017
## 48 48 7 1 0 0 1 -1 -1 1 0.260 0.048
## 49 49 7 1 0 0 1 -1 -1 1 0.108 0.283
## 50 50 7 1 0 0 1 -1 -1 1 0.384 0.078
## a3 a4 a5 a6 a7 a8 a9 a10
## 1 0.068 0.053 0.251 0.094 0.283 0.102 0.014 0.014
## 2 0.013 0.049 0.071 0.238 0.094 0.256 0.102 0.018
## 3 0.086 0.008 0.051 0.034 0.224 0.095 0.243 0.092
## 4 0.088 0.064 0.008 0.048 0.056 0.195 0.064 0.191
## 5 0.067 0.090 0.065 0.005 0.045 0.045 0.246 0.090
## 6 0.316 0.071 0.085 0.073 0.007 0.054 0.039 0.210
## 7 0.023 0.237 0.044 0.077 0.057 0.012 0.044 0.036
## 8 0.095 0.015 0.239 0.052 0.058 0.052 0.005 0.041
## 9 0.391 0.096 0.023 0.231 0.053 0.062 0.047 0.009
## 10 0.033 0.181 0.055 0.009 0.102 0.022 0.029 0.027
## 11 0.004 0.011 0.046 0.009 0.002 0.024 0.007 0.006
## 12 0.117 0.001 0.004 0.033 0.007 0.000 0.023 0.000
## 13 0.559 0.094 0.002 0.012 0.032 0.008 0.000 0.009
## 14 0.027 0.147 0.023 0.000 0.001 0.005 0.002 0.000
## 15 0.030 0.033 0.114 0.018 0.001 0.000 0.002 0.002
## 16 0.735 0.040 0.038 0.112 0.016 0.000 0.000 0.007
## 17 0.023 0.381 0.022 0.019 0.082 0.008 0.000 0.000
## 18 0.005 0.022 0.292 0.016 0.010 0.052 0.004 0.000
## 19 0.297 0.002 0.018 0.237 0.010 0.004 0.034 0.011
## 20 0.191 0.260 0.006 0.014 0.203 0.007 0.009 0.018
## 21 0.112 0.111 0.171 0.002 0.010 0.134 0.009 0.007
## 22 0.067 0.083 0.096 0.114 0.002 0.009 0.099 0.003
## 23 0.256 0.037 0.067 0.064 0.110 0.002 0.003 0.071
## 24 0.165 0.218 0.038 0.067 0.056 0.089 0.001 0.004
## 25 0.094 0.099 0.136 0.026 0.031 0.038 0.060 0.002
## 26 0.039 0.032 0.026 0.038 0.010 0.011 0.006 0.013
## 27 0.074 0.027 0.008 0.017 0.029 0.004 0.010 0.009
## 28 0.518 0.057 0.014 0.012 0.017 0.023 0.003 0.008
## 29 0.094 0.464 0.056 0.017 0.009 0.011 0.012 0.005
## 30 0.222 0.072 0.448 0.049 0.016 0.007 0.006 0.009
## 31 0.059 0.170 0.051 0.300 0.041 0.015 0.013 0.007
## 32 0.055 0.069 0.156 0.049 0.282 0.031 0.008 0.004
## 33 0.245 0.041 0.064 0.146 0.040 0.219 0.024 0.005
## 34 0.056 0.213 0.045 0.049 0.107 0.040 0.191 0.022
## 35 0.123 0.054 0.175 0.029 0.047 0.090 0.041 0.152
## 36 0.128 0.131 0.060 0.211 0.033 0.046 0.107 0.034
## 37 0.202 0.115 0.140 0.061 0.182 0.032 0.047 0.104
## 38 0.027 0.107 0.063 0.088 0.038 0.084 0.016 0.019
## 39 0.037 0.037 0.102 0.051 0.061 0.032 0.095 0.013
## 40 0.469 0.038 0.024 0.091 0.059 0.065 0.024 0.077
## 41 0.132 0.488 0.042 0.029 0.108 0.052 0.052 0.023
## 42 0.028 0.124 0.491 0.033 0.029 0.102 0.061 0.064
## 43 0.049 0.030 0.116 0.504 0.030 0.034 0.080 0.037
## 44 0.010 0.019 0.014 0.072 0.235 0.015 0.012 0.045
## 45 0.051 0.007 0.025 0.014 0.053 0.204 0.012 0.010
## 46 0.560 0.053 0.009 0.015 0.013 0.082 0.215 0.011
## 47 0.039 0.554 0.057 0.011 0.020 0.009 0.058 0.202
## 48 0.012 0.023 0.525 0.048 0.005 0.020 0.009 0.050
## 49 0.026 0.013 0.023 0.475 0.037 0.011 0.018 0.006
## 50 0.163 0.022 0.005 0.020 0.293 0.027 0.002 0.006Once we have the natlength matrix, we can munge the data into the proper CAAL and length bin format for SS
#Check length compositions match age compositions
#ss_length_stdsurv$natlength %>%
# group_by(time,agecl) %>%
# summarise(sum(atoutput))
len_comp_data <- ss_length_stdsurv$natlength
fish_len_comp_data <- catch_lengthwt_samp$natlength
if(fstepperyr>1){
fish_len_comp_anndata <- fish_len_comp_data %>%
mutate(yr = floor(time/fstepperyr)) %>%
group_by(species, agecl, lower.bins, upper.bins, time=as.integer(yr)) %>%
summarise(annnatlength=sum(atoutput)) %>%
rename(atoutput = annnatlength)
} else {
fish_len_comp_anndata <- fish_len_comp_data
}
caal_comp_flat <- reformat_compositions(len_comp_data, round.places=4,
comp_type="caalcomp")
#remove burnin
caal_comp_final <- filter(caal_comp_flat,
time %in% survey_years)
#Add over age classes to get sample size
len_comp_flat <- reformat_compositions(len_comp_data,
round.places = 0,
comp_type="lencomp")
#remove burnin
len_comp_final <- filter(len_comp_flat,
time %in% survey_years)
length_bins <- as.integer(names(len_comp_final))## Warning: NAs introduced by coercionlength_bins <- length_bins[!is.na(length_bins)]
# fishery length comps are still 5 timesteps per year
# need to aggregate to annual (done above)
# also, make effN annual goal/fstepsperyr (done above)
fish_len_comp_flat <- reformat_compositions(fish_len_comp_anndata,
round.places = 0,
comp_type="lencomp")
#remove burnin works after adjustment above
fish_len_comp_final <- filter(fish_len_comp_flat,
time %in% fish_years)
notbins <- c("time", "nsamp")
# fish_length_bins <- as.integer(names(fish_len_comp_final))
# fish_length_bins <- fish_length_bins[!is.na(fish_length_bins)]
# need to fill empty length bins with 0s to have same bins as survey for SS_write_comps
missing.lengths <- setdiff(length_bins, names(fish_len_comp_final)[!names(fish_len_comp_final) %in% notbins])
fish_len_comp_final[as.character(missing.lengths)] <- 0 # Add them, filled with '0's
fish_len_comp_final <- fish_len_comp_final[c("time", length_bins, "nsamp")]
# fishery age comps also 5 timesteps per year
if(fstepperyr>1){
fish_age_comp_anndata <- catch_numsss_samp %>%
mutate(yr = floor(time/fstepperyr)) %>%
group_by(species, agecl, time=as.integer(yr)) %>%
summarise(annnatage=sum(atoutput)) %>%
rename(atoutput = annnatage)
} else {
fish_age_comp_anndata <- catch_numsss_samp
}
fish_age_comp_flat <- reformat_compositions(fish_age_comp_anndata,
comp_type="agecomp")
#remove burnin (not necessary?fish comps made with fish_years only)
fish_age_comp_final <- filter(fish_age_comp_flat,
time %in% fish_years)
# #SS_write_comps breaking because fishery age bins start with 2 not 1; extracting bins from fish file may help?
# fish_age_bins <- names(fish_age_comp_flat)[!names(fish_age_comp_flat) %in% notbins]
# that leaves an empty column in data file, so instead fill with 0s
missing.ages <- setdiff(age_bins, names(fish_age_comp_final)[!names(fish_age_comp_final) %in% notbins])
fish_age_comp_final[missing.ages] <- 0 # Add them, filled with '0's
fish_age_comp_final <- fish_age_comp_final[c("time", age_bins, "nsamp")]
comp_list <- list(caal_comp_final,len_comp_final, fish_age_comp_final, fish_len_comp_final)
apply_month <- list(rep(survey_month, nrow(comp_list[[1]])),
rep(survey_month, nrow(comp_list[[2]])),
rep(fishing_month,nrow(comp_list[[3]])),
rep(fishing_month,nrow(comp_list[[4]])))
# This now runs by ensuring that survey and fishery compositions have the same bins
# (filled with 0s for missing bins in fishery relative to survey)
# Write CAAL and length composition data
stocksynthesis.data <- SS_write_comps(ss_data_list = stocksynthesis.data,
comp_matrix = comp_list,
data_rows = list((comp_list[[1]]$time-survey_sample_time)/timestep + 1 , (survey_years-survey_sample_time)/timestep + 1,fish_years,fish_years),
sampling_month = apply_month,
data_type = rep(c("agecomp", "lencomp"),2),
fleet_number = c(2,2,1,1),
bins = list(age_bins,
length_bins,
age_bins,
length_bins),
caal_bool = c(TRUE, rep(FALSE,3)))
head(stocksynthesis.data$lencomp)## Yr Seas FltSvy Gender Part Nsamp l12 l13 l14 l15 l16 l17 l18 l19 l20 l21
## 1 30 7 2 0 0 988 1 5 14 24 25 19 12 7 5 5
## 2 31 7 2 0 0 995 1 4 10 17 19 17 12 13 14 15
## 3 32 7 2 0 0 991 0 3 7 13 15 14 12 13 17 21
## 4 33 7 2 0 0 993 2 9 26 47 55 46 29 18 13 13
## 5 34 7 2 0 0 988 0 1 3 5 7 8 14 26 42 55
## 6 35 7 2 0 0 992 1 4 12 23 28 24 16 9 8 10
## l22 l23 l24 l25 l26 l27 l28 l29 l30 l31 l32 l33 l34 l35 l36 l37 l38 l39
## 1 6 10 11 14 17 21 27 32 39 48 54 60 63 65 64 62 57 50
## 2 15 11 10 8 8 11 12 17 23 30 36 45 53 58 65 66 66 63
## 3 22 22 20 18 16 15 13 13 12 17 19 25 31 41 47 55 60 61
## 4 13 14 16 17 18 18 17 17 15 16 16 19 22 26 31 38 42 46
## 5 59 53 44 34 27 24 22 21 20 21 21 22 20 23 26 28 33 35
## 6 15 23 34 44 50 52 48 44 37 31 29 27 25 25 27 28 29 30
## l40 l41 l42 l43 l44 l45 l46 l47 l48 l49 l50 l51 l52 l53 l54 l55 l56 l57
## 1 42 35 27 22 15 11 8 5 3 2 1 0 0 0 0 0 0 0
## 2 57 50 42 36 27 20 15 11 7 5 2 2 2 0 0 0 0 0
## 3 61 58 54 47 38 31 24 18 13 10 6 4 2 2 1 0 0 0
## 4 48 47 44 42 37 29 24 20 14 10 7 5 3 2 1 1 0 0
## 5 38 39 39 35 33 27 23 18 13 10 7 5 3 2 2 0 0 0
## 6 31 32 32 29 27 24 20 17 13 11 8 6 3 2 2 1 1 0
## l58 l59
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 0 0head(stocksynthesis.data$agecomp)## Yr Seas FltSvy Gender Part Ageerr Lbin_lo Lbin_hi Nsamp a1 a2 a3 a4
## 1 30 7 2 0 0 1 12 13 0.0010 0.0010 0 0 0
## 2 30 7 2 0 0 1 13 14 0.0051 0.0051 0 0 0
## 3 30 7 2 0 0 1 14 15 0.0142 0.0142 0 0 0
## 4 30 7 2 0 0 1 15 16 0.0243 0.0243 0 0 0
## 5 30 7 2 0 0 1 16 17 0.0253 0.0253 0 0 0
## 6 30 7 2 0 0 1 17 18 0.0192 0.0192 0 0 0
## a5 a6 a7 a8 a9 a10
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0#Change length bin structure to match atlantis data
stocksynthesis.data$lbin_vector <- length_bins
#Get correct number of length bins
stocksynthesis.data$N_lbins <- length(length_bins)
#Set lbin_method to 1 - this makes the population length bins match the data bins
#When lbin_method==1, we just comment out the binwidth, minimum, and maximum size arguments since they aren't used
stocksynthesis.data$lbin_method <- 1
stocksynthesis.data$binwidth <- "#"
stocksynthesis.data$minimum_size <- "#"
stocksynthesis.data$maximum_size <- "#"
#Change minimum sample size to 0.001 for CAAL data (SS won't let it go lower than this)
stocksynthesis.data$age_info$minsamplesize <- rep(0.001,2)
SS_writedat_3.30(stocksynthesis.data, outfile = paste0("./inst/extdata/",model_dir,
datfile_name),
overwrite=TRUE)## running SS_writedat_3.30## opening connection to ./inst/extdata/Sardine_SS_files/sardEM_3_3.dat## file written to ./inst/extdata/Sardine_SS_files/sardEM_3_3.datComing up with dynamic life history parameters (partial from Christine)
Below we get the weight-length relationship parameters from atlantis, calculate h and R0 from Atlantis \(\alpha\) and \(\beta\) parameters, use the function to back out natural mortality (M) from Atlantis, and use the survey length-at-age from Atlantis to estimate a growth curve.
The following is under construction from CreateStockSynthesis.Rmd and not evaluated. Existing control file was modified by hand to make the number of recruitment settlement assignments 1 instead of two, with GPattern month area age vector 1 1 1 1.
#Get biological parameters
#Load needed inputs for biological parameters
source(here("config/CC3Config.R"))
#biological parameters
biolprm <- load_biolprm(dir = d.name, file_biolprm = biol.prm.file)#"C:/Users/chris/Documents/GitHub/atlantisom",
#d.name <- "C:/Users/chris/Documents/GitHub/atlantisom"
#functional groups
fgs <- truth$fgs
runprm <- load_runprm(dir = d.name, file_runprm = run.prm.file)
YOY <- load_yoy(d.name, paste0("output", scenario.name, "YOY.txt"))
truenums_ss <- truth$nums[results$nums$species==species,]
YOY_ss <- YOY %>%
select(Time, "SAR.0")
fullresZ <- calc_Z(yoy = YOY_ss,
nums = truenums_ss,
fgs = fgs,
biolprm = biolprm,
toutinc = runprm$toutinc)
meanwt_spp <- ss_length_stdsurv$muweight %>%
filter(time>burnin) %>%
group_by(agecl) %>%
summarize(meanwt = mean(atoutput))
ctlfile_name <- "sardEM_3_3.ctl"
sardine.ctl <-r4ss::SS_readctl_3.30(paste0("./inst/extdata/",
model_dir,
ctlfile_name))
sardine.ctl <- SS_write_biol(sardine.ctl, biolprm, "SAR", Z = mean(fullresZ$atoutput), wtsage=meanwt_spp)
li_a_use <- biolprm$wl[match(fgs$Code[match(species,fgs$Name)],biolprm$wl[, 1]), 2]/1000
li_b_use <- biolprm$wl[match(fgs$Code[match(species,fgs$Name)],biolprm$wl[, 1]), 3]
#maturity ogive
#check what R0 is
len_nonburn <- ss_length_stdsurv$mulen %>%
filter(time > burnin*timestep)
plot(len_nonburn$atoutput~len_nonburn$agecl, col=len_nonburn$time)
length.data <- data.frame("Year"=(len_nonburn$time-survey_sample_time)/timestep+1, length=len_nonburn$atoutput, Weight=NA, Sex="Female", age=as.integer(len_nonburn$agecl))
#require(bbmle)
#vb_estimates <-sample_fit_vbgf(length.data, 25, 45, 0.4,
# 0.1, 0.1, 20, 40, 0.05, 0.01, 0.01,
# 30, 50, 0.6, 0.3, 0.3, 0.5, 10)Run Stock Synthesis
After the data and control file modifications, ensure you copy the local SS 3.3 executable into the folder containing the model files. Then run the following chunk.
#Run stock synthesis model
# default runs with admb_options = "-noest"
#run_stocksynthesis(species = species, model_dir = model_dir, show_output = TRUE)
# for estimation
run_stocksynthesis(species = species, model_dir = model_dir, admb_options = "", show_output = TRUE)## Running ss EM for species: Pacific_sardine