inst/scripts/99_example_snowcrab.R

In jae0/ecofishmod: Ecosystem-based fishery stock assessment models using conditional autoregressive models to represent spatial and temporal variabilty.

# 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")
  }


}}}