title: "Snow crab Areal unit modelling" author: "Jae S. Choi" toc: true number-sections: true highlight-style: pygments editor: render-on-save: false format: html: code-fold: true html-math-method: katex self-contained: true pdf: pdf-engine: lualatex docx: default
The snow crab biomass index is derived from the convolution of three separate Bayesian spatiotemporal model solutions via posterior simulation. This is completed through the carstm front-end to INLA, used to perform "non-separable" spatial Conditional autocorrelation (CAR) and temporal (AR1) models.
The convolution of all three after posterior simulation is also known as a Hurdle or Delta model.
require(aegis)
require(bio.snowcrab) # loadfunctions("bio.snowcrab")
year.assessment = 2023
yrs = 1999:year.assessment
spec_bio = bio.taxonomy::taxonomy.recode( from="spec", to="parsimonious", tolookup=2526 )
snowcrab_filter_class = "fb" # fishable biomass (including soft-shelled ) "m.mat" "f.mat" "imm"
runlabel= paste( "1999_present", snowcrab_filter_class, sep="_" )
# params for number
pN = snowcrab_parameters(
project_class="carstm",
yrs=yrs,
areal_units_type="tesselation",
carstm_model_label= runlabel,
selection = list(
type = "number",
biologicals=list( spec_bio=spec_bio ),
biologicals_using_snowcrab_filter_class=snowcrab_filter_class
)
)
# params for mean size .. mostly the same as pN
pW = snowcrab_parameters(
project_class="carstm",
yrs=yrs,
areal_units_type="tesselation",
carstm_model_label= runlabel,
selection = list(
type = "meansize",
biologicals=list( spec_bio=spec_bio ),
biologicals_using_snowcrab_filter_class=snowcrab_filter_class
)
)
# params for probability of observation
pH = snowcrab_parameters(
project_class="carstm",
yrs=yrs,
areal_units_type="tesselation",
carstm_model_label= runlabel,
selection = list(
type = "presence_absence",
biologicals=list( spec_bio=spec_bio ),
biologicals_using_snowcrab_filter_class=snowcrab_filter_class
)
)
# areal units upon which carstm will operate ... this is made in 01.snowcrab.r
sppoly=areal_units( p=pN )
pN$space_name = sppoly$AUID
pN$space_id = 1:nrow(sppoly) # must match M$space
pN$time_name = as.character(pN$yrs)
pN$time_id = 1:pN$ny
pN$cyclic_name = as.character(pN$cyclic_levels)
pN$cyclic_id = 1:pN$nw
pW$space_name = sppoly$AUID
pW$space_id = 1:nrow(sppoly) # must match M$space
pW$time_name = as.character(pW$yrs)
pW$time_id = 1:pW$ny
pW$cyclic_name = as.character(pW$cyclic_levels)
pW$cyclic_id = 1:pW$nw
pH$space_name = sppoly$AUID
pH$space_id = 1:nrow(sppoly) # must match M$space
pH$time_name = as.character(pH$yrs)
pH$time_id = 1:pH$ny
pH$cyclic_name = as.character(pH$cyclic_levels)
pH$cyclic_id = 1:pH$nw
# create model data inputs and the output fields
M = snowcrab.db( p=pN, DS="carstm_inputs", sppoly=sppoly, redo=TRUE ) # will redo if not found
# for mapping below, some bathymetry and polygons
additional_features = snowcrab_mapping_features(pN)
With all required parameters defined, the modelling is straightforward. Here we reused previous solution modes ("theta") to speed up solutions.
Each variable is modelled with the same covariates.
# total numbers
sppoly = areal_units( p=pN )
M = snowcrab.db( p=pN, DS="carstm_inputs", sppoly=sppoly ) # will redo if not found
io = which(M$tag=="observations")
ip = which(M$tag=="predictions")
iq = unique( c( which( M$totno > 0), ip ) )
iw = unique( c( which( M$totno > 5), ip ) ) # need a good sample to estimate mean size
# number
res = NULL; gc()
res = carstm_model( p=pN, data=M[ iq, ], sppoly=sppoly,
theta=c( 2.409, 1.874, 0.772, 2.092, -1.490, 5.145, 4.509, 2.178, 5.453, 0.182, 2.742, 0.525, 0.051, 0.779 ),
nposteriors=5000,
posterior_simulations_to_retain=c( "summary", "random_spatial", "predictions"),
family = "poisson",
verbose=TRUE,
# redo_fit=FALSE,
# debug = "summary",
# debug = "predictions",
num.threads="4:3"
)
# mean size
res = NULL; gc()
res = carstm_model( p=pW, data=M[ iw, ], sppoly = sppoly,
theta=c( 6.108, 8.632, 0.883, 2.946, 9.801, 7.265, 10.726, 12.214, 11.849, 9.826, 6.556, 3.456, 5.832, 2.939, 1.625 ),
nposteriors=5000,
posterior_simulations_to_retain=c( "summary", "random_spatial", "predictions"),
family = "gaussian",
verbose=TRUE,
# redo_fit=FALSE,
# debug = "summary",
# control.inla = list( strategy="laplace", int.strategy="eb" ),
num.threads="4:3"
)
# model pa using all data
res = NULL; gc()
res = carstm_model( p=pH, data=M, sppoly=sppoly,
theta = c( 0.926, 1.743, -0.401, 0.705, -2.574, 1.408, 2.390, 3.459, 3.321, -2.138, 3.083, -1.014, 3.558, 2.703 ),
nposteriors=5000,
posterior_simulations_to_retain=c( "summary", "random_spatial", "predictions"),
family = "binomial", # "binomial", # "nbinomial", "betabinomial", "zeroinflatedbinomial0" , "zeroinflatednbinomial0"
verbose=TRUE,
#redo_fit=FALSE,
# debug = "summary",
# control.family=list(control.link=list(model="logit")), # default for binomial .. no need to specify
# control.inla = list( strategy="laplace", int.strategy="eb" ),
num.threads="4:3"
)
# end spatiotemporal model
# some maps and plots
for (vns in c( "number", "meansize", "habitat") ) {
#vns ="number"
#vns ="meansize"
#vns ="habitat"
if ( vns=="number" ) {
p=pN
ylab = "Number"
outputdir = file.path( p$modeldir, p$carstm_model_label, "predicted.numerical.densities" )
if ( !file.exists(outputdir)) dir.create( outputdir, recursive=TRUE, showWarnings=FALSE )
fn_root_prefix = "Predicted_numerical_abundance"
fn_root = "Predicted_numerical_abundance_persistent_spatial_effect"
outfilename = file.path( outputdir, paste(fn_root, "png", sep=".") )
title= paste( snowcrab_filter_class, "Number; no./m^2" )
}
if ( vns=="meansize") {
p=pW
ylab = "Mean weight"
outputdir = file.path( p$modeldir, p$carstm_model_label, "predicted.meansize" )
if ( !file.exists(outputdir)) dir.create( outputdir, recursive=TRUE, showWarnings=FALSE )
fn_root_prefix = "Predicted_meansize"
fn_root = "Predicted_meansize_persistent_spatial_effect"
outfilename = file.path( outputdir, paste(fn_root, "png", sep=".") )
title= paste( snowcrab_filter_class, "Mean weight; kg" )
}
if ( vns=="habitat" ) {
p=pH
ylab = "Probability"
outputdir = file.path( p$modeldir, p$carstm_model_label, "predicted.presence_absence" )
if ( !file.exists(outputdir)) dir.create( outputdir, recursive=TRUE, showWarnings=FALSE )
fn_root_prefix = "Predicted_presence_absence"
fn_root = "Predicted_presence_absence_persistent_spatial_effect"
outfilename = file.path( outputdir, paste(fn_root, "png", sep=".") )
title= paste( snowcrab_filter_class, "Probability")
# # to compute habitat prob
# sims = carstm_posterior_simulations( pH=pH, pa_threshold=0.05, qmax=0.95 )
# SM = aggregate_simulations(
# sims=sims,
# sppoly=sppoly,
# fn=carstm_filenames( pN, returnvalue="filename", fn="aggregated_timeseries" ),
# yrs=pN$yrs,
# method="mean",
# redo=TRUE
# )
# outputdir = file.path( carstm_filenames( pN, returnvalue="output_directory"), "aggregated_habitat_timeseries" )
# ylabel = "Habitat probability"
# fn_ts = "habitat_M0.png"
# vn = paste("habitat", "predicted", sep=".")
# outputdir2 = file.path( carstm_filenames( pN, returnvalue="output_directory"), "predicted_habitat" )
}
res = carstm_model( p=p, DS="carstm_modelled_summary", sppoly = sppoly ) # to load currently saved results
# maps
vn = c( "random", "space", "combined" )
toplot = carstm_results_unpack( res, vn )
brks = pretty( quantile(toplot[,"mean"], probs=c(0,0.975), na.rm=TRUE ) )
carstm_map( res=res, vn=vn,
sppoly = sppoly,
breaks = brks,
palette="-RdYlBu",
plot_elements="",
additional_features=additional_features,
legend.position=c( 0.08, 0.865 ),
annotation="spatial effect",
outfilename=outfilename
)
vn="predictions"
toplot = carstm_results_unpack( res, vn )
brks = pretty( quantile(toplot[,,"mean"], probs=c(0,0.975), na.rm=TRUE ) )
for (y in res$time_id ){
tmatch = as.character(y)
fn_root = paste(fn_root_prefix, paste0(tmatch, collapse="-"), sep="_")
outfilename = file.path( outputdir, paste(fn_root, "png", sep=".") )
carstm_map( res=res, vn=vn, tmatch=tmatch,
sppoly = sppoly,
breaks =brks,
palette="-RdYlBu",
plot_elements="",
# plot_elements=c( "compass", "scale_bar", "legend" ),
additional_features=additional_features,
legend.position=c( 0.08, 0.865 ),
annotation=y,
# title=paste(fn_root_prefix, snowcrab_filter_class, paste0(tmatch, collapse="-") )
outfilename=outfilename
)
}
# plots with 95% PI
oeffdir = file.path( outputdir, fn_root_prefix, "effects" )
if ( !file.exists(oeffdir)) dir.create( oeffdir, recursive=TRUE, showWarnings=FALSE )
(fn = file.path( oeffdir, "time.png"))
png( filename=fn, width=1024, height=1024, pointsize=12, res=196 )
carstm_plotxy( res, vn=c( "res", "random", "time" ),
type="b", xlab="Year", ylab=ylab, h=0, cex=1.25, cex.axis=1.25, cex.lab=1.25 )
dev.off()
(fn = file.path( oeffdir, "cyclic.png"))
png( filename=fn, width=1024, height=1024, pointsize=12, res=196 )
carstm_plotxy( res, vn=c( "res", "random", "cyclic" ),
type="b", col="slategray", pch=19, lty=1, lwd=2.5,
xlab="Season", ylab=ylab, cex=1.25, cex.axis=1.25, cex.lab=1.25 )
dev.off()
(fn = file.path( oeffdir, "temperature.png"))
png( filename=fn, width=1024, height=1024, pointsize=12, res=196 )
carstm_plotxy( res, vn=c( "res", "random", "inla.group(t, method = \"quantile\", n = 13)" ),
type="b", col="slategray", pch=19, lty=1, lwd=2.5 ,
xlab="Bottom temperature (degrees Celsius)", ylab=ylab, cex=1.25, cex.axis=1.25, cex.lab=1.25 )
dev.off()
(fn = file.path( oeffdir, "pca1.png"))
png( filename=fn, width=1024, height=1024, pointsize=12, res=196 )
carstm_plotxy( res, vn=c( "res", "random", "inla.group(pca1, method = \"quantile\", n = 13)" ),
type="b", col="slategray", pch=19, lty=1, lwd=2.5 ,
xlab="PCA1", ylab=ylab, cex=1.25, cex.axis=1.25, cex.lab=1.25 )
dev.off()
(fn = file.path( oeffdir, "pca2.png"))
png( filename=fn, width=1024, height=1024, pointsize=12, res=196 )
carstm_plotxy( res, vn=c( "res", "random", "inla.group(pca2, method = \"quantile\", n = 13)" ),
type="b", col="slategray", pch=19, lty=1, lwd=2.5 ,
xlab="PCA2", ylab=ylab, cex=1.25, cex.axis=1.25, cex.lab=1.25 )
dev.off()
(fn = file.path( oeffdir, "depth.png"))
png( filename=fn, width=1024, height=1024, pointsize=12, res=196 )
carstm_plotxy( res, vn=c( "res", "random", "inla.group(z, method = \"quantile\", n = 13)" ),
type="b", col="slategray", pch=19, lty=1, lwd=2.5 ,
xlab="Depth (m)", ylab=ylab, cex=1.25, cex.axis=1.25, cex.lab=1.25 )
dev.off()
# (fn = file.path( outputdir, "substrate.png"))
# png( filename=fn, width=1024, height=1024, pointsize=12, res=196 )
# carstm_plotxy( res, vn=c( "res", "random", "inla.group(substrate.grainsize, method = \"quantile\", n = 11)" ),
# type="b", col="slategray", pch=19, lty=1, lwd=2.5 ,
# xlab="Substrate grain size (mm)", ylab=ylab, cex=1.25, cex.axis=1.25, cex.lab=1.25 )
# dev.off()
res = NULL; gc()
# posterior predictive check
#vns ="number"
#vns ="meansize"
#vns ="habitat"
if ( vns=="number" ) {
MM = M[iq,]
iobs = which(MM$tag == "observations")
vn = "totno"
p = pN
outputdir = file.path( p$modeldir, p$carstm_model_label, "predicted.numerical.densities" )
}
if ( vns=="meansize" ) {
MM = M[iw,]
iobs = which(MM$tag == "observations")
vn = "meansize"
p = pW
outputdir = file.path( p$modeldir, p$carstm_model_label, "predicted.meansize" )
}
if ( vns=="habitat") {
MM = M
iobs = which(MM$tag == "observations")
vn = "pa"
p = pH
outputdir = file.path( p$modeldir, p$carstm_model_label, "predicted.presence_absence" )
}
fit = NULL; gc()
fit = carstm_model( p=p, DS="carstm_modelled_fit") #, sppoly = sppoly )
pld = data.table(
observed =MM[iobs , vn] ,
fitted = fit$summary.fitted.values[["mean"]] [iobs]
)
anno1 = paste( "Pearson correlation: ", round( cor( pld$fitted, pld$observed, use="pairwise.complete.obs" ), 3))
# cor( fitted, observed, use="pairwise.complete.obs", "spearman" )
out = ggplot(pld, aes(x = jitter(observed), y = fitted )) +
geom_abline(slope=1, intercept=0, color="darkgray", lwd=1.4 ) +
geom_point(color="slategray") +
labs(caption=anno1, color="slateblue") +
theme( plot.caption = element_text(hjust = 0, size=12 ) )# move caption to the left
fn = file.path(outputdir, "posterior_predictive_check.png" )
ggsave(filename=fn, plot=out, device="png", width=12, height = 8)
print(out)
fit = MM = NULL; gc()
}
Convolution is straightforward as it is operating upon joint posterior simulations. Add more maps/figures as required.
# wgts_max = 1.1 # kg, hard upper limit
# N_max = NULL
# # quantile( M$totno[ipositive]/M$data_offset[ipositive], probs=0.95, na.rm=TRUE )
# posterior sims
require(ggplot2)
library(ggbreak)
regions = c("cfanorth", "cfasouth", "cfa4x" )
region_label = c("N-ENS", "S-ENS", "4X")
color_map = c("#E69F00", "#56B4E9", "#CC79A7" )
sppoly = areal_units( p=pN ) # to reload
for ( vns in c("abundance", "habitat") ) {
if (vns=="abundance") {
sims = carstm_posterior_simulations( pN=pN, pW=pW, pH=pH, pa_threshold=0.05, qmax=0.95 )
sims = sims / 10^6 # units: kg ; div (10^6) -> kt ;;
# sims[ which(!is.finite(sppoly$npts)),, ] = 0
SM = aggregate_simulations(
sims=sims,
sppoly=sppoly,
fn=carstm_filenames( pN, returnvalue="filename", fn="aggregated_timeseries" ),
yrs=pN$yrs,
method="sum",
redo=TRUE
)
# units: kt
# note: using pN, even though this is biomass
outputdir = file.path( carstm_filenames( pN, returnvalue="output_directory"), "aggregated_biomass_timeseries" )
ylabel = "Biomass index (kt)"
fn_ts = "biomass_M0.png"
vn = paste("biomass", "predicted", sep=".")
outputdir2 = file.path( carstm_filenames( pN, returnvalue="output_directory"), "predicted_biomass_densities" )
}
if (vns=="habitat") {
sims = carstm_posterior_simulations( pH=pH, pa_threshold=0.05, qmax=0.95 )
SM = aggregate_simulations(
sims=sims,
sppoly=sppoly,
fn=carstm_filenames( pN, returnvalue="filename", fn="aggregated_timeseries" ),
yrs=pN$yrs,
method="mean",
redo=TRUE
)
# units: probability
# note: using pN, even though this is habitat
outputdir = file.path( carstm_filenames( pN, returnvalue="output_directory"), "aggregated_habitat_timeseries" )
ylabel = "Habitat probability"
fn_ts = "habitat_M0.png"
vn = paste("habitat", "predicted", sep=".")
outputdir2 = file.path( carstm_filenames( pN, returnvalue="output_directory"), "predicted_habitat" )
}
RES= SM$RES
if ( !file.exists(outputdir)) dir.create( outputdir, recursive=TRUE, showWarnings=FALSE )
# plot effects
( fn = file.path( outputdir, "cfa_all.png") )
png( filename=fn, width=3072, height=2304, pointsize=12, res=300 )
plot( cfaall ~ yrs, data=RES, lty="solid", lwd=4, pch=20, col="slateblue", type="b", ylab=ylabel, xlab="")
lines( cfaall_lb ~ yrs, data=RES, lty="dotted", lwd=2, col="slategray" )
lines( cfaall_ub ~ yrs, data=RES, lty="dotted", lwd=2, col="slategray" )
dev.off()
( fn = file.path( outputdir, "cfa_south.png") )
png( filename=fn, width=3072, height=2304, pointsize=12, res=300 )
plot( cfasouth ~ yrs, data=RES, lty="solid", lwd=4, pch=20, col="slateblue", type="b", ylab=ylabel, xlab="")
lines( cfasouth_lb ~ yrs, data=RES, lty="dotted", lwd=2, col="slategray" )
lines( cfasouth_ub ~ yrs, data=RES, lty="dotted", lwd=2, col="slategray" )
dev.off()
( fn = file.path( outputdir, "cfa_north.png") )
png( filename=fn, width=3072, height=2304, pointsize=12, res=300 )
plot( cfanorth ~ yrs, data=RES, lty="solid", lwd=4, pch=20, col="slateblue", type="b", ylab=ylabel, xlab="")
lines( cfanorth_lb ~ yrs, data=RES, lty="dotted", lwd=2, col="slategray" )
lines( cfanorth_ub ~ yrs, data=RES, lty="dotted", lwd=2, col="slategray" )
dev.off()
( fn = file.path( outputdir, "cfa_4x.png") )
png( filename=fn, width=3072, height=2304, pointsize=12, res=300 )
plot( cfa4x ~ yrs, data=RES, lty="solid", lwd=4, pch=20, col="slateblue", type="b", ylab=ylabel, xlab="")
lines( cfa4x_lb ~ yrs, data=RES, lty="dotted", lwd=2, col="slategray" )
lines( cfa4x_ub ~ yrs, data=RES, lty="dotted", lwd=2, col="slategray" )
dev.off()
a = cbind( "cfanorth", RES[,c("yrs", "cfanorth", "cfanorth_lb", "cfanorth_ub")] )
b = cbind( "cfasouth", RES[,c("yrs", "cfasouth", "cfasouth_lb", "cfasouth_ub")] )
c = cbind( "cfa4x", RES[,c("yrs", "cfa4x", "cfa4x_lb", "cfa4x_ub")] )
names(a) = names(b) = names(c) = c("region", "year", "mean", "lb", "ub")
tdb = rbind(a, b, c)
tdb$region = factor(tdb$region, levels=regions, labels =region_label)
tdb = tdb[(which(!is.na(tdb$region))), ]
fn = file.path( outputdir, fn_ts )
if (vns=="abundance") {
out = ggplot(tdb, aes(x=year, y=mean, fill=region, colour=region)) +
geom_line( alpha=0.9, linewidth=1.2 ) +
geom_point(aes(shape=region), size=3, alpha=0.7 ) +
geom_errorbar(aes(ymin=lb,ymax=ub), linewidth=0.8, alpha=0.8, width=0.3) +
labs(x="Year/Année", y="Biomass index (kt) / Indice de biomasse (kt)", size = rel(1.5)) +
scale_colour_manual(values=color_map) +
scale_fill_manual(values=color_map) +
scale_shape_manual(values = c(15, 17, 19)) +
theme_light( base_size = 22) +
theme( legend.position=c(0.75, 0.9), legend.title=element_blank()) +
scale_y_break(c(14, 28), scales = 1)
# scale_y_continuous( limits=c(0, 300) )
ggsave(filename=fn, plot=out, device="png", width=12, height = 8)
print(out)
# map it ..mean density
if ( !file.exists(outputdir2)) dir.create( outputdir2, recursive=TRUE, showWarnings=FALSE )
B = apply( sims, c(1,2), mean ) # sims units (kt);
B[ which(!is.finite(B)) ] = NA
# brks = pretty( log10( quantile( B[], probs=c(0.05, 0.95), na.rm=TRUE )* 10^6) )
sa = units::drop_units(sppoly$au_sa_km2)
brks = pretty( ( quantile( log(B * 10^6 / sa), probs=c(0.05, 0.95), na.rm=TRUE )) )
for (i in 1:length(pN$yrs) ) {
y = as.character( pN$yrs[i] )
# u = log10( B[,y]* 10^6 ) ## Total kt->kg: log10( kg )
u = log( B[,y]* 10^6 / sa) # ;; density log10( kg /km^2 )
sppoly[,vn] = u
outfilename = file.path( outputdir2 , paste( "biomass", y, "png", sep=".") )
carstm_map( sppoly=sppoly, vn=vn,
breaks=brks,
additional_features=additional_features,
legend.position=c( 0.1, 0.9 ),
annotation=y,
# annotation=paste( "log_10( Predicted biomass density; kg/km^2 )", y ),
colors=rev(RColorBrewer::brewer.pal(5, "RdYlBu")),
outfilename=outfilename
)
} # end year loop
}
if (vns=="habitat") {
out = ggplot(tdb, aes(x=year, y=mean, fill=region, colour=region)) +
geom_line( alpha=0.9, linewidth=1.2 ) +
geom_point(aes(shape=region), size=3, alpha=0.7 ) +
geom_errorbar(aes(ymin=lb,ymax=ub), linewidth=0.8, alpha=0.8, width=0.3) +
labs(x="Year/Année", y="Viable habitat (probability) /\nHabitat viable (probabilité)", size = rel(1.0)) +
scale_colour_manual(values=color_map) +
scale_fill_manual(values=color_map) +
scale_shape_manual(values = c(15, 17, 19)) +
theme_light( base_size = 22) +
theme( legend.position=c(0.75, 0.9), legend.title=element_blank())
# scale_y_continuous( limits=c(0, 300) )
ggsave(filename=fn, plot=out, device="png", width=12, height = 8)
print(out)
if ( !file.exists(outputdir2)) dir.create( outputdir2, recursive=TRUE, showWarnings=FALSE )
B = apply( sims, c(1,2), mean ) # sims units (kt);
B[ which(!is.finite(B)) ] = NA
# brks = pretty( log10( quantile( B[], probs=c(0.05, 0.95), na.rm=TRUE )* 10^6) )
sa = units::drop_units(sppoly$au_sa_km2)
brks = pretty( c(0, 1) )
for (i in 1:length(pN$yrs) ) {
y = as.character( pN$yrs[i] )
# u = log10( B[,y]* 10^6 ) ## Total kt->kg: log10( kg )
u = B[,y]
sppoly[,vn] = u
outfilename = file.path( outputdir2 , paste( "habitat", y, "png", sep=".") )
carstm_map( sppoly=sppoly, vn=vn,
breaks=brks,
additional_features=additional_features,
legend.position=c( 0.1, 0.9 ),
annotation=y,
# annotation=paste( "log_10( Predicted biomass density; kg/km^2 )", y ),
colors=rev(RColorBrewer::brewer.pal(5, "RdYlBu")),
outfilename=outfilename
)
} # end year loop
}
} # end vns loop
# end assimilate size and numbers
The fishery model used for stock assessment is the biomass dynamics model. It is operated through Julia. This step creates the data required.
## note the output directory .. this is used for the next script
fishery_model_label = "turing1"
carstm_results_directory = file.path( pN$modeldir, pN$carstm_model_label, "fishery_model_results", fishery_model_label ) # default
fishery_model_data_inputs( year.assessment=year.assessment, type="biomass_dynamics", for_julia=TRUE,
save_location=carstm_results_directory, fishery_model_label=fishery_model_label )
# for development, save at alternate locations
# carstm_results_directory = file.path( homedir, "projects", "dynamical_model", "snowcrab", "data" )
# fishery_model_data_inputs( year.assessment=year.assessment, type="biomass_dynamics", for_julia=TRUE, save_location=carstm_results_directory )
Rdata files are ready to load through Julia. Continue to step 04.snowcrab_fishery_model_turing.md to complete the assessment.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.