# snow crab using alt AUs
# -------------------------------------------------
# set up paramerters
yrs = 1999:2018
for ( au_type in c("lattice", "snowcrab_polygons_inla_mesh", "snowcrab_polygons_tesselation") ) {
au_type = "snowcrab_polygons_tesselation"
for ( areal_units_constraint_nmin in c( 3, 10, 15 ) ) {
# areal_units_constraint_nmin = trunc(length(yrs) / 3) # = 6
for ( dist_km in c( 1, 2, 5, 7.5, 10, 15, 20, 25, 50, 75, 100 ) ) {
# dist_km = 1
p = bio.snowcrab::snowcrab_carstm(
DS="parameters",
assessment.years=1999:2018,
modeldir = project.datadirectory("bio.snowcrab", "modelled", "testing" ), ## <--- important: alter save location for this .. default is "*/modelled"
carstm_model_label = paste( "testing", au_type, dist_km, areal_units_constraint_nmin, sep="_" ),
aegis_internal_resolution_km = 1,
boundingbox = list( xlim = c(-70.5, -56.5), ylim=c(39.5, 47.5)), # bounding box for plots using spplot
areal_units_proj4string_planar_km = projection_proj4string("utm20"), # set up default map projection
areal_units_constraint = "snowcrab",
areal_units_constraint_nmin = areal_units_constraint_nmin,
areal_units_source= au_type,
areal_units_resolution_km = dist_km,
sa_threshold_km2 = 5,
inla_num.threads = 4,
inla_blas.num.threads = 4
)
if (0) coastLayout = aegis.coastline::coastline_layout( p=p, redo=TRUE )
p = c(p, aegis.coastline::coastline_layout( p=p ) )
p$mypalette = RColorBrewer::brewer.pal(9, "YlOrRd")
# p$areal_units_proj4string_planar_km = projection_proj4string("omerc_nova_scotia") # oblique mercator, centred on Scotian Shelf rotated by 325 degrees
# -------------------------------------------------
# create polygons
if (0) {
# ensure if polys exist and create if required
for (au in c("cfanorth", "cfasouth", "cfa4x", "cfaall" )) {
plot(polygons_managementarea( species="snowcrab", au))
}
}
sppoly = areal_units( p=p )
if (0) {
# verify that params are correct
sppoly = areal_units( p=p, redo=TRUE ) # to create
plot(sppoly[,"AUID"], col="orange")
plot( sppoly@nb, coords=st_centroid(st_geometry( as(sppoly, "sf")) ), col="green", add=T )
x11(); spplot( sppoly, "au_sa_km2", main="AUID", sp.layout=p$coastLayout )
}
# -------------------------------------------------
# Part 3 -- create covariate field for bathymetry
# bathymetry -- ensure the data assimilation in bathymetry is first completed :: 01.bathymetry_data.R
pB = bathymetry_carstm( p=p, DS="parameters", variabletomodel="z" )
M = bathymetry.db( p=pB, DS="aggregated_data" , redo=TRUE )
M = bathymetry_carstm( p=pB, DS="carstm_inputs", redo=TRUE ) # will redo if not found
M = NULL; gc()
fit = carstm_model( p=pB, M='bathymetry_carstm( p=pB, DS="carstm_inputs" )', DS="redo" ) # run model and obtain predictions
if (0) {
# to use a saved instance
fit = carstm_model( p=pB, DS="carstm_modelled_fit" ) # extract currently saved model fit
summary(fit)
res = carstm_summary( p=pB ) # load summary
# maps of some of the results
vn = paste(pB$variabletomodel, "predicted", sep=".")
zplot=carstm_plot( p=pB, res=res, vn=vn )
print(zplot)
#to save map of predicted
if (0) {
fn=paste("z.predicted", year.assesment,"pdf", sep="." )
pdf(zplot, file=paste(plot.dir, fn, sep="/"))
}
vn = paste(pB$variabletomodel, "predicted_se", sep=".")
zplot=carstm_plot( p=pB, res=res, vn=vn )
print(zplot)
vn = paste(pB$variabletomodel, "random_auid_nonspatial", sep=".")
carstm_plot( p=pB, res=res, vn=vn )
vn = paste(pB$variabletomodel, "random_auid_spatial", sep=".")
carstm_plot( p=pB, res=res, vn=vn )
plot(fit, plot.prior=TRUE, plot.hyperparameters=TRUE, plot.fixed.effects=FALSE, single=TRUE )
}
# -------------------------------------------------
# Part 4 -- create covariate field for substrate
# ensure the data assimilation in substrate is first completed :: 01.substrate_data.R
pS = substrate_carstm( p=p, DS="parameters", variabletomodel="substrate.grainsize" )
M = substrate.db( p=pS, DS="aggregated_data", redo=TRUE ) # used for data matching/lookup in other aegis projects that use substrate
M = substrate_carstm( p=pS, DS="carstm_inputs", redo=TRUE ) # will redo if not found
M = NULL; gc()
fit = carstm_model( p=pS, M='substrate_carstm( p=pS, DS="carstm_inputs")', DS="redo" ) # run model and obtain predictions
if (0) {
vn = paste(pS$variabletomodel, "predicted", sep=".")
carstm_plot( p=pS, res=res, vn=vn ) # maps of some of the results
vn = paste(pS$variabletomodel, "predicted_se", sep=".")
carstm_plot( p=pS, res=res, vn=vn )
vn = paste(pS$variabletomodel, "random_auid_nonspatial", sep=".")
carstm_plot( p=pS, res=res, vn=vn )
vn = paste(pS$variabletomodel, "random_auid_spatial", sep=".")
carstm_plot( p=pS, res=res, vn=vn )
plot(fit, plot.prior=TRUE, plot.hyperparameters=TRUE, plot.fixed.effects=FALSE, single=TRUE )
}
# -------------------------------------------------
# Part 5 -- create covariate field for temperature
# ensure the data assimilation in temperature is first completed :: 01.temperature_data.R
pT = temperature_carstm( p=p, DS="parameters", variabletomodel="t" )
M = temperature.db( p=pT, DS="aggregated_data", redo=TRUE ) # used for data matching/lookup in other aegis projects that use temperature
M = temperature_carstm( p=pT, DS="carstm_inputs", redo=TRUE ) # will redo if not found
M = NULL; gc()
fit = carstm_model( p=pT, M='temperature_carstm( p=pT, DS="carstm_inputs")', DS="redo" ) # run model and obtain predictions
if (0) {
# to use a saved instance
fit = carstm_model( p=pT, DS="carstm_modelled_fit" ) # extract currently saved model fit
summary(fit)
res = carstm_summary( p=pT )
vn = paste(pT$variabletomodel, "predicted", sep=".")
carstm_plot( p=pT, res=res, vn=vn, time_match=list(year="2000", dyear="0.85" ) ) # maps of some of the results
vn = paste(pT$variabletomodel, "predicted_se", sep=".")
carstm_plot( p=pT, res=res, vn=vn, time_match=list(year="2000", dyear="0.85" ) ) # maps of some of the results
vn = paste(pT$variabletomodel, "random_auid_nonspatial", sep=".")
carstm_plot( p=pT, res=res, vn=vn, time_match=list(year="2000" ) ) # maps of some of the results
vn = paste(pT$variabletomodel, "random_auid_spatial", sep=".")
carstm_plot( p=pT, res=res, vn=vn, time_match=list(year="2000" ) ) # maps of some of the results
plot(fit, plot.prior=TRUE, plot.hyperparameters=TRUE, plot.fixed.effects=FALSE, single=TRUE )
# recent=as.character((year.assessment-6): year.assessment)
# vn = paste(pT$variabletomodel, "predicted", sep=".")
# for (x in recent){
# fn=paste(x,"t", "pdf", sep=".")
# outfile=paste(plot.dir, fn, sep="/")
# each.plot= carstm_plot( p=pT, res=res, vn=vn, time_match=list(year=x, dyear="0.85" ) )
# pdf(outfile)
# print(each.plot)
# dev.off()
# }
}
# -------------------------------------------------
# Part 6 -- create covariate field for species composition 1
# ensure that survey data is assimilated : bio.snowcrab::01snowcb_data.R, aegis.survey::01.surveys.data.R , etc.
pPC1 = speciescomposition_carstm( p=p, DS="parameters", variabletomodel="pca1" )
M = speciescomposition_carstm( p=pPC1, DS="carstm_inputs", redo=TRUE ) # will redo if not found
M = NULL; gc()
fit = carstm_model( p=pPC1, M='speciescomposition_carstm( p=p, DS="carstm_inputs" )', DS="redo" ) # run model and obtain predictions
if (0) {
# to use a saved instance
fit = carstm_model( p=pPC1, DS="carstm_modelled_fit" ) # extract currently saved model fit
summary(fit)
res = carstm_summary( p=pPC1 )
plot( fit, plot.prior=TRUE, plot.hyperparameters=TRUE, plot.fixed.effects=FALSE )
vn = paste(pPC1$variabletomodel, "predicted", sep=".")
carstm_plot( p=pPC1, res=res, vn=vn, time_match=list(year="2017" ), dyear="0.85" ) # maps of some of the results
vn = paste(pPC1$variabletomodel, "predicted_se", sep=".")
carstm_plot( p=pPC1, res=res, vn=vn, time_match=list(year="2017" ), dyear="0.85" ) # maps of some of the results
vn = paste(pPC1$variabletomodel, "random_auid_nonspatial", sep=".")
carstm_plot( p=pPC1, res=res, vn=vn, time_match=list(year="2017" ), dyear="0.85" ) # maps of some of the results , dyear="0.85"
vn = paste(pPC1$variabletomodel, "random_auid_spatial", sep=".")
carstm_plot( p=pPC1, res=res, vn=vn, time_match=list(year="2017" ), dyear="0.85" ) # maps of some of the results , dyear="0.85"
}
# -------------------------------------------------
# Part 7 -- create covariate field for species composition 2
# ensure that survey data is assimilated : bio.snowcrab::01snowcb_data.R, aegis.survey::01.surveys.data.R ,
pPC2 = speciescomposition_carstm( p=p, DS="parameters", variabletomodel="pca2")
M = speciescomposition_carstm( p=pPC2, DS="carstm_inputs", redo=TRUE ) # will redo if not found
M = NULL; gc()
fit = carstm_model( p=pPC2, M='speciescomposition_carstm( p=p, DS="carstm_inputs" )', DS="redo" ) # run model and obtain predictions
if (0) {
# to use a saved instance
fit = carstm_model( p=pPC2, DS="carstm_modelled_fit" ) # extract currently saved model fit
summary(fit)
res = carstm_summary( p=pPC2 )
plot( fit, plot.prior=TRUE, plot.hyperparameters=TRUE, plot.fixed.effects=FALSE )
vn = paste(pPC2$variabletomodel, "predicted", sep=".")
carstm_plot( p=pPC2, res=res, vn=vn, time_match=list(year="2017" ), dyear="0.85" ) # maps of some of the results
vn = paste(pPC2$variabletomodel, "predicted_se", sep=".")
carstm_plot( p=pPC2, res=res, vn=vn, time_match=list(year="2017" ), dyear="0.85" ) # maps of some of the results
vn = paste(pPC2$variabletomodel, "random_auid_nonspatial", sep=".")
carstm_plot( p=pPC2, res=res, vn=vn, time_match=list(year="2017" ), dyear="0.85" ) # maps of some of the results , dyear="0.85"
vn = paste(pPC2$variabletomodel, "random_auid_spatial", sep=".")
carstm_plot( p=pPC2, res=res, vn=vn, time_match=list(year="2017" ), dyear="0.85" ) # maps of some of the results , dyear="0.85"
}
# finished covariates ... move onto abundance index estimation
# -------------------------------------------------
# Part 8 -- Snow crab anbundance -- main mode used for production
M = snowcrab_carstm( p=p, DS="carstm_inputs", redo=TRUE ) # will redo if not found
M = NULL; gc()
fit = carstm_model( p=p, M='snowcrab_carstm( p=p, DS="carstm_inputs" )' ) # 151 configs and long optim .. 19 hrs
fit = carstm_model( p=p, DS="carstm_modelled_fit" ) # extract currently saved model fit
summary(fit)
res = carstm_summary( p=p )
RES = snowcrab_carstm(p=p, DS="carstm_output_compute" )
if (0) {
vn = paste(p$variabletomodel, "predicted", sep=".")
carstm_plot( p=p, res=res, vn=vn, time_match=list(year="2000" ) ) # maps of some of the results
plot(fit)
plot(fit, plot.prior=TRUE, plot.hyperparameters=TRUE, plot.fixed.effects=FALSE, single=TRUE )
s = summary(fit)
s$dic$dic
s$dic$p.eff
# maps of some of the results
vn = paste(p$variabletomodel, "random_sample_iid", sep=".")
if (exists(vn, res)) carstm_plot( p=p, res=res, vn=vn, time_match=list(year="1950", dyear="0.05") )
vn = paste(p$variabletomodel, "random_auid_nonspatial", sep=".")
if (exists(vn, res)) {
res_dim = dim( res[[vn]] )
if (res_dim == 1 ) time_match = NULL
if (res_dim == 2 ) time_match = list(year="2000")
if (res_dim == 3 ) time_match = list(year="2000", dyear="0.85" )
carstm_plot( p=p, res=res, vn=vn, time_match=time_match )
}
vn = paste(p$variabletomodel, "random_auid_spatial", sep=".")
if (exists(vn, res)) {
res_dim = dim( res[[vn]] )
if (res_dim == 1 ) time_match = NULL
if (res_dim == 2 ) time_match = list(year="2000")
if (res_dim == 3 ) time_match = list(year="2000", dyear="0.85" )
carstm_plot( p=p, res=res, vn=vn, time_match=time_match )
}
RES = snowcrab_carstm(p=p, DS="carstm_output_timeseries" )
bio = snowcrab_carstm(p=p, DS="carstm_output_spacetime_biomass" )
num = snowcrab_carstm(p=p, DS="carstm_output_spacetime_number" )
plot( cfaall ~ yrs, data=RES, lty="solid", lwd=4, pch=20, col="slateblue", type="b", ylab="Biomass (kt)", 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" )
plot( cfasouth ~ yrs, data=RES, lty="solid", lwd=4, pch=20, col="slateblue", type="b", ylab="Biomass (kt)", 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" )
plot( cfanorth ~ yrs, data=RES, lty="solid", lwd=4, pch=20, col="slateblue", type="b", ylab="Biomass (kt)", 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" )
plot( cfa4x ~ yrs, data=RES, lty="solid", lwd=4, pch=20, col="slateblue", type="b", ylab="Biomass (kt)", 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" )
p$boundingbox = list( xlim=p$corners$lon, ylim=p$corners$lat) # bounding box for plots using spplot
p$coastLayout = aegis.coastline::coastline_layout(p=p)
p$mypalette=RColorBrewer::brewer.pal(9, "YlOrRd")
sppoly = areal_units( p=p ) # to reload
# map it ..mean density
vn = "pred"
sppoly@data[,vn] = bio[,"2017"]
brks = interval_break(X= sppoly[[vn]], n=length(p$mypalette), style="quantile")
spplot( sppoly, vn, col.regions=p$mypalette, main=vn, at=brks, sp.layout=p$coastLayout, col="transparent" )
plot( fit, plot.prior=TRUE, plot.hyperparameters=TRUE, plot.fixed.effects=FALSE )
plot( fit$marginals.hyperpar$"Phi for auid", type="l") # posterior distribution of phi nonspatial dominates
plot( fit$marginals.hyperpar$"Precision for auid", type="l")
plot( fit$marginals.hyperpar$"Precision for setno", type="l")
}
}}}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.