Generate data for NOBA cod ensemble modeling
Sarah Gaichas
05 June, 2020
Standard atlantisom output test (age class not age)
Let’s use the wrapper functions to get initial stock assessment model inputs.
Then we need to adjust the base atlantisom functions for true age classes, and develop options to call those within the wrappers.
We need config files for the NOBA model, our surveys and fisheries. Surveys and fisheries are based on the ones we used for sardine in the California Current as a test.
Special note: Atlantis output files formerly had the prefix “output” and atlantisom was coded to expect this. NOBA output files (and newer Atlantis models) do not universally use the “output” prefix naming convention. Annage branch atlantisom functions have been updated for user-specified output filenames.
Special note 2: NOBA’s nordic_groups_v04.csv file has an inFishery column, rather than the atlantisom expected InFishery column. I’ve changed the column name in the input file, but need to add robustness to atlantisom code to remove case sensitivity.
Config files
NOBA2config.R looks like this (adjusted from Alfonso’s original):
d.name <- here("atlantisoutput","NOBA_march_2020")
functional.groups.file <- "nordic_groups_v04.csv"
biomass.pools.file <- "nordic_biol_v23.nc"
biol.prm.file <- "nordic_biol_incl_harv_v_007_3.prm"
box.file <- "Nordic02.bgm"
initial.conditions.file <- "nordic_biol_v23.nc"
run.prm.file <- "nordic_run_v01.xml"
scenario.name <- "nordic_runresults_01"
bioind.file <- "nordic_runresults_01BiomIndx.txt"
catch.file <- "nordic_runresults_01Catch.txt"
annage <- TRUE
fisheries.file <- "NoBAFisheries.csv"omdimensions.R standardizes timesteps, etc.:
#survey species inherited from omlist_ss
survspp <- omlist_ss$species_ss
# survey season and other time dimensioning parameters
# generalized timesteps all models
noutsteps <- omlist_ss$runpar$tstop/omlist_ss$runpar$outputstep
timeall <- c(0:noutsteps)
stepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutinc
midptyr <- round(median(seq(0,stepperyr)))
# model areas, subset in surveyconfig
allboxes <- c(0:(omlist_ss$boxpars$nbox - 1))
# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc
# survey selectivity (agecl based)
# should return all age classes fully sampled (Atlantis output is 10 age groups per spp)
n_age_classes <- omlist_ss$funct.group_ss$NumCohorts
# changed below for multiple species
age_classes <- sapply(n_age_classes, seq)
names(age_classes)<-survspp
n_annages <- n_age_classes * omlist_ss$funct.group_ss$NumAgeClassSize
# changed below for multiple species
annages <- sapply(n_annages, seq)
names(annages)<-survsppcodsurvey.R configures the fishery independent survey and looks like this:
# Default survey configuration here is a census
# Need to define survey season, area, efficiency, selectivity
#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 #2; defined in omdimensions.R
#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)
survtime <- survey_sample_full
# survey area
# should return all model areas
survboxes <- allboxes
# survey efficiency (q)
# should return a perfectly efficient survey
surveffic <- data.frame(species=survspp,
efficiency=rep(1.0,length(survspp)))
# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(survspp, each=n_age_classes),
# agecl=rep(c(1:n_age_classes),length(survspp)),
# selex=rep(1.0,length(survspp)*n_age_classes))
# for annage output
survselex <- data.frame(species=rep(survspp, n_annages), #
agecl=unlist(sapply(n_annages,seq)),
selex=rep(1.0,sum(n_annages)))
# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))
# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))
# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1
# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200codfishery.R configures the fishery dependent data and looks like this:
# Default fishery configuration here is a census
# Inherits species from input omlist_ss
fishspp <- omlist_ss$species_ss
#Number of years of data to pull
nyears <- 50
#Atlantis initialization period in years
burnin <- 30
# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc
# same time dimensioning parameters as in surveycensus.R
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample) #total_sample defined in sardinesurvey.R
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[fstepperyr], 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
fishtime <- fish_times
# fishery sampling area
# should return all model areas, this assumes you see everything that it caught
fishboxes <- c(0:(omlist_ss$boxpars$nbox - 1))
# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
fisheffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))
#this adjusts for subannual fishery output so original effN is for whole year
fisheffN$effN <- fisheffN$effN/fstepperyr
# fishery catch cv can be used in sample_survey_biomass
# perfect observation
fish_cv <- data.frame(species=survspp, cv=rep(0.01,length(survspp)))Using the wrappers (now updated for annual age)
All the config files go in the config folder and we’ll try running this:
#library(here) now above
library(tidyverse)
library(atlantisom)
library(ggthemes)
NOBAom <- om_init(here("config/NOBA2config.R"))
NOBAom_cod <- om_species(c("North_atl_cod"), NOBAom)
NOBAom_cod_ind <- om_index(usersurvey = here("config/codsurvey.R"),
userfishery = here("config/codfishery.R"),
omlist_ss = NOBAom_cod,
n_reps = 1,
save = TRUE)
NOBAom_cod_comp <- om_comps(usersurvey = here("config/codsurvey.R"),
userfishery = here("config/codfishery.R"),
omlist_ss = NOBAom_cod,
n_reps = 1,
save = TRUE)Wrapper test results
Biomass index
omlist_ss <- NOBAom_cod
source(here("config/omdimensions.R"))
#read time series data
survObsBiom <- readRDS(file.path(d.name, paste0(scenario.name, "surveyB.rds")))
survObsBiom <- survObsBiom[[1]]
plotB <-ggplot() +
geom_line(data=survObsBiom, aes(x=time/stepperyr,y=atoutput, color="survey Biomass"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
plotB +
facet_wrap(~species, scales="free") 
Catch time series
#read time series data
catchbio_ss <- readRDS(file.path(d.name, paste0(scenario.name, "fishCatch.rds")))
catchbio_ss <- catchbio_ss[[1]]
plotC <-ggplot() +
geom_line(data=catchbio_ss, aes(x=time/365,y=atoutput, color="observed Catch"),
alpha = 10/10) +
theme_tufte() +
theme(legend.position = "top") +
labs(colour=scenario.name)
plotC +
facet_wrap(~species, scales="free") 
Survey length composition
#length comps
len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsLenComp.rds")))
fish_len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsLenComp.rds")))
len_comp_data <- len_comp_data[[1]]
fish_len_comp_data <- fish_len_comp_data[[1]]
#add this to om_indices function so that this has years when read in
fish_len_comp_data$time <- as.integer(floor(fish_len_comp_data$time/fstepperyr))
len <- filter(len_comp_data, time %in% c(55:175))
lfplot <- ggplot(len, aes(upper.bins)) +
geom_bar(aes(weight = atoutput)) +
theme_tufte() +
labs(subtitle = paste(scenario.name,
len$species))
lfplot + facet_wrap(~time/stepperyr, ncol=6, scales="free_y")
Survey age composition (age classes)
#read in comp data
age_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsAgeComp.rds")))
age_comp_data <- age_comp_data[[1]]
Natage <- filter(age_comp_data, time %in% c(150:270))
Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste(scenario.name,
Natage$species))
Natageplot + facet_wrap(~time/stepperyr, ncol=6, scales="free_y")
Survey weight at age (age classes)
wtage <- readRDS(file.path(d.name, paste0(scenario.name, "survObsWtAtAge.rds")))
wtage <- wtage[[1]]
wageplot <- ggplot(wtage, 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")
Fishery age composition (age classes)
#read in comp data
fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]
#add this to om_indices function so that this has years when read in
fish_age_comp$time <- fish_age_comp$time/fstepperyr
Natage <- filter(fish_age_comp, time %in% c(30:53))
Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste(scenario.name,
Natage$species))
Natageplot + facet_wrap(~time, ncol=6, scales="free_y")
Wrappers still work for NOBA cod using standard age classes.
Visualize new annual age output for numbers at age, catch at age, and weight at age (interpolated):
Wrapper test results for annual ages
Survey age composition (annual ages)
#read in comp data
annage_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsFullAgeComp.rds")))
annage_comp_data <- annage_comp_data[[1]]
#fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
#fish_age_comp <- fish_age_comp[[1]]
#add this to om_indices function so that this has years when read in
#fish_age_comp$time <- as.integer(floor(fish_age_comp$time/fstepperyr))
Natage <- filter(annage_comp_data, time %in% c(150:270))
Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste(scenario.name,
Natage$species))
Natageplot + facet_wrap(~time/stepperyr, ncol=6, scales="free_y")
Survey iterpolated weight at age (annual ages)
wtage <- readRDS(file.path(d.name, paste0(scenario.name, "survObsFullWtAtAge.rds")))
wtage <- wtage[[1]] #this still has second list component, diagnostic plot
wtage <- wtage[[1]]
wageplot <- ggplot(wtage, 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")
Fishery catch at age (annual ages)
fish_annage_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullAgeComp.rds")))
fish_annage_comp <- fish_annage_comp[[1]]
#add this to om_indices function so that this has years when read in
fish_annage_comp$time <- fish_annage_comp$time/fstepperyr
Catage <- filter(fish_annage_comp, time %in% c(30:53))
Catageplot <- ggplot(Catage, aes(x=agecl, y=atoutput)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste(scenario.name,
Catage$species))
Catageplot + facet_wrap(~time, ncol=6, scales="free_y")
Fishery iterpolated weight at age (annual ages)
wtage <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullWtAtAge.rds")))
wtage <- wtage[[1]] #this still has second list component, diagnostic plot
wtage <- wtage[[1]]
wageplot <- ggplot(wtage, 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")
Note: NOBA doesn’t seem to have a selective fishery for cod. Lots of age 1 in the catch!
(Annage development details below)
Now to modify the underlying code to get full age structure.
Updating functions for full age structure
General strategy:
- add annual age option to run_truth function call, default it to FALSE
- if true, looks for ANNAGEBIO and ANNAGECATCH in designated directory, error out if missing
- if true, read these files in addition to existing set rather than replace existing
- read them in with new load_nc_annage function
- add them to the output object in addition to existing set
- write load_nc_annage
- need to understand annage nc structures
- if same as others just with full age structure, replace previous nc
- only numbers output?
- how to match with resn and structn for lengths/wts?
- reconcile with bio parameters such as maturity ogive?
- only for some groups?
What is in the ANNAGEBIO and ANNAGECATCH .nc files?
# taken from existing load_nc
file.nc <- file.path(d.name, "outputnordic_runresults_01ANNAGEBIO.nc")
# Load ATLANTIS output!
at_out <- RNetCDF::open.nc(con = file.nc)
# 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
RNetCDF::close.nc(at_out)ANNAGEBIO.nc variable names (first hundred): volume, nominal_dz, dz, topk, numlayers, Polar_bear1_Nums, Polar_bear1_Weight, Polar_bear2_Nums, Polar_bear2_Weight, Polar_bear3_Nums, Polar_bear3_Weight, Polar_bear4_Nums, Polar_bear4_Weight, Polar_bear5_Nums, Polar_bear5_Weight, Polar_bear6_Nums, Polar_bear6_Weight, Polar_bear7_Nums, Polar_bear7_Weight, Polar_bear8_Nums, Polar_bear8_Weight, Polar_bear9_Nums, Polar_bear9_Weight, Polar_bear10_Nums, Polar_bear10_Weight, Polar_bear11_Nums, Polar_bear11_Weight, Polar_bear12_Nums, Polar_bear12_Weight, Polar_bear13_Nums, Polar_bear13_Weight, Polar_bear14_Nums, Polar_bear14_Weight, Polar_bear15_Nums, Polar_bear15_Weight, Polar_bear16_Nums, Polar_bear16_Weight, Polar_bear17_Nums, Polar_bear17_Weight, Polar_bear18_Nums, Polar_bear18_Weight, Polar_bear19_Nums, Polar_bear19_Weight, Polar_bear20_Nums, Polar_bear20_Weight, Killer_whale1_Nums, Killer_whale1_Weight, Killer_whale2_Nums, Killer_whale2_Weight, Killer_whale3_Nums, Killer_whale3_Weight, Killer_whale4_Nums, Killer_whale4_Weight, Killer_whale5_Nums, Killer_whale5_Weight, Killer_whale6_Nums, Killer_whale6_Weight, Killer_whale7_Nums, Killer_whale7_Weight, Killer_whale8_Nums, Killer_whale8_Weight, Killer_whale9_Nums, Killer_whale9_Weight, Killer_whale10_Nums, Killer_whale10_Weight, Killer_whale11_Nums, Killer_whale11_Weight, Killer_whale12_Nums, Killer_whale12_Weight, Killer_whale13_Nums, Killer_whale13_Weight, Killer_whale14_Nums, Killer_whale14_Weight, Killer_whale15_Nums, Killer_whale15_Weight, Killer_whale16_Nums, Killer_whale16_Weight, Killer_whale17_Nums, Killer_whale17_Weight, Killer_whale18_Nums, Killer_whale18_Weight, Killer_whale19_Nums, Killer_whale19_Weight, Killer_whale20_Nums, Killer_whale20_Weight, Killer_whale21_Nums, Killer_whale21_Weight, Killer_whale22_Nums, Killer_whale22_Weight, Killer_whale23_Nums, Killer_whale23_Weight, Killer_whale24_Nums, Killer_whale24_Weight, Killer_whale25_Nums, Killer_whale25_Weight, Killer_whale26_Nums, Killer_whale26_Weight, Killer_whale27_Nums, Killer_whale27_Weight, Killer_whale28_Nums
ANNABEBIO.nc variable names (last hundred): Norwegian_ssh13_Nums, Norwegian_ssh13_Weight, Norwegian_ssh14_Nums, Norwegian_ssh14_Weight, Norwegian_ssh15_Nums, Norwegian_ssh15_Weight, Norwegian_ssh16_Nums, Norwegian_ssh16_Weight, Norwegian_ssh17_Nums, Norwegian_ssh17_Weight, Norwegian_ssh18_Nums, Norwegian_ssh18_Weight, Norwegian_ssh19_Nums, Norwegian_ssh19_Weight, Norwegian_ssh20_Nums, Norwegian_ssh20_Weight, North_atl_cod1_Nums, North_atl_cod1_Weight, North_atl_cod2_Nums, North_atl_cod2_Weight, North_atl_cod3_Nums, North_atl_cod3_Weight, North_atl_cod4_Nums, North_atl_cod4_Weight, North_atl_cod5_Nums, North_atl_cod5_Weight, North_atl_cod6_Nums, North_atl_cod6_Weight, North_atl_cod7_Nums, North_atl_cod7_Weight, North_atl_cod8_Nums, North_atl_cod8_Weight, North_atl_cod9_Nums, North_atl_cod9_Weight, North_atl_cod10_Nums, North_atl_cod10_Weight, North_atl_cod11_Nums, North_atl_cod11_Weight, North_atl_cod12_Nums, North_atl_cod12_Weight, North_atl_cod13_Nums, North_atl_cod13_Weight, North_atl_cod14_Nums, North_atl_cod14_Weight, North_atl_cod15_Nums, North_atl_cod15_Weight, North_atl_cod16_Nums, North_atl_cod16_Weight, North_atl_cod17_Nums, North_atl_cod17_Weight, North_atl_cod18_Nums, North_atl_cod18_Weight, North_atl_cod19_Nums, North_atl_cod19_Weight, North_atl_cod20_Nums, North_atl_cod20_Weight, Polar_cod1_Nums, Polar_cod1_Weight, Polar_cod2_Nums, Polar_cod2_Weight, Polar_cod3_Nums, Polar_cod3_Weight, Polar_cod4_Nums, Polar_cod4_Weight, Polar_cod5_Nums, Polar_cod5_Weight, Polar_cod6_Nums, Polar_cod6_Weight, Polar_cod7_Nums, Polar_cod7_Weight, Polar_cod8_Nums, Polar_cod8_Weight, Polar_cod9_Nums, Polar_cod9_Weight, Polar_cod10_Nums, Polar_cod10_Weight, Capelin1_Nums, Capelin1_Weight, Capelin2_Nums, Capelin2_Weight, Capelin3_Nums, Capelin3_Weight, Capelin4_Nums, Capelin4_Weight, Capelin5_Nums, Capelin5_Weight, Snow_crab1_Nums, Snow_crab1_Weight, Snow_crab2_Nums, Snow_crab2_Weight, Snow_crab3_Nums, Snow_crab3_Weight, Snow_crab4_Nums, Snow_crab4_Weight, Snow_crab5_Nums, Snow_crab5_Weight, Snow_crab6_Nums, Snow_crab6_Weight, NA, NA, NA
ANNAGEBIO.nc timesteps: 721
ANNAGEBIO.nc boxes: 60
ANNAGEBIO.nc layers: 8
Annual age output is in numbers and weight at each timestep, layer, and box, which is really good news. We should be able to get weight at true age more directly this way than the interpolation I was planning. Also, we may be able to do length estimation differently?
To modify load_nc_annage we need to figure out how many true ages each species has from the biological parameter file using load_biolprm:
biolprm <- load_biolprm(d.name, biol.prm.file)## Warning in `[<-.data.frame`(`*tmp*`, , 3:12, value = structure(list(`1` =
## structure(c(124L, : provided 20 variables to replace 10 variablesfgs <- load_fgs(d.name, functional.groups.file)
#max age should be ages per cohort times the number of age classes (<10 for some NOBA spp)
names(biolprm$agespercohort) <- c("code", "ageperagecl")
maxage <- merge(biolprm$maturityogive, biolprm$agespercohort, all.x = T) %>%
select(code, nagecl, ageperagecl) %>%
mutate(maxage = nagecl * ageperagecl) %>%
mutate(name = fgs$Name[match(code, fgs$Code)])
select_groups <- c('Mackerel', 'North_atl_cod', 'Capelin')
for (i in seq_along(select_groups)) {
select_group_ages <- 1:maxage$maxage[match(select_groups[i], maxage$name)]
}Try this load_nc_annage function:
load_nc_annage <- function(dir = getwd(), file_nc, bps, fgs, biolprm, select_groups,
select_variable =
c("Nums", "Weight"),
check_acronyms = TRUE, bboxes = c(0), verbose = FALSE) {
# NOTE: The extraction procedure may look a bit complex...
# A different approach would be to
# create a dataframe for each variable (e.g. GroupAge_Nums)
# and combine all dataframes at the end.
# This alternative approach requires a lot more storage
# and the code wouldn't be vectorized.
# Check input of the nc file
if (tail(strsplit(file_nc, "\\.")[[1]], 1) != "nc") {
stop("The argument for file_nc,", file_nc, "does not end in nc")
}
if (is.null(dir)) {
file.nc <- file_nc
} else {
file.nc <- file.path(dir, file_nc)
}
# Check input of select_variable as only one value is allowed
select_variable <- match.arg(select_variable, several.ok = FALSE)
# Check input structure!
if (check_acronyms) {
active_groups <- as.vector(subset(fgs, IsTurnedOn == 1)$Name)
inactive_groups <- select_groups[which(
!is.element(select_groups, active_groups))]
if (length(inactive_groups) >= 1) {
select_groups <- select_groups[!is.element(select_groups, inactive_groups)]
warning(paste(paste("Some selected groups are not active in the model run.",
"Check 'IsTurnedOn' in fgs\n"),
paste(inactive_groups, collapse = "\n")))
}
if (all(!is.element(select_groups, active_groups))) {
stop(paste("None of the species selected are active in the model run.",
"Check spelling and Check 'IsTurnedOn' in fgs"))
}
}
# Deal with file structures
# Load ATLANTIS output!
at_out <- RNetCDF::open.nc(con = file.nc)
on.exit(RNetCDF::close.nc(at_out))
if (select_variable != "N" & all(is.element(select_groups, bps))) {
stop("The only output for Biomasspools is N.")
} else{
print(paste("Read", file.nc, "successfully"))
}
# 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
# Extract data from the ncdf file
# Create a vector of all potential variable names
# Only use names which are available in the ncdf-file as an
# extraction of missing variables is not possible
# Create vector of available species at the end using search_clean
# This is needed to create species-names later on
# CHANGED hardcoded cohort number to reflect true ages from
# To make the creation of variables as robust as possible
# we introduce different combinations of groups, variable, and cohort
# Only combinations present in the ncdf are used later on
# Loop over select_groups to use the same ordering
names(biolprm$agespercohort) <- c("code", "ageperagecl")
maxage <- merge(biolprm$maturityogive, biolprm$agespercohort, all.x = T) %>%
select(code, nagecl, ageperagecl) %>%
mutate(maxage = nagecl * ageperagecl) %>%
mutate(name = fgs$Name[match(code, fgs$Code)])
search <- list()
for (i in seq_along(select_groups)) {
select_group_ages <- 1:maxage$maxage[match(select_groups[i], maxage$name)]
search[[i]] <- c(
unlist(lapply(paste0(select_groups[i], select_group_ages),
paste0, select_variable)), # GroupCohortVariable
unlist(lapply(paste0(select_groups[i], select_variable),
paste0, select_group_ages)), # GroupVariableCohort
unlist(lapply(paste0(select_groups[i], select_group_ages),
paste, select_variable, sep = "_")), # GroupCohort_Variable
unlist(lapply(paste(select_groups[i], select_variable, sep = "_"),
paste0, select_group_ages)), # Group_VariableCohort
unlist(lapply(paste(select_groups[i], select_group_ages, sep = "_"),
paste, select_variable, sep = "_")), # Group_Cohort_Variable
unlist(lapply(paste(select_groups[i], select_variable, sep = "_"),
paste, select_group_ages, sep = "_")), # Group_Variable_Cohort
unlist(paste0(select_groups[i], select_variable)), # GroupVariable
unlist(paste(select_groups[i], select_variable,
sep = "_")) # Group_Variable
)
search[[i]] <- search[[i]][is.element(search[[i]], var_names_ncdf)]
search[[i]] <- unique(search[[i]])
}
search_clean <- do.call(c, search)
# If the combination of select_groups and select_variable ends up not being found.
if (length(search_clean) == 0) return(0)
at_data <- lapply(search_clean, RNetCDF::var.get.nc, ncfile = at_out)
# Get final species and number of ageclasses per species
final_species <- select_groups[sapply(
lapply(select_groups, grepl, x = search_clean), any)]
final_agecl <- maxage$maxage[
sapply(final_species, function(x) which(x == maxage$name))]
num_layers <- RNetCDF::var.get.nc(ncfile = at_out, variable = "numlayers")[, 1]
# add sediment layer!
num_layers <- num_layers + ifelse(num_layers == 0, 0, 1)
# Create an array of layerids.
# Every entry in the array indicates if a layer is present (= 1) or not (= 0).
# Boxes without layers (= islands) have only 0s as id,
# used later on to remove data from non-existent layers!
# By default output should be 0 in those layers.
# Layers in boundary boxes are set to 0 if bboxes is anything other than NULL!
# Applying a boolean array to an array results in a vector!
for (i in seq_along(num_layers)) {
if (i == 1) layerid <- array(dim = c(n_layers, n_boxes))
if (num_layers[i] == 0) {
layerid[, i] <- 0
} else {
if (!is.null(bboxes) & is.element((i - 1), bboxes)) {
layerid[, i] <- 0
} else {
layerid[, i] <- c(rep(0, times = n_layers - num_layers[i]),
rep(1, times = num_layers[i]))
}
}
}
# Create vectors for polygons and layers
# Each vector has the length equal to one time-step
# All data from islands and non-existent layers is removed
# Therefore the length of these
# vectors is equal for each extracted variable
boxes <- 0:(n_boxes - 1)
# Remove islands and boundary boxes
island_ids <- num_layers == 0
if (!is.null(bboxes)) {
boundary_ids <- is.element(boxes, bboxes)
island_ids <- island_ids | boundary_ids
}
boxes <- boxes[!island_ids]
num_layers <- num_layers[!island_ids]
polygons <- rep(boxes, times = num_layers)
layers <- sapply(num_layers[num_layers != 0] - 2,
function(x) c(seq(x, from = 0, by = 1), n_layers - 1))
if (any(sapply(layers, length) != num_layers[num_layers != 0])) {
stop("Number of layers incorrect.")
}
layers <- do.call(c, layers)
if (length(polygons) != length(layers)) {
stop("Number of polygons and layers do not match.")
}
# In the following section the data is transformed to a long dataframe
# I haven't found any solution to vectorize the creation of the dataframe
# columns (species, age, polygons,...)
# when data from 2d and 3d arrays
# (e.g. select_variable = "N" all biomasspools are only present in the
# sediment layer.) are read in simultaneously.
# Therefore the current "messy" solution splits the data
# in 2 subpopulations: 2d-data and 3d-data
at_data3d <- at_data[which(sapply(at_data, function(x) length(dim(x))) == 3)]
at_data2d <- at_data[which(sapply(at_data, function(x) length(dim(x))) == 2)]
int_fs <- final_species
int_fa <- final_agecl
if (length(at_data3d) >= 1) {
# Remove biomasspools if selected variable is "N"!
if (select_variable == "N") {
int_fs <- final_species[!is.element(final_species, bps)]
int_fa <- final_agecl[!is.element(final_species, bps)]
# Note this only works if age-structured vertebrates have 10 ageclasses
int_fa[int_fa == 10] <- 1
}
for (i in seq_along(at_data3d)) {# for loop over all variables
if (i == 1) result3d <- list()
for (j in 1:n_timesteps) {# loop over timesteps
if (j == 1) values <- array(dim = c(length(layers), n_timesteps))
values[, j] <- at_data3d[[i]][,, j][layerid == 1]
}
result3d[[i]] <- as.vector(values)
}
result3d <- data.frame(
species = unlist(sapply(
X = mapply(FUN = rep, x = int_fs, each = int_fa, SIMPLIFY = FALSE,
USE.NAMES = FALSE),
FUN = rep, each = length(layers) * n_timesteps, simplify = FALSE)),
agecl = unlist(sapply(
X = sapply(X = int_fa, FUN = seq, from = 1, by = 1, simplify = FALSE,
USE.NAMES = FALSE),
FUN = rep, each = length(layers) * n_timesteps, simplify = FALSE)),
polygon = unlist(sapply(
X = n_timesteps * int_fa, FUN = rep, x = polygons, simplify = F,
USE.NAMES = FALSE)),
layer = unlist(sapply(
X = n_timesteps * int_fa, FUN = rep, x = layers, simplify = FALSE,
USE.NAMES = FALSE)),
time = unlist(sapply(
X = int_fa, FUN = rep, x = rep(0:(n_timesteps - 1), each = length(layers)),
simplify = FALSE, USE.NAMES = FALSE)),
atoutput = do.call(c, result3d),
stringsAsFactors = FALSE)
}
if (length(at_data2d) >= 1) {
# Only select biomasspools if selected variable is "N"!
if (select_variable == "N") {
int_fs <- final_species[is.element(final_species, bps)]
int_fa <- final_agecl[is.element(final_species, bps)]
}
# age-structured invert groups are combined in ncdf file!
if (select_variable == "Grazing") int_fa <- 1
for (i in seq_along(at_data2d)) {# for loop over all variables!
if (i == 1) result2d <- list()
for (j in 1:n_timesteps) {# loop over timesteps
if (j == 1) values <- array(dim = c(length(boxes), n_timesteps))
values[, j] <- at_data2d[[i]][, j][boxes + 1]
}
result2d[[i]] <- as.vector(values)
}
# Order of the data in value column = "atoutput".
# 1. species --> rep each with the number of
# ageclasses and n_timesteps * boxes
# 2. age --> rep each (1:10 for each species) with n_timesteps * boxes
# 3. timestep --> rep each timestep (1:n_timesteps)
# with the number of boxes and final_agecl
# (num cohorts per species)
# 4. polygon --> rep boxes times n_timesteps * final_agecl
# (num cohorts per species)
result2d <- data.frame(species = unlist(sapply(
X = mapply(FUN = rep, x = int_fs, each = int_fa, SIMPLIFY = FALSE,
USE.NAMES = FALSE),
FUN = rep, each = length(boxes) * n_timesteps, simplify = FALSE)),
agecl = unlist(sapply(X = sapply(X = int_fa, FUN = seq, from = 1,
by = 1, simplify = FALSE, USE.NAMES = FALSE),
FUN = rep, each = length(boxes) * n_timesteps, simplify = FALSE)),
polygon = unlist(sapply(X = n_timesteps * int_fa,
FUN = rep, x = boxes, simplify = FALSE, USE.NAMES = FALSE)),
time = unlist(sapply(X = int_fa, FUN = rep, x = rep(0:(n_timesteps - 1),
each = length(boxes)), simplify = FALSE, USE.NAMES = FALSE)),
atoutput = do.call(c, result2d),
stringsAsFactors = F)
if (select_variable == "N") result2d$layer <- n_layers - 1
}
# Combine dataframes if necessary!
if (all(sapply(lapply(at_data, dim), length) == 3) & select_variable != "N") {
result <- result3d
}
if (all(sapply(lapply(at_data, dim), length) == 2) & select_variable != "N") {
result <- result2d
}
if (select_variable == "N") {
if (length(at_data2d) >= 1 & length(at_data3d) == 0) result <- result2d
if (length(at_data2d) == 0 & length(at_data3d) >= 1) result <- result3d
if (length(at_data2d) >= 1 & length(at_data3d) >= 1) {
result <- rbind(result2d, result3d)
}
}
# Remove min_pools if existent (well, there always are min pools... ;)).
min_pools <- is.element(result$atoutput, c(0, 1e-08, 1e-16))
if (length(min_pools) > 0) {
# exclude 1st timestep and sediment layer from calculation
print_min_pools <- sum(min_pools) -
length(result[min_pools & result$time == 1, 1]) -
length(result[min_pools & result$time > 1 & result$layer == 7, 1])
if (print_min_pools > 0 & verbose) {
warning(paste0(round(print_min_pools/dim(result)[1] * 100),
"% of ", select_variable, " are true min-pools (0, 1e-08, 1e-16)"))
}
result <- result[!min_pools, ]
}
# Remove non-existent layers.
# WARNING: Biomass is build up (very few) in sediment layer for
# NON sediment groups (e.g. baleen whales)
# Therefore, I subset all data from that layer for non biomass groups and
# groups which cannot penetrate into the sediment!
# UPDATE: Doesn't work with layers as species are not distributed through
# the whole water column and do not appear in
# every polygon.
# Sum up N for invert cohorts if invert cohorts are present!
# NOTE: invert cohorts of size 10 are not considered!
if (select_variable == "N" & any(final_agecl != 10 & final_agecl > 1)) {
result <- result %>%
dplyr::group_by(polygon, layer, species, time) %>%
dplyr::summarise(atoutput = sum(atoutput))
}
return(result)
}Does it work? This snippet will go into run_truth if it does:
file_nc <- ("outputnordic_runresults_01ANNAGEBIO.nc")
bps <- load_bps(dir = d.name, fgs = functional.groups.file, file_init = biomass.pools.file)
allboxes <- load_box(dir = d.name, file_bgm = box.file)
boxes <- get_boundary(allboxes)
numsage <- load_nc_annage(dir = d.name,
file_nc = file_nc,
bps = bps,
fgs = fgs,
biolprm = biolprm,
select_groups = select_groups,
select_variable = "Nums",
check_acronyms = TRUE,
bboxes = boxes,
verbose = TRUE)## [1] "Read /Users/sgaichas/Documents/0_Data/Atlantis/Poseidon/poseidon-dev/atlantisoutput/NOBA_march_2020/outputnordic_runresults_01ANNAGEBIO.nc successfully"## Warning in load_nc_annage(dir = d.name, file_nc = file_nc, bps = bps, fgs =
## fgs, : 1% of Nums are true min-pools (0, 1e-08, 1e-16) #if(verbose) message("Numbers read in from ANNAGEBIO.")
weightage <- load_nc_annage(dir = d.name,
file_nc = file_nc,
bps = bps,
fgs = fgs,
biolprm = biolprm,
select_groups = select_groups,
select_variable = "Weight",
check_acronyms = TRUE,
bboxes = boxes,
verbose = TRUE)## [1] "Read /Users/sgaichas/Documents/0_Data/Atlantis/Poseidon/poseidon-dev/atlantisoutput/NOBA_march_2020/outputnordic_runresults_01ANNAGEBIO.nc successfully"## Warning in load_nc_annage(dir = d.name, file_nc = file_nc, bps = bps, fgs =
## fgs, : 81% of Weight are true min-pools (0, 1e-08, 1e-16) #if(verbose) message("Weight read in from ANNAGEBIO.")Need to aggreagate these to see if they look reasonable; this code snippet from om_comps with numsage output subsituted for nums.
omlist_ss<-NOBAom_cod
source(here("/config/NOBA2Config.R"))
source(here("/config/omdimensions.R"))
source(here("/config/codsurvey.R"))
source(here("/config/codfishery.R"))
n_reps <- 1
#numbers based fishery independent survey for age and length comps
# same user specifications as indices
survey_N <- atlantisom::create_survey(dat = numsage,
time = survtime,
species = survspp,
boxes = survboxes,
effic = surveffic,
selex = survselex)
#Sample fish for age composition
# if we want replicates for obs error this sample function will generate them
age_comp_data <- list()
for(i in 1:n_reps){
age_comp_data[[i]] <- atlantisom::sample_fish(survey_N, surveffN)
}
# save age comps
#if(save){
saveRDS(age_comp_data, file.path(d.name, paste0(scenario.name, "survObsFullAgeComp.rds")))
#}Plot full age comps, have 20 years for cod now, no longer age class but annual ages.
#read in comp data
age_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsFullAgeComp.rds")))
age_comp_data <- age_comp_data[[1]]
#fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
#fish_age_comp <- fish_age_comp[[1]]
#add this to om_indices function so that this has years when read in
#fish_age_comp$time <- as.integer(floor(fish_age_comp$time/fstepperyr))
Natage <- filter(age_comp_data, time %in% c(55:175))
Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste(scenario.name,
Natage$species))
Natageplot + facet_wrap(~time/stepperyr, ncol=6, scales="free_y")
What is the weight output? Documentation only mentions numbers output. Looks like average weight at age for each time and polygon? If so, need to do weighted average weight using the numbers at age in each polygon?
This output is described in Atlantis code, file atlantis/atphysics/atagetracerIO.c, function writeBMAnnAgeBioData:
if (propit == out_nums_id) { // Numbers case
for (b = 0; b < bm->nbox; b++) {
val[b] = 0.0;
for (k = 0; k < bm->wcnz; k++) { // Sum over all individuals in the box
val[b] += bm->wctr[b][k][id] * FunctGroupArray[sp].boxPopRatio[b][k][cohort][ai];
}
}
} else { // Size data - get weighted average of values per layer if numbers > min_dens
for (b = 0; b < bm->nbox; b++) {
val[b] = 0.0;
tot = 0.0;
for (k = 0; k < bm->wcnz; k++) { // Sum over all individuals in the box
if (bm->wctr[b][k][id] > bm->min_dens) {
val[b] += bm->wctr[b][k][id] * (bm->wctr[b][k][sn] + bm->wctr[b][k][rn]) * 5.7 * mg_2_g * FunctGroupArray[sp].boxPopRatio[b][k][cohort][ai];
tot += bm->wctr[b][k][id];
}
}
if (tot > 0.0)
val[b] /= tot;
}
}So weight output appears to be average weight at age in grams?
annwtage <- weightage %>%
group_by(species, agecl, time) %>%
summarise(wtage = mean(atoutput)) %>%
filter(time %in% (survtime))
boxwtage <- weightage %>%
mutate(wtage = atoutput) %>%
filter(species %in% c("North_atl_cod"),
polygon %in% c(1,7,24,30,38,45),
agecl %in% c(16:20))
wtagedat <- boxwtage
wageplot <- ggplot(wtagedat, aes(time, wtage)) +
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, " "))
wageplot + facet_wrap(c("polygon"), scales="free_y")## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
So, to modify run_truth we need to add the annage output but not replace the regular output in case users want output from non-age structured groups (inverts, etc.). Also need to see if lengths can be generaged from the annage outputs, if not need to keep the full set of outputs from the .nc file for all species to generate them.
UPDATE: there is a difference between the annage weight at age and that calculated using ‘calc_age2length()’. While the results from the standard nc output look reasonable for NOBA cod, the annage weight at age is half that for the standard (and half what real cod would be). More worrisome, the oldest age classes have smaller average weight at age than some younger ones. I cannot think of why this would be. Until we are sure the annage weight output is correct, we will instead interpolate the weight at age output of ’calc_age2length" for use as assessment input.
Looks like the annagecatch output is only numbers and not weights. We don’t get weight at age from fisheries usually so this is fine.
# taken from existing load_nc
file.nc <- file.path(d.name, "outputnordic_runresults_01ANNAGECATCH.nc")
# Load ATLANTIS output!
at_out <- RNetCDF::open.nc(con = file.nc)
# 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
RNetCDF::close.nc(at_out)ANNAGECATCH.nc variable names (first hundred): volume, nominal_dz, dz, topk, numlayers, Redfish_other1_Catch_dtrawlNCO, Redfish_other1_Discard_dtrawlNCO, Redfish_other2_Catch_dtrawlNCO, Redfish_other2_Discard_dtrawlNCO, Redfish_other3_Catch_dtrawlNCO, Redfish_other3_Discard_dtrawlNCO, Redfish_other4_Catch_dtrawlNCO, Redfish_other4_Discard_dtrawlNCO, Redfish_other5_Catch_dtrawlNCO, Redfish_other5_Discard_dtrawlNCO, Redfish_other6_Catch_dtrawlNCO, Redfish_other6_Discard_dtrawlNCO, Redfish_other7_Catch_dtrawlNCO, Redfish_other7_Discard_dtrawlNCO, Redfish_other8_Catch_dtrawlNCO, Redfish_other8_Discard_dtrawlNCO, Redfish_other9_Catch_dtrawlNCO, Redfish_other9_Discard_dtrawlNCO, Redfish_other10_Catch_dtrawlNCO, Redfish_other10_Discard_dtrawlNCO, Redfish_other11_Catch_dtrawlNCO, Redfish_other11_Discard_dtrawlNCO, Redfish_other12_Catch_dtrawlNCO, Redfish_other12_Discard_dtrawlNCO, Redfish_other13_Catch_dtrawlNCO, Redfish_other13_Discard_dtrawlNCO, Redfish_other14_Catch_dtrawlNCO, Redfish_other14_Discard_dtrawlNCO, Redfish_other15_Catch_dtrawlNCO, Redfish_other15_Discard_dtrawlNCO, Redfish_other16_Catch_dtrawlNCO, Redfish_other16_Discard_dtrawlNCO, Redfish_other17_Catch_dtrawlNCO, Redfish_other17_Discard_dtrawlNCO, Redfish_other18_Catch_dtrawlNCO, Redfish_other18_Discard_dtrawlNCO, Redfish_other19_Catch_dtrawlNCO, Redfish_other19_Discard_dtrawlNCO, Redfish_other20_Catch_dtrawlNCO, Redfish_other20_Discard_dtrawlNCO, Redfish_other21_Catch_dtrawlNCO, Redfish_other21_Discard_dtrawlNCO, Redfish_other22_Catch_dtrawlNCO, Redfish_other22_Discard_dtrawlNCO, Redfish_other23_Catch_dtrawlNCO, Redfish_other23_Discard_dtrawlNCO, Redfish_other24_Catch_dtrawlNCO, Redfish_other24_Discard_dtrawlNCO, Redfish_other25_Catch_dtrawlNCO, Redfish_other25_Discard_dtrawlNCO, Redfish_other26_Catch_dtrawlNCO, Redfish_other26_Discard_dtrawlNCO, Redfish_other27_Catch_dtrawlNCO, Redfish_other27_Discard_dtrawlNCO, Redfish_other28_Catch_dtrawlNCO, Redfish_other28_Discard_dtrawlNCO, Redfish_other29_Catch_dtrawlNCO, Redfish_other29_Discard_dtrawlNCO, Redfish_other30_Catch_dtrawlNCO, Redfish_other30_Discard_dtrawlNCO, Redfish_other31_Catch_dtrawlNCO, Redfish_other31_Discard_dtrawlNCO, Redfish_other32_Catch_dtrawlNCO, Redfish_other32_Discard_dtrawlNCO, Redfish_other33_Catch_dtrawlNCO, Redfish_other33_Discard_dtrawlNCO, Redfish_other34_Catch_dtrawlNCO, Redfish_other34_Discard_dtrawlNCO, Redfish_other35_Catch_dtrawlNCO, Redfish_other35_Discard_dtrawlNCO, Redfish_other36_Catch_dtrawlNCO, Redfish_other36_Discard_dtrawlNCO, Redfish_other37_Catch_dtrawlNCO, Redfish_other37_Discard_dtrawlNCO, Redfish_other38_Catch_dtrawlNCO, Redfish_other38_Discard_dtrawlNCO, Redfish_other39_Catch_dtrawlNCO, Redfish_other39_Discard_dtrawlNCO, Redfish_other40_Catch_dtrawlNCO, Redfish_other40_Discard_dtrawlNCO, Green_halibut1_Catch_dtrawlGRH, Green_halibut1_Discard_dtrawlGRH, Green_halibut2_Catch_dtrawlGRH, Green_halibut2_Discard_dtrawlGRH, Green_halibut3_Catch_dtrawlGRH, Green_halibut3_Discard_dtrawlGRH, Green_halibut4_Catch_dtrawlGRH, Green_halibut4_Discard_dtrawlGRH, Green_halibut5_Catch_dtrawlGRH, Green_halibut5_Discard_dtrawlGRH, Green_halibut6_Catch_dtrawlGRH, Green_halibut6_Discard_dtrawlGRH, Green_halibut7_Catch_dtrawlGRH, Green_halibut7_Discard_dtrawlGRH, Green_halibut8_Catch_dtrawlGRH
ANNAGECATCH.nc variable names (last hundred): Norwegian_ssh8_Catch_pseineSSH, Norwegian_ssh8_Discard_pseineSSH, Norwegian_ssh9_Catch_pseineSSH, Norwegian_ssh9_Discard_pseineSSH, Norwegian_ssh10_Catch_pseineSSH, Norwegian_ssh10_Discard_pseineSSH, Norwegian_ssh11_Catch_pseineSSH, Norwegian_ssh11_Discard_pseineSSH, Norwegian_ssh12_Catch_pseineSSH, Norwegian_ssh12_Discard_pseineSSH, Norwegian_ssh13_Catch_pseineSSH, Norwegian_ssh13_Discard_pseineSSH, Norwegian_ssh14_Catch_pseineSSH, Norwegian_ssh14_Discard_pseineSSH, Norwegian_ssh15_Catch_pseineSSH, Norwegian_ssh15_Discard_pseineSSH, Norwegian_ssh16_Catch_pseineSSH, Norwegian_ssh16_Discard_pseineSSH, Norwegian_ssh17_Catch_pseineSSH, Norwegian_ssh17_Discard_pseineSSH, Norwegian_ssh18_Catch_pseineSSH, Norwegian_ssh18_Discard_pseineSSH, Norwegian_ssh19_Catch_pseineSSH, Norwegian_ssh19_Discard_pseineSSH, Norwegian_ssh20_Catch_pseineSSH, Norwegian_ssh20_Discard_pseineSSH, North_atl_cod1_Catch_dtrawlNCO, North_atl_cod1_Discard_dtrawlNCO, North_atl_cod2_Catch_dtrawlNCO, North_atl_cod2_Discard_dtrawlNCO, North_atl_cod3_Catch_dtrawlNCO, North_atl_cod3_Discard_dtrawlNCO, North_atl_cod4_Catch_dtrawlNCO, North_atl_cod4_Discard_dtrawlNCO, North_atl_cod5_Catch_dtrawlNCO, North_atl_cod5_Discard_dtrawlNCO, North_atl_cod6_Catch_dtrawlNCO, North_atl_cod6_Discard_dtrawlNCO, North_atl_cod7_Catch_dtrawlNCO, North_atl_cod7_Discard_dtrawlNCO, North_atl_cod8_Catch_dtrawlNCO, North_atl_cod8_Discard_dtrawlNCO, North_atl_cod9_Catch_dtrawlNCO, North_atl_cod9_Discard_dtrawlNCO, North_atl_cod10_Catch_dtrawlNCO, North_atl_cod10_Discard_dtrawlNCO, North_atl_cod11_Catch_dtrawlNCO, North_atl_cod11_Discard_dtrawlNCO, North_atl_cod12_Catch_dtrawlNCO, North_atl_cod12_Discard_dtrawlNCO, North_atl_cod13_Catch_dtrawlNCO, North_atl_cod13_Discard_dtrawlNCO, North_atl_cod14_Catch_dtrawlNCO, North_atl_cod14_Discard_dtrawlNCO, North_atl_cod15_Catch_dtrawlNCO, North_atl_cod15_Discard_dtrawlNCO, North_atl_cod16_Catch_dtrawlNCO, North_atl_cod16_Discard_dtrawlNCO, North_atl_cod17_Catch_dtrawlNCO, North_atl_cod17_Discard_dtrawlNCO, North_atl_cod18_Catch_dtrawlNCO, North_atl_cod18_Discard_dtrawlNCO, North_atl_cod19_Catch_dtrawlNCO, North_atl_cod19_Discard_dtrawlNCO, North_atl_cod20_Catch_dtrawlNCO, North_atl_cod20_Discard_dtrawlNCO, Capelin1_Catch_pseineCAP, Capelin1_Discard_pseineCAP, Capelin2_Catch_pseineCAP, Capelin2_Discard_pseineCAP, Capelin3_Catch_pseineCAP, Capelin3_Discard_pseineCAP, Capelin4_Catch_pseineCAP, Capelin4_Discard_pseineCAP, Capelin5_Catch_pseineCAP, Capelin5_Discard_pseineCAP, Snow_crab1_Catch_trapKCR, Snow_crab1_Catch_trapSCR, Snow_crab1_Discard_trapKCR, Snow_crab1_Discard_trapSCR, Snow_crab2_Catch_trapKCR, Snow_crab2_Catch_trapSCR, Snow_crab2_Discard_trapKCR, Snow_crab2_Discard_trapSCR, Snow_crab3_Catch_trapKCR, Snow_crab3_Catch_trapSCR, Snow_crab3_Discard_trapKCR, Snow_crab3_Discard_trapSCR, Snow_crab4_Catch_trapKCR, Snow_crab4_Catch_trapSCR, Snow_crab4_Discard_trapKCR, Snow_crab4_Discard_trapSCR, Snow_crab5_Catch_trapKCR, Snow_crab5_Catch_trapSCR, Snow_crab5_Discard_trapKCR, Snow_crab5_Discard_trapSCR, Snow_crab6_Catch_trapKCR, Snow_crab6_Catch_trapSCR, Snow_crab6_Discard_trapKCR, Snow_crab6_Discard_trapSCR
ANNAGECATCH.nc timesteps: 721
ANNAGECATCH.nc boxes: 60
ANNAGECATCH.nc layers: 8
Catch file contains “IsFished” species discard at age and catch at age by fishery in numbers. So now modifying load_nc_annage to get the fishery specific outputs, which requires additional atlantis fishery setup files and additions to the config files above.
New function load_fisheries() needed to load the fishery specification file (a .csv), similar to load_fgs():
load_fisheries <- function(dir = getwd(), file_fish){
if (is.null(dir)) {
file.fish <- file_fish
} else {
file.fish <- file.path(dir, file_fish)
}
result <- read.table(file = file.fish, sep = ",",
header = TRUE, stringsAsFactors = FALSE)
return(result)
}Call load_fisheries() within load_nc_annage() to rename select_variable as used in ANNAGECATCH.nc:
select_variable <- c("Catch", "Discard")
dir <- here("atlantisoutput","NOBA_march_2020")
file_fish <- "NoBAFisheries.csv"
#starting at line 115 in current load_nc_annage
if(select_variable %in% c("Catch", "Discard")){
#read in fleet names
fleetnames <- load_fisheries(dir = dir, file_fish = file_fish)
svfish <- c()
for(i in select_variable){
sv <-paste0(i, "_", fleetnames$Code)
svfish <- c(svfish, sv)
}
select_variable <- svfish
}Also added if statements to add fleets for fisheries outputs, if statements to separate output for annage and annagecatch.
Test new load_nc_annage():
# source(here("R/load_nc_annage.R"))
# source(here("R/load_fisheries.R"))
source(here("config/NOBA2Config.R"))
NOBAom_cod <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))
fgs <- NOBAom_cod$funct.group_ss
biolprm <- NOBAom_cod$biol
#select_groups <- c('Mackerel', 'North_atl_cod', 'Capelin')
select_groups <- 'North_atl_cod'
bps <- load_bps(dir = d.name, fgs = functional.groups.file, file_init = biomass.pools.file)
allboxes <- load_box(dir = d.name, file_bgm = box.file)
boxes <- get_boundary(allboxes)
file_nc <- "outputnordic_runresults_01ANNAGECATCH.nc"
file_fish <- "NoBAFisheries.csv"
#for testing line by line
bboxes <- boxes
dir <- d.name
select_variable <- "Catch"
fishcatchage <- load_nc_annage(dir = d.name,
file_nc = file_nc,
file_fish = file_fish,
bps = bps,
fgs = fgs,
biolprm = biolprm,
select_groups = select_groups,
select_variable = "Catch",
check_acronyms = TRUE,
bboxes = boxes,
verbose = TRUE)## [1] "Read /Users/sgaichas/Documents/0_Data/Atlantis/Poseidon/poseidon-dev/atlantisoutput/NOBA_march_2020/outputnordic_runresults_01ANNAGECATCH.nc successfully"## Warning in if (select_variable == "N") {: the condition has length > 1 and only
## the first element will be used## Warning in if (select_variable == "Grazing") int_fa <- 1: the condition has
## length > 1 and only the first element will be used## Warning in if (all(sapply(lapply(at_data, dim), length) == 3) &
## select_variable != : the condition has length > 1 and only the first element
## will be used## Warning in if (all(sapply(lapply(at_data, dim), length) == 2) &
## select_variable != : the condition has length > 1 and only the first element
## will be used## Warning in if (select_variable == "N") {: the condition has length > 1 and only
## the first element will be used## Warning in load_nc_annage(dir = d.name, file_nc = file_nc, file_fish =
## file_fish, : 72% of Catch_pseineSSH are true min-pools (0, 1e-08, 1e-16)72% of
## Catch_pseineBWH are true min-pools (0, 1e-08, 1e-16)72% of Catch_pseineMAC are
## true min-pools (0, 1e-08, 1e-16)72% of Catch_pseineCAP are true min-pools (0,
## 1e-08, 1e-16)72% of Catch_dtrawlNCO are true min-pools (0, 1e-08, 1e-16)72% of
## Catch_dtrawlHAD are true min-pools (0, 1e-08, 1e-16)72% of Catch_dtrawlSAI are
## true min-pools (0, 1e-08, 1e-16)72% of Catch_dtrawlGRH are true min-pools (0,
## 1e-08, 1e-16)72% of Catch_dtrawlPWN are true min-pools (0, 1e-08, 1e-16)72% of
## Catch_dlineNCO are true min-pools (0, 1e-08, 1e-16)72% of Catch_dlineHAD are
## true min-pools (0, 1e-08, 1e-16)72% of Catch_dlineSAI are true min-pools (0,
## 1e-08, 1e-16)72% of Catch_dlineGRH are true min-pools (0, 1e-08, 1e-16)72% of
## Catch_dseineNCO are true min-pools (0, 1e-08, 1e-16)72% of Catch_dseineHAD are
## true min-pools (0, 1e-08, 1e-16)72% of Catch_dseineSAI are true min-pools (0,
## 1e-08, 1e-16)72% of Catch_dseineGRH are true min-pools (0, 1e-08, 1e-16)72%
## of Catch_netNCO are true min-pools (0, 1e-08, 1e-16)72% of Catch_netHAD are
## true min-pools (0, 1e-08, 1e-16)72% of Catch_netSAI are true min-pools (0,
## 1e-08, 1e-16)72% of Catch_netGRH are true min-pools (0, 1e-08, 1e-16)72% of
## Catch_trapKCR are true min-pools (0, 1e-08, 1e-16)72% of Catch_trapSCR are
## true min-pools (0, 1e-08, 1e-16)72% of Catch_cullHOS are true min-pools (0,
## 1e-08, 1e-16)72% of Catch_cullHAS are true min-pools (0, 1e-08, 1e-16)72% of
## Catch_cullMWH are true min-pools (0, 1e-08, 1e-16)72% of Catch_REC are true min-
## pools (0, 1e-08, 1e-16)## Warning in if (select_variable == "N" & any(final_agecl != 10 & final_agecl > :
## the condition has length > 1 and only the first element will be usedfishcatchage[1:100,]## species agecl polygon fleet time atoutput
## 68 North_atl_cod 1 19 dtrawlNCO 1 1.566299e+03
## 70 North_atl_cod 1 21 dtrawlNCO 1 1.726160e+03
## 71 North_atl_cod 1 23 dtrawlNCO 1 3.664902e+03
## 72 North_atl_cod 1 24 dtrawlNCO 1 3.724425e+03
## 73 North_atl_cod 1 25 dtrawlNCO 1 8.758350e+03
## 74 North_atl_cod 1 26 dtrawlNCO 1 3.920000e+03
## 75 North_atl_cod 1 27 dtrawlNCO 1 1.241475e+03
## 76 North_atl_cod 1 28 dtrawlNCO 1 1.113926e+02
## 77 North_atl_cod 1 29 dtrawlNCO 1 5.144468e+02
## 78 North_atl_cod 1 30 dtrawlNCO 1 2.278872e+03
## 79 North_atl_cod 1 31 dtrawlNCO 1 3.078178e+02
## 81 North_atl_cod 1 33 dtrawlNCO 1 3.571366e+03
## 82 North_atl_cod 1 34 dtrawlNCO 1 3.469327e+03
## 83 North_atl_cod 1 35 dtrawlNCO 1 3.146204e+03
## 84 North_atl_cod 1 37 dtrawlNCO 1 4.285639e+03
## 85 North_atl_cod 1 38 dtrawlNCO 1 4.889370e+03
## 86 North_atl_cod 1 39 dtrawlNCO 1 4.889370e+03
## 87 North_atl_cod 1 40 dtrawlNCO 1 4.370672e+03
## 88 North_atl_cod 1 41 dtrawlNCO 1 4.719305e+03
## 89 North_atl_cod 1 42 dtrawlNCO 1 3.171713e+03
## 90 North_atl_cod 1 43 dtrawlNCO 1 5.816225e+03
## 92 North_atl_cod 1 45 dtrawlNCO 1 8.214142e+03
## 94 North_atl_cod 1 47 dtrawlNCO 1 6.675054e+03
## 3568 North_atl_cod 1 19 dtrawlNCO 71 2.207436e+06
## 3570 North_atl_cod 1 21 dtrawlNCO 71 2.432734e+06
## 3571 North_atl_cod 1 23 dtrawlNCO 71 5.165065e+06
## 3572 North_atl_cod 1 24 dtrawlNCO 71 5.248953e+06
## 3573 North_atl_cod 1 25 dtrawlNCO 71 1.234343e+07
## 3574 North_atl_cod 1 26 dtrawlNCO 71 5.524582e+06
## 3575 North_atl_cod 1 27 dtrawlNCO 71 1.749651e+06
## 3576 North_atl_cod 1 28 dtrawlNCO 71 1.569892e+05
## 3577 North_atl_cod 1 29 dtrawlNCO 71 7.250266e+05
## 3578 North_atl_cod 1 30 dtrawlNCO 71 3.211688e+06
## 3579 North_atl_cod 1 31 dtrawlNCO 71 4.338175e+05
## 3581 North_atl_cod 1 33 dtrawlNCO 71 5.033242e+06
## 3582 North_atl_cod 1 34 dtrawlNCO 71 4.889435e+06
## 3583 North_atl_cod 1 35 dtrawlNCO 71 4.434047e+06
## 3584 North_atl_cod 1 37 dtrawlNCO 71 6.039891e+06
## 3585 North_atl_cod 1 38 dtrawlNCO 71 6.890748e+06
## 3586 North_atl_cod 1 39 dtrawlNCO 71 6.890748e+06
## 3587 North_atl_cod 1 40 dtrawlNCO 71 6.159730e+06
## 3588 North_atl_cod 1 41 dtrawlNCO 71 6.651070e+06
## 3589 North_atl_cod 1 42 dtrawlNCO 71 4.469998e+06
## 3590 North_atl_cod 1 43 dtrawlNCO 71 8.196994e+06
## 3592 North_atl_cod 1 45 dtrawlNCO 71 1.157646e+07
## 3594 North_atl_cod 1 47 dtrawlNCO 71 9.407369e+06
## 3603 North_atl_cod 1 3 dtrawlNCO 72 3.934830e+06
## 3604 North_atl_cod 1 4 dtrawlNCO 72 1.433358e+06
## 3605 North_atl_cod 1 5 dtrawlNCO 72 6.338402e+06
## 3621 North_atl_cod 1 23 dtrawlNCO 72 3.733815e+06
## 3622 North_atl_cod 1 24 dtrawlNCO 72 3.794865e+06
## 3624 North_atl_cod 1 26 dtrawlNCO 72 4.001732e+06
## 3625 North_atl_cod 1 27 dtrawlNCO 72 3.620414e+06
## 3626 North_atl_cod 1 28 dtrawlNCO 72 3.244843e+05
## 3627 North_atl_cod 1 29 dtrawlNCO 72 1.503507e+06
## 3628 North_atl_cod 1 30 dtrawlNCO 72 6.661252e+06
## 3629 North_atl_cod 1 31 dtrawlNCO 72 8.980941e+05
## 3631 North_atl_cod 1 33 dtrawlNCO 72 1.044068e+07
## 3633 North_atl_cod 1 35 dtrawlNCO 72 3.208086e+06
## 3635 North_atl_cod 1 38 dtrawlNCO 72 8.560537e+06
## 3636 North_atl_cod 1 39 dtrawlNCO 72 6.997901e+06
## 3639 North_atl_cod 1 42 dtrawlNCO 72 3.234250e+06
## 3640 North_atl_cod 1 43 dtrawlNCO 72 5.934050e+06
## 3668 North_atl_cod 1 19 dtrawlNCO 73 8.074563e+06
## 3670 North_atl_cod 1 21 dtrawlNCO 73 9.174135e+06
## 3671 North_atl_cod 1 23 dtrawlNCO 73 2.077905e+07
## 3672 North_atl_cod 1 24 dtrawlNCO 73 2.109960e+07
## 3673 North_atl_cod 1 25 dtrawlNCO 73 4.636747e+07
## 3674 North_atl_cod 1 26 dtrawlNCO 73 2.223648e+07
## 3679 North_atl_cod 1 31 dtrawlNCO 73 2.043659e+06
## 3681 North_atl_cod 1 33 dtrawlNCO 73 2.372728e+07
## 3683 North_atl_cod 1 35 dtrawlNCO 73 1.783748e+07
## 3684 North_atl_cod 1 37 dtrawlNCO 73 2.272006e+07
## 3685 North_atl_cod 1 38 dtrawlNCO 73 2.965904e+07
## 3686 North_atl_cod 1 39 dtrawlNCO 73 2.885959e+07
## 3687 North_atl_cod 1 40 dtrawlNCO 73 2.315061e+07
## 3688 North_atl_cod 1 41 dtrawlNCO 73 2.500531e+07
## 3689 North_atl_cod 1 42 dtrawlNCO 73 1.801271e+07
## 3690 North_atl_cod 1 43 dtrawlNCO 73 3.298128e+07
## 3718 North_atl_cod 1 19 dtrawlNCO 74 9.301161e+06
## 3720 North_atl_cod 1 21 dtrawlNCO 74 1.056777e+07
## 3721 North_atl_cod 1 23 dtrawlNCO 74 2.235636e+07
## 3722 North_atl_cod 1 24 dtrawlNCO 74 2.269972e+07
## 3723 North_atl_cod 1 25 dtrawlNCO 74 5.341110e+07
## 3724 North_atl_cod 1 26 dtrawlNCO 74 2.392054e+07
## 3729 North_atl_cod 1 31 dtrawlNCO 74 1.877019e+06
## 3731 North_atl_cod 1 33 dtrawlNCO 74 2.178410e+07
## 3733 North_atl_cod 1 35 dtrawlNCO 74 1.918984e+07
## 3734 North_atl_cod 1 37 dtrawlNCO 74 2.617144e+07
## 3735 North_atl_cod 1 38 dtrawlNCO 74 2.994836e+07
## 3736 North_atl_cod 1 39 dtrawlNCO 74 2.994836e+07
## 3737 North_atl_cod 1 40 dtrawlNCO 74 2.666740e+07
## 3738 North_atl_cod 1 41 dtrawlNCO 74 2.880384e+07
## 3739 North_atl_cod 1 42 dtrawlNCO 74 1.938060e+07
## 3740 North_atl_cod 1 43 dtrawlNCO 74 3.548023e+07
## 3768 North_atl_cod 1 19 dtrawlNCO 75 4.097839e+06
## 3770 North_atl_cod 1 21 dtrawlNCO 75 4.654757e+06
## 3771 North_atl_cod 1 23 dtrawlNCO 75 9.847523e+06
## 3772 North_atl_cod 1 24 dtrawlNCO 75 9.998832e+06
## 3773 North_atl_cod 1 25 dtrawlNCO 75 2.352656e+07fishdiscage <- load_nc_annage(dir = d.name,
file_nc = file_nc,
file_fish = file_fish,
bps = bps,
fgs = fgs,
biolprm = biolprm,
select_groups = select_groups,
select_variable = "Discard",
check_acronyms = TRUE,
bboxes = boxes,
verbose = TRUE)## [1] "Read /Users/sgaichas/Documents/0_Data/Atlantis/Poseidon/poseidon-dev/atlantisoutput/NOBA_march_2020/outputnordic_runresults_01ANNAGECATCH.nc successfully"## Warning in if (select_variable == "N") {: the condition has length > 1 and only
## the first element will be used## Warning in if (select_variable == "Grazing") int_fa <- 1: the condition has
## length > 1 and only the first element will be used## Warning in if (all(sapply(lapply(at_data, dim), length) == 3) &
## select_variable != : the condition has length > 1 and only the first element
## will be used## Warning in if (all(sapply(lapply(at_data, dim), length) == 2) &
## select_variable != : the condition has length > 1 and only the first element
## will be used## Warning in if (select_variable == "N") {: the condition has length > 1 and only
## the first element will be used## Warning in load_nc_annage(dir = d.name, file_nc = file_nc, file_fish
## = file_fish, : 100% of Discard_pseineSSH are true min-pools (0, 1e-08,
## 1e-16)100% of Discard_pseineBWH are true min-pools (0, 1e-08, 1e-16)100% of
## Discard_pseineMAC are true min-pools (0, 1e-08, 1e-16)100% of Discard_pseineCAP
## are true min-pools (0, 1e-08, 1e-16)100% of Discard_dtrawlNCO are true min-
## pools (0, 1e-08, 1e-16)100% of Discard_dtrawlHAD are true min-pools (0, 1e-08,
## 1e-16)100% of Discard_dtrawlSAI are true min-pools (0, 1e-08, 1e-16)100% of
## Discard_dtrawlGRH are true min-pools (0, 1e-08, 1e-16)100% of Discard_dtrawlPWN
## are true min-pools (0, 1e-08, 1e-16)100% of Discard_dlineNCO are true min-
## pools (0, 1e-08, 1e-16)100% of Discard_dlineHAD are true min-pools (0, 1e-08,
## 1e-16)100% of Discard_dlineSAI are true min-pools (0, 1e-08, 1e-16)100% of
## Discard_dlineGRH are true min-pools (0, 1e-08, 1e-16)100% of Discard_dseineNCO
## are true min-pools (0, 1e-08, 1e-16)100% of Discard_dseineHAD are true min-
## pools (0, 1e-08, 1e-16)100% of Discard_dseineSAI are true min-pools (0, 1e-08,
## 1e-16)100% of Discard_dseineGRH are true min-pools (0, 1e-08, 1e-16)100% of
## Discard_netNCO are true min-pools (0, 1e-08, 1e-16)100% of Discard_netHAD are
## true min-pools (0, 1e-08, 1e-16)100% of Discard_netSAI are true min-pools (0,
## 1e-08, 1e-16)100% of Discard_netGRH are true min-pools (0, 1e-08, 1e-16)100% of
## Discard_trapKCR are true min-pools (0, 1e-08, 1e-16)100% of Discard_trapSCR are
## true min-pools (0, 1e-08, 1e-16)100% of Discard_cullHOS are true min-pools (0,
## 1e-08, 1e-16)100% of Discard_cullHAS are true min-pools (0, 1e-08, 1e-16)100% of
## Discard_cullMWH are true min-pools (0, 1e-08, 1e-16)100% of Discard_REC are true
## min-pools (0, 1e-08, 1e-16)## Warning in if (select_variable == "N" & any(final_agecl != 10 & final_agecl > :
## the condition has length > 1 and only the first element will be usedfishdiscage## [1] species agecl polygon fleet time atoutput
## <0 rows> (or 0-length row.names)Do we want catch by fishery? Yes. No fishery column here, check what it is doing. fishdiscage is empty, make sure there is no discard output from NOBA for these species (confirmed).
Add a check for multiple fisheries per species, if only one then fleet irrelevant.
NOBA has multiple fleets for cod (but they are not in output!), so add fleet column to output.
Update: there is now a fleet column in output which should index correctly if a species is caught by multiple fleets.
Plot catch at age; aggregate over polygons, apply sample_fish:
omlist_ss<-NOBAom_cod
source(here("/config/NOBA2Config.R"))
source(here("/config/omdimensions.R"))
source(here("/config/codsurvey.R"))
source(here("/config/codfishery.R"))
n_reps <- 1
#numbers based fishery independent survey for age and length comps
# same user specifications as indices
#fishery catch at age each observed timestep summed over observed polygons
# catch at age by area and timestep
catch_numbers <- atlantisom::create_fishery_subset(dat = fishcatchage,
time = fishtime,
species = select_groups,
boxes = fishboxes)
# if we want replicates for obs error this sample function will generate them
catch_age_comp <- list()
for(i in 1:n_reps){
catch_age_comp[[i]] <- atlantisom::sample_fish(catch_numbers, fisheffN)
}
# save fishery age comps
#if(save){
saveRDS(catch_age_comp, file.path(d.name, paste0(scenario.name, "fishObsFullAgeComp.rds")))
#}NOTE that the fleet column is lost in the atlantisom::sample_fish function!
Look at only every 10 years output (this is a sample, effN = 1000 for the year, which is 200 per timestep here; see calc below confirming):
#read in comp data
# age_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsFullAgeComp.rds")))
# age_comp_data <- age_comp_data[[1]]
fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]
#add this to om_indices function so that this has years when read in
fish_age_comp$time <- fish_age_comp$time/fstepperyr
Catage <- filter(fish_age_comp, time %in% c(30.0, 40.0, 50.0, 60.0, 70.0, 79.8))
Catageplot <- ggplot(Catage, aes(x=agecl, y=atoutput)) +
geom_point() +
theme_tufte() +
labs(subtitle = paste(scenario.name,
Catage$species))
Catageplot + facet_wrap(~time, ncol=3, scales="free_y")
fish_age_comp %>%
filter(time<31) %>%
summarise(sum(atoutput))## sum(atoutput)
## 1 1000Appears to be working for single fleet per species now.
Quick test of run_truth() with annage option turned on for NOBA:
source(here("config/NOBA2Config.R"))
#Load functional groups
funct.groups <- atlantisom::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
# # load biolprm file
#biol <- atlantisom::load_biolprm(d.name, biol.prm.file)
#
# # load run_prm
# runpar <- atlantisom::load_runprm(d.name, run.prm.file)
#
# # load box info file
# boxpars <- atlantisom::load_box(d.name, box.file)
testtruth <- atlantisom::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,
file_fish = fisheries.file,
verbose = TRUE,
annage = TRUE,
save = TRUE
)After much debugging to properly index multiple fleets per species in catch at annual age output, this produced a run_truth.Rdata file and took about 40 minutes to run.
Now all of the wrappers need to read from annage to produce composition data. om_init() now takes annage flag from the config file. om_species() looks for annage bio and catch output in the truth file and includes it in the species file (adding NULL if it isn’t there). om_comps() looks for annage bio and catch in the species file and adds them to the comps list if there (adding NULL if they aren’t there).
Also still need the weight at age interpolation function. Incorporate this at the om_comps() step. Next rmd!
Running om_init(), om_species(), om_index(), om_comps() in chunk above to make sure the original output is preserved, then plot the annage output above.
Write input data for specific assessment models
We will need to alter this for SAM, SS3, etc. Right now just makes SS3.
THIS IS STILL THE CALIFORNIA CURRENT INFO, TO BE UPDATED
require(r4ss)
omlist_ss<-CC3om_sardine
source(here("/config/CC3Config_sardLWcorr.R"))
source(here("/config/omdimensions.R"))
source(here("/config/sardinesurvey.R"))
source(here("/config/sardinefishery.R"))
#Directory with SS files
model_dir <- "Sardine_SS_files/"
#Name of SS data file
datfile_name <- "sardEM_3_3.dat"
#CVs for length at age, catch, and survey
CVs <- list("lenage"=0.1, "fishery"=0.01, "survey"=0.1)
#Month of survey/fishing
survey_month <- 7
fishing_month <- 1
#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
#read time series data
survObsBiom <- readRDS(file.path(d.name, paste0(scenario.name, "surveyB.rds")))
truecatchbio_ss <- readRDS(file.path(d.name, paste0(scenario.name, "fishCatch.rds")))
survObsBiom <- survObsBiom[[1]]
truecatchbio_ss <- truecatchbio_ss[[1]]
#add this to om_indices function so that this has years when read in
truecatchbio_ss$time <- as.integer(truecatchbio_ss$time/365)
#load dummy dat file
stocksynthesis.data <- r4ss::SS_readdat_3.30(paste0("./inst/extdata/",
model_dir,
datfile_name))
#Test writing CPUE in biomass
stocksynthesis.data <- atlantisom::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))
#test
stocksynthesis.data$CPUE
stocksynthesis.data$catch
#read in comp data
age_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsAgeComp.rds")))
age_comp_data <- age_comp_data[[1]]
fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]
#Get the age bins
age_bin_names <- names(stocksynthesis.data$agecomp)[10:length(names(stocksynthesis.data$agecomp))]
age_bins <- sub("a","",age_bin_names)
require(maditr)
# add dependency on maditr::dcast to atlantisom
age_comp_flat <- atlantisom::reformat_compositions(age_comp_data,
round.places = 4,
comp_type = "agecomp")
## Write age composition data for survey
stocksynthesis.data <- atlantisom::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
#length comps
len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsLenComp.rds")))
fish_len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsLenComp.rds")))
len_comp_data <- len_comp_data[[1]]
fish_len_comp_data <- fish_len_comp_data[[1]]
#add this to om_indices function so that this has years when read in
fish_len_comp_data$time <- as.integer(floor(fish_len_comp_data$time/fstepperyr))
fish_age_comp$time <- as.integer(floor(fish_age_comp$time/fstepperyr))
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)) %>%
group_by(species, agecl, lower.bins, upper.bins, time) %>%
summarise(annnatlength=sum(atoutput)) %>%
rename(atoutput = annnatlength)
} else {
fish_len_comp_anndata <- fish_len_comp_data
}
caal_comp_flat <- atlantisom::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 <- atlantisom::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))
length_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 <- atlantisom::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 <- fish_age_comp %>%
#mutate(yr = floor(time/fstepperyr)) %>%
#group_by(species, agecl, time=as.integer(yr)) %>%
group_by(species, agecl, time) %>%
summarise(annnatage=sum(atoutput)) %>%
rename(atoutput = annnatage)
} else {
fish_age_comp_anndata <- fish_age_comp
}
fish_age_comp_flat <- atlantisom::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 <- atlantisom::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)
head(stocksynthesis.data$agecomp)
#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)