Update docs and examples (#69)

* credit Cecilia O'Leary for convert_shapefile(.)

* update docs

* update docs

* update docs

* Adding CalCOFI extrapolation-grid option...

... using coordinates defining shapefile from Andrew Thompson

* add CalCOFI+IMECOCAL summer extrapolation-grid option

* improve docs

* improved docs

* add juvenile rockfish shapefiles and improve docs

* Update plot_anisotropy.R

* update plot_factors(.) ...

... to plot standard errors in factor loadings when `SD` is provided

* improve plot_factors to extract SE for response-maps

* adding main Hawaiian Islands extrapolation-grid...

... thanks to Ben Richards

* fixing plot_results for new covariate interface

thanks to E. Yasumiishi for pointing out the issue

* small update to docs

* update version number

* adding plots for spatialy-varying density response

* Update simulate_data.R

to allow multiple DLLs to be loaded when using simulate_data.  h/t @dshanisko and @kellijohnson-NOAA for flagging the issue and discussing a (temporary) fix.

* streamline interface...

... when Hessian calculation fails, to allow results to still be outputted.

* small doc improvements

* adding interface for one-step-ahead residuals...

... using TMBhelper::oneStepPredict_deltaModel.  Only available using CPP version >= 13.0.0

* fix rotate_factors(.)...

... to improve denominator in error-testing for use when zero-centering loadings matrices for use in EOF

* small updates

* h/t Arnaud in VAST #55 and #56

* Update make_settings documentation

* Add purpose to make_kmeans to save extrapolation kmeans file to avoid recalculating it

* Clean up file paths a bit

* Fix tiny typo bug

* Rename new argument and files per JTT request

* Streamline console output from make_kmeans

* fix file-path specification

* adding informative warning..

... about change in make_kmeans

* Update in the PESC example

Update in the PESC example for which we now provide two datasets instead of just one: a predator biomass catch rate dataset and a stomach content dataset.

* adding AFSC shapefiles

* fix bug in PESC update to `load_example`

* update docs

* fix small bug in improvements to `make_kmeans` ...

... where previous change made `make_extrapolation_info` made kmeans object in working-directory that sometimes doesn't have write access; now uses user-specified `working_dir`

* Fix plotting bug in DHARMa quantile residuals...

... wherein the color-scale was accidentally inverted (i.e., a residual >0.5 previously indicated data below the median of the predictive distribution, in contrast to standard definition of quantile residuals)

* additional plotting option in plot_timeseries

* integrating package effects to visualize covariate responses

* Response to #64

* change plot defaults to use projargs specified during fitting

* update predict.fit_model ...

... which seems to be working as intended and interfacing correctly with package pdp for partial dependence plots

* add Cold-Pool-Extent to pollock example ...

... for use in new Wiki demo

* fix bug in recent improvements to predict.fit_model

* small updates to reference docs

* a bunch of small doc improvements

* improve docs

* Add plot_factors within plot_results and improve docs

* adding purpose = "EOF3" as default for EOF

* fix bug in plot_biomass_index when Ftarg is known without error...

... i.e., in data-poor production models.  Also improve docs.

* improve plot_variable to allow mapping points (instead of rasters)

* fix bug which was causing error in convert_shapefile using R 4.0.2

* switch to plotting maps using viridis colorscheme

* Add better error checking in get_latest_version

* Small bugs

* Add path argument to get_latest_version to run local tests

* Change what happens if user path fails

* Fix small bug

* switch rotate_factors to using SVD...

... which avoids problems in Cholesky decomposition when using a rank-deficient covariance

* update interface for plot_range_edge

* update to facilitate seasonal model

* adding seasonal example data

* small doc fixes

* update DESCRIPTION

Co-authored-by: Cole-Monnahan-NOAA <[email protected]>
Co-authored-by: Arnaud Grüss <[email protected]>

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: FishStatsUtils
 Type: Package
 Title: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox
-Version: 2.8.0
-Date: 2020-09-22
+Version: 2.9.0
+Date: 2021-03-12
 Authors@R: 
     c(person(given = "James",
              family = "Thorson",
@@ -42,9 +42,10 @@ Imports:
     rnaturalearthdata,
     formatR,
     splancs,
-    DHARMa
+    DHARMa,
+    viridisLite
 Depends:
-    R (>= 3.1.0)
+    R (>= 3.5.0)
 Suggests:
     testthat
 Remotes: 

NAMESPACE

@@ -16,6 +16,7 @@ export(PlotLoadings)
 export(PlotMap_Fn)
 export(Prepare_AI_Extrapolation_Data_Fn)
 export(Prepare_BC_Coast_Extrapolation_Data_Fn)
+export(Prepare_BFISH_MHI_Extrapolation_Data_Fn)
 export(Prepare_BSslope_Extrapolation_Data_Fn)
 export(Prepare_EBS_Extrapolation_Data_Fn)
 export(Prepare_GOA_Extrapolation_Data_Fn)

R/Prepare_XXXX_Extrapolation_Data_Fn.R

@@ -485,6 +485,7 @@ function( strata.limits, observations_LL, projargs=NA, grid_dim_km=c(2,2), maxim
     Data_Extrap = expand.grid( "E_km"=E_seq, "N_km"=N_seq, "Area_km2"=prod(grid_dim_km) )
 
     # Add LL
+    # NOTE:  inputs to project_coordinates are reversed due to projecting from projargs back to origargs
     #TmpUTM = rename_columns( Data_Extrap[,c("E_km","N_km")], newname=c("X","Y"))
     #attr(TmpUTM, "projection") = "UTM"
     #attr(TmpUTM, "zone") = attr(TmpUTM,"zone")
@@ -743,3 +744,42 @@ function( strata.limits=NULL, projargs=NA, zone=NA, flip_around_dateline=FALSE,
     "flip_around_dateline"=flip_around_dateline, "Area_km2_x"=Area_km2_x)
   return( Return )
 }
+
+#' @export
+Prepare_BFISH_MHI_Extrapolation_Data_Fn <-
+function( strata.limits=NULL, projargs=NA, zone=NA, flip_around_dateline=FALSE, ... ){
+  # Infer strata
+  if( is.null(strata.limits)){
+    strata.limits = data.frame('STRATA'="All_areas")
+  }
+  message("Using strata ", strata.limits)
+
+  # Read extrapolation data
+  utils::data( BFISH_MHI, package="FishStatsUtils" )
+  Data_Extrap <- BFISH_MHI
+
+  # Survey areas
+  Area_km2_x = Data_Extrap[,'Area_km2']
+
+  # Augment with strata for each extrapolation cell
+  Tmp = cbind("BEST_DEPTH_M"=Data_Extrap[,'Depth_MEAN_m'], "BEST_LAT_DD"=Data_Extrap[,'lat_deg'], "BEST_LON_DD"=Data_Extrap[,'lon_deg'])
+  a_el = as.data.frame(matrix(NA, nrow=nrow(Data_Extrap), ncol=length(strata.limits), dimnames=list(NULL,names(strata.limits))))
+  for(l in 1:ncol(a_el)){
+    a_el[,l] = apply(Tmp, MARGIN=1, FUN=match_strata_fn, strata_dataframe=strata.limits[l,,drop=FALSE])
+    a_el[,l] = ifelse( is.na(a_el[,l]), 0, Area_km2_x)
+  }
+
+  # Convert extrapolation-data to an Eastings-Northings coordinate system
+  tmpUTM = project_coordinates( X=Data_Extrap[,'lon_deg'], Y=Data_Extrap[,'lat_deg'], projargs=projargs, zone=zone, flip_around_dateline=flip_around_dateline)                                                         #$
+
+  # Extra junk
+  Data_Extrap = cbind( Data_Extrap, 'Include'=1)
+  Data_Extrap[,c('E_km','N_km')] = tmpUTM[,c('X','Y')]
+  colnames(Data_Extrap)[1:2] = c("Lon","Lat")
+
+  # Return
+  Return = list( "a_el"=a_el, "Data_Extrap"=Data_Extrap, "zone"=attr(tmpUTM,"zone"), "projargs"=attr(tmpUTM,"projargs"),
+    "flip_around_dateline"=flip_around_dateline, "Area_km2_x"=Area_km2_x)
+  return( Return )
+}
+

R/calculate_proportion.R

@@ -15,10 +15,26 @@
 #'   \item{Neff_tl}{Effective sample size (median across categories) for each time t and stratum l}
 #' }
 #'
+#' @references For details regarding multivariate index standardization and expansion see \url{https://cdnsciencepub.com/doi/full/10.1139/cjfas-2018-0015}
 #' @export
-calculate_proportion = function( TmbData, Index, Expansion_cz=NULL, Year_Set=NULL, Years2Include=NULL, strata_names=NULL, category_names=NULL,
-  plot_legend=ifelse(TmbData$n_l>1,TRUE,FALSE), DirName=paste0(getwd(),"/"), PlotName="Proportion.png", PlotName2="Average.png",
-  interval_width=1, width=6, height=6, xlab="Category", ylab="Proportion", ... ){
+calculate_proportion <-
+function( TmbData,
+          Index,
+          Expansion_cz = NULL,
+          year_labels = NULL,
+          years_to_plot = NULL,
+          strata_names = NULL,
+          category_names = NULL,
+          plot_legend = ifelse(TmbData$n_l>1,TRUE,FALSE),
+          DirName = paste0(getwd(),"/"),
+          PlotName = "Proportion.png",
+          PlotName2 = "Average.png",
+          interval_width = 1,
+          width = 6,
+          height = 6,
+          xlab = "Category",
+          ylab = "Proportion",
+          ... ){
 
   # Warnings and errors
   if( !all(TmbData[['FieldConfig']] %in% c(-3,-2,-1)) ){
@@ -58,7 +74,7 @@ calculate_proportion = function( TmbData, Index, Expansion_cz=NULL, Year_Set=NUL
 
   # Plot
   if( !is.na(PlotName) ){
-    plot_index( Index_ctl=Prop_ctl, sd_Index_ctl=sqrt(var_Prop_ctl), Year_Set=Year_Set, Years2Include=Years2Include,
+    plot_index( Index_ctl=Prop_ctl, sd_Index_ctl=sqrt(var_Prop_ctl), year_labels=year_labels, years_to_plot=years_to_plot,
       strata_names=strata_names, category_names=category_names, plot_legend=plot_legend,
       DirName=DirName, PlotName=PlotName, interval_width=interval_width, width=width, height=height,
       xlab=xlab, ylab=ylab, scale="uniform", ... )
@@ -73,7 +89,7 @@ calculate_proportion = function( TmbData, Index, Expansion_cz=NULL, Year_Set=NUL
 
   # Plot
   if( !is.na(PlotName2) ){
-    plot_index( Index_ctl=1%o%Mean_tl, sd_Index_ctl=1%o%sd_Mean_tl, Year_Set=Year_Set, Years2Include=Years2Include,
+    plot_index( Index_ctl=1%o%Mean_tl, sd_Index_ctl=1%o%sd_Mean_tl, year_labels=year_labels, years_to_plot=years_to_plot,
       strata_names=strata_names, category_names=category_names, plot_legend=plot_legend,
       DirName=DirName, PlotName=PlotName2, interval_width=interval_width, width=width, height=height,
       xlab=xlab, ylab="Category", scale="uniform", Yrange=c(NA,NA), ... )     # , Yrange=c(1,dim(var_Prop_ctl)[1])

R/convert_shapefile.R

@@ -18,18 +18,26 @@
 #' \dontrun{
 #'  convert_shapefile( file_path="C:/Users/James.Thorson/Desktop/Work files/AFSC/2020-03 -- Add ICES grids/IBTS grids/BITS/Shapefile.shp", make_plots=TRUE )
 #' }
-
+#'
+#' @author Cecilia O'Leary, James Thorson
 #' @export
-convert_shapefile = function( file_path, projargs=NULL, grid_dim_km=c(2,2), projargs_for_shapefile=NULL,
-  make_plots=FALSE, quiet=TRUE, area_tolerance=0.05, ... ){
-
-  shapefile_orig = rgdal::readOGR( file_path, verbose=FALSE, p4s=projargs_for_shapefile )
+convert_shapefile = function( file_path,
+    projargs = NULL,
+    grid_dim_km = c(2,2),
+    projargs_for_shapefile = NULL,
+    make_plots = FALSE,
+    quiet = TRUE,
+    area_tolerance = 0.05,
+    ... ){
+
+  shapefile_input = rgdal::readOGR( file_path, verbose=FALSE, p4s=projargs_for_shapefile )
+  message("Reading shapefile with projargs: ", print(shapefile_input@proj4string))
   # raster::shapefile(.) has simplified read-write interface for future reference
-  if( !(shapefile_orig@class %in% c("SpatialPolygonsDataFrame","SpatialPolygons")) ){
+  if( !(shapefile_input@class %in% c("SpatialPolygonsDataFrame","SpatialPolygons")) ){
     warning( "object at `file_path` doesn't appear to be a shapefile")
   }
   proj_orig = "+proj=longlat +ellps=WGS84 +no_defs"
-  shapefile_orig = sp::spTransform(shapefile_orig, CRSobj=sp::CRS(proj_orig) )
+  shapefile_orig = sp::spTransform(shapefile_input, CRSobj=sp::CRS(proj_orig) )
   #shapefile_orig@proj4string = sp::CRS(proj_orig)
 
   # Infer projargs if missing, and project
@@ -89,6 +97,7 @@ convert_shapefile = function( file_path, projargs=NULL, grid_dim_km=c(2,2), proj
   }
 
   # Compare areas
+  area_shapefile_input = sum( raster::area(shapefile_input) )
   area_shapefile_orig = sum( raster::area(shapefile_orig) ) / 1000^2  # Convert to square-kiometers
   area_shapefile_proj = sum( raster::area(shapefile_proj) )
   area_grid_proj = sum( extrapolation_grid[,'Area_km2'] )