inst/scripts/03.biomass_index_carstm.r
In jae0/snowcrab: Snow crab assessment suite for aegis* and associated projects.
---
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
---
<!-- Preamble
This is a Markdown document designed to be viewed by a web browser.
We do not convert to a Quarto/Rmarkdown report as the step involved are complex and may require tweaking.
-->
## Purpose
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](https://github.com/jae0/carstm) front-end to [INLA](https://www.r-inla.org/), used to perform "non-separable" spatial Conditional autocorrelation (CAR) and temporal (AR1) models.
- Poisson model of positive valued numbers offset by swept area in space and time
- Mean body size (mass; kg) in space and time
- Presence-absence of organisms in space and time, where presence or absence is defined on a quantile threshold of numerical density presence_absence
The convolution of all three after posterior simulation is also known as a Hurdle or Delta model.
## Part 1 -- construct basic parameter list defining the main characteristics of the study
```r
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)
```
## Part 2 -- spatiotemporal statistical model
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.
```r
# 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()
}
```
## Part 3: assimilation of models
Convolution is straightforward as it is operating upon joint posterior simulations. Add more maps/figures as required.
```r
# 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
```
## Part 4: prep data for integration into the fishery model (discrete version of the logistic model)
The fishery model used for stock assessment is the biomass dynamics model. It is operated through Julia. This step creates the data required.
```r
## 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](04.snowcrab_fishery_model_turing.md) to complete the assessment.
# end
jae0/snowcrab documentation built on Feb. 27, 2024, 2:42 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
---
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
---
<!-- Preamble
This is a Markdown document designed to be viewed by a web browser.
We do not convert to a Quarto/Rmarkdown report as the step involved are complex and may require tweaking.
-->
## Purpose
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](https://github.com/jae0/carstm) front-end to [INLA](https://www.r-inla.org/), used to perform "non-separable" spatial Conditional autocorrelation (CAR) and temporal (AR1) models.
- Poisson model of positive valued numbers offset by swept area in space and time
- Mean body size (mass; kg) in space and time
- Presence-absence of organisms in space and time, where presence or absence is defined on a quantile threshold of numerical density presence_absence
The convolution of all three after posterior simulation is also known as a Hurdle or Delta model.
## Part 1 -- construct basic parameter list defining the main characteristics of the study
```r
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)
```
## Part 2 -- spatiotemporal statistical model
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.
```r
# 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()
}
```
## Part 3: assimilation of models
Convolution is straightforward as it is operating upon joint posterior simulations. Add more maps/figures as required.
```r
# 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
```
## Part 4: prep data for integration into the fishery model (discrete version of the logistic model)
The fishery model used for stock assessment is the biomass dynamics model. It is operated through Julia. This step creates the data required.
```r
## 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](04.snowcrab_fishery_model_turing.md) to complete the assessment.
# end
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.