Improve plots and interface

Improving plots and interface by:

Using sp package to plot variables
Simplifying passing optional arguments to fit_model
Adding formula interface for density covariates.

DESCRIPTION

@@ -1,12 +1,12 @@
 Package: FishStatsUtils
 Type: Package
 Title: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox
-Version: 2.2.0
-Date: 2019-07-25
+Version: 2.3.0
+Date: 2019-09-18
 Authors@R: person("James","Thorson", email="[email protected]", role=c("aut","cre"))
 Maintainer: James Thorson <[email protected]>
 Description: FishStatsUtils contains utilities (shared code and data) used by multiple 
-  packages (VAST, SpatialDeltaGLMM, MIST, SPatial_FA, SpatialDFA, surplus_production) that are designed
+  packages (VAST, SpatialDeltaGLMM, MIST, Spatial_FA, SpatialDFA, surplus_production, EOFR) that are designed
   for spatio-temporal analysis of ecological data.  
 Imports:
     graphics,
@@ -24,17 +24,18 @@ Imports:
     devtools,
     mixtools,
     sp,
+    plotKML,
+    plotrix,
     TMB,
     MatrixModels,
     rgdal,
     ThorsonUtilities,
     abind,
     corpcor,
     pander,
+    rnaturalearth,
     formatR
 Depends:
-    maps,
-    mapdata,
     R (>= 3.1.0)
 Suggests:
     testthat

NAMESPACE

@@ -43,6 +43,7 @@ export(fit_model)
 export(format_covariates)
 export(get_latest_version)
 export(load_example)
+export(make_covariates)
 export(make_extrapolation_info)
 export(make_map_info)
 export(make_settings)
@@ -62,9 +63,11 @@ export(plot_loadings)
 export(plot_maps)
 export(plot_overdispersion)
 export(plot_quantile_diagnostic)
+export(plot_range_edge)
 export(plot_range_index)
 export(plot_residuals)
 export(plot_results)
 export(plot_timeseries)
+export(plot_variable)
 export(rotate_factors)
 export(summarize_covariance)

R/PlotMap_Fn.R

@@ -23,6 +23,9 @@ function(MappingDetails, Mat, PlotDF, MapSizeRatio=c('Width(in)'=4,'Height(in)'=
          Legend=list("use"=FALSE, "x"=c(10,30), "y"=c(10,30)), mfrow=c(1,1), plot_legend_fig=TRUE, land_color="grey", ignore.na=FALSE,
          map_style="rescale", ...){
 
+  # Warning
+  warning( "`PlotMap_Fn` is soft-deprecated, please use `plot_variable` instead for many improvements")
+
   # Check for problems
   if( length(Year_Set) != ncol(Mat) ){
     warning( "Year_Set and `ncol(Mat)` don't match: Changing Year_Set'")

R/Prepare_AI_Extrapolation_Data_Fn.R

@@ -1,6 +1,6 @@
 #' @export
 Prepare_AI_Extrapolation_Data_Fn <-
-function( strata.limits=NULL, zone=NA, ... ){
+function( strata.limits=NULL, zone=NA, flip_around_dateline=TRUE, ... ){
   # Infer strata
   if( is.null(strata.limits)){
     strata.limits = data.frame('STRATA'="All_areas")
@@ -24,13 +24,13 @@ function( strata.limits=NULL, zone=NA, ... ){
   }
              #
   # Convert extrapolation-data to an Eastings-Northings coordinate system
-  tmpUTM = Convert_LL_to_UTM_Fn( Lon=Data_Extrap[,'Lon'], Lat=Data_Extrap[,'Lat'], zone=zone, flip_around_dateline=TRUE)                                                         #$
+  tmpUTM = Convert_LL_to_UTM_Fn( Lon=Data_Extrap[,'Lon'], Lat=Data_Extrap[,'Lat'], 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')]
 
   # Return
-  Return = list( "a_el"=a_el, "Data_Extrap"=Data_Extrap, "zone"=attr(tmpUTM,"zone"), "flip_around_dateline"=TRUE, "Area_km2_x"=Area_km2_x)
+  Return = list( "a_el"=a_el, "Data_Extrap"=Data_Extrap, "zone"=attr(tmpUTM,"zone"), "flip_around_dateline"=flip_around_dateline, "Area_km2_x"=Area_km2_x)
   return( Return )
 }

R/Prepare_BC_Coast_Extrapolation_Data_Fn.r

@@ -1,6 +1,6 @@
 #' @export
 Prepare_BC_Coast_Extrapolation_Data_Fn <-
-function( strata.limits=NULL, strata_to_use=c('SOG','WCVI','QCS','HS','WCHG'), zone=NA, ... ){
+function( strata.limits=NULL, strata_to_use=c('SOG','WCVI','QCS','HS','WCHG'), zone=NA, flip_around_dateline=FALSE, ... ){
   # Infer strata
   if( is.null(strata.limits)){
     strata.limits = data.frame('STRATA'="All_areas")
@@ -24,13 +24,13 @@ function( strata.limits=NULL, strata_to_use=c('SOG','WCVI','QCS','HS','WCHG'), z
   }
 
   # Convert extrapolation-data to an Eastings-Northings coordinate system
-  tmpUTM = Convert_LL_to_UTM_Fn( Lon=Data_Extrap[,'Lon'], Lat=Data_Extrap[,'Lat'], zone=zone)
+  tmpUTM = Convert_LL_to_UTM_Fn( Lon=Data_Extrap[,'Lon'], Lat=Data_Extrap[,'Lat'], 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')]
 
   # Return
-  Return = list( "a_el"=a_el, "Data_Extrap"=Data_Extrap, "zone"=attr(tmpUTM,"zone"), "flip_around_dateline"=FALSE, "Area_km2_x"=Area_km2_x)
+  Return = list( "a_el"=a_el, "Data_Extrap"=Data_Extrap, "zone"=attr(tmpUTM,"zone"), "flip_around_dateline"=flip_around_dateline, "Area_km2_x"=Area_km2_x)
   return( Return )
 }

R/Prepare_NZ_Extrapolation_Data_Fn.R

@@ -1,6 +1,6 @@
 #' @export
 Prepare_NZ_Extrapolation_Data_Fn <-
-function( strata.limits=NULL, zone=NA, survey="Chatham_rise", ... ){
+function( strata.limits=NULL, zone=NA, survey="Chatham_rise", flip_around_dateline=FALSE, ... ){
   # Infer strata
   if( is.null(strata.limits)){
     strata.limits = data.frame('STRATA'="All_areas")
@@ -30,7 +30,7 @@ function( strata.limits=NULL, zone=NA, survey="Chatham_rise", ... ){
 
   # Convert extrapolation-data to an Eastings-Northings coordinate system
   if( is.numeric(zone) ){
-    tmpUTM = Convert_LL_to_UTM_Fn( Lon=Data_Extrap[,'Lon'], Lat=Data_Extrap[,'Lat'], zone=zone, flip_around_dateline=FALSE)                                                         #$
+    tmpUTM = Convert_LL_to_UTM_Fn( Lon=Data_Extrap[,'Lon'], Lat=Data_Extrap[,'Lat'], zone=zone, flip_around_dateline=flip_around_dateline)                                                         #$
     colnames(tmpUTM) = c('E_km','N_km')
   }else{
     tmpUTM = Convert_LL_to_EastNorth_Fn( Lon=Data_Extrap[,'Lon'], Lat=Data_Extrap[,'Lat'], crs=zone )
@@ -41,6 +41,6 @@ function( strata.limits=NULL, zone=NA, survey="Chatham_rise", ... ){
   Data_Extrap[,c('E_km','N_km')] = tmpUTM[,c('E_km','N_km')]
 
   # Return
-  Return = list( "a_el"=a_el, "Data_Extrap"=Data_Extrap, "zone"=attr(tmpUTM,"zone"), "flip_around_dateline"=FALSE, "Area_km2_x"=Area_km2_x)
+  Return = list( "a_el"=a_el, "Data_Extrap"=Data_Extrap, "zone"=attr(tmpUTM,"zone"), "flip_around_dateline"=flip_around_dateline, "Area_km2_x"=Area_km2_x)
   return( Return )
 }

R/calculate_proportion.R

@@ -40,6 +40,7 @@ calculate_proportion = function( TmbData, Index, Year_Set=NULL, Years2Include=NU
     #var_Prop_ctl[cI,tI,lI] = Index_ctl[cI,tI,lI]^2/Index_tl[tI,lI]^2 * (SE_Index_ctl[cI,tI,lI]^2/Index_ctl[cI,tI,lI]^2  + SE_Index_tl[tI,lI]^2/Index_tl[tI,lI]^2 )
     # Slightly extended version
     var_Prop_ctl[cI,tI,lI] = Index_ctl[cI,tI,lI]^2/Index_tl[tI,lI]^2 * (SE_Index_ctl[cI,tI,lI]^2/Index_ctl[cI,tI,lI]^2 - 2*SE_Index_ctl[cI,tI,lI]^2/(Index_ctl[cI,tI,lI]*Index_tl[tI,lI]) + SE_Index_tl[tI,lI]^2/Index_tl[tI,lI]^2 )
+    var_Prop_ctl[cI,tI,lI] = ifelse( Index_ctl[cI,tI,lI]==0, 0, var_Prop_ctl[cI,tI,lI] )  # If dividing by zero, replace with 0
     # Covert to effective sample size
     Neff_ctl[cI,tI,lI] = Prop_ctl[cI,tI,lI] * (1-Prop_ctl[cI,tI,lI]) / var_Prop_ctl[cI,tI,lI]
   }}}

R/combine_extrapolation_info.R

@@ -8,7 +8,7 @@
 #' @return Identical output from \code{FishStatsUtils::make_extrapolation_info}, but combined from each input
 
 #' @export
-combine_extrapolation_info = function( ... ){
+combine_extrapolation_info = function( ..., create_strata_per_region=FALSE ){
 
   input_list = list( ... )
 
@@ -24,17 +24,47 @@ combine_extrapolation_info = function( ... ){
   }
 
   # Combine stuff
-  a_el = Data_Extrap = Area_km2_x = NULL
+  Data_Extrap = Area_km2_x = NULL
+  a1_el = matrix(0, nrow=0, ncol=1)
+  a2_el = matrix(0, nrow=0, ncol=0)
   #assign( x="input_list", value=input_list, envir = .GlobalEnv )
 
   for( lI in 1:length(input_list) ){
-    Tmp = input_list[[lI]]$Data_Extrap
-    colnames(Tmp) = ifelse( colnames(Tmp)=="Area_in_survey_km2", "Area_km2", colnames(Tmp) )
-    Data_Extrap = rbind( Data_Extrap, Tmp[,c('E_km','N_km','Lon','Lat','Include','Area_km2')] )
-    a_el = rbind( a_el, input_list[[lI]]$a_el )
+    #Tmp = input_list[[lI]]$Data_Extrap
+    #colnames(Tmp) = ifelse( colnames(Tmp)=="Area_in_survey_km2", "Area_km2", colnames(Tmp) )
+    #Data_Extrap = rbind( Data_Extrap, Tmp[,c('E_km','N_km','Lon','Lat','Include','Area_km2')] )
+
+    # Warnings
+    if( ncol(input_list[[lI]]$a_el)>1 ){
+      if( !(create_strata_per_region==TRUE & lI==1) ){
+        stop("`combine_extrapolation_info` isn't designed to combine regions with multiple identified strata, except when `create_strata_per_region=TRUE`")
+      }
+    }
+
+    # Combine Data_Extrap
+    Data_Extrap = rbind( Data_Extrap, input_list[[lI]]$Data_Extrap[,c('E_km','N_km','Lon','Lat','Include')] )
+
+    # Combine area vector
     Area_km2_x = c( Area_km2_x, input_list[[lI]]$Area_km2_x )
+
+    # Combine strata definitions
+    a1_el = rbind( as.matrix(a1_el), as.matrix(input_list[[lI]]$a_el[,1,drop=FALSE]) )
+
+    # Make one stratum per region
+    a2_el = rbind(
+      as.matrix(cbind(a2_el, matrix(0,nrow=nrow(a2_el),ncol=ncol(input_list[[lI]]$a_el)))),
+      as.matrix(cbind(matrix(0,nrow=nrow(input_list[[lI]]$a_el),ncol=ncol(a2_el)), input_list[[lI]]$a_el))
+    )
+  }
+
+  # Only pass back the
+  if( create_strata_per_region==TRUE ){
+    a_el = a2_el
+  }else{
+    a_el = a1_el
   }
 
+  # Return
   Return = list( "a_el"=a_el, "Data_Extrap"=Data_Extrap, "zone"=Zone[1], "flip_around_dateline"=Flip[1], "Area_km2_x"=Area_km2_x)
 }
 

R/deprecated.R

@@ -732,3 +732,50 @@ Plot_States_in_UTM_Fn = function( MappingDetails, Rotate=0, fillcol=NA, zone=NA,
     polygon(tmp@coords[indx,'x'], tmp@coords[indx,'y'], col=fillcol)
   }
 }
+
+
+Plot_range_quantiles = function( Data_Extrap, Report, TmbData, a_xl, NN_Extrap, Year_Set=NULL, Prob_vec=c(0.10,0.5,0.90), FileName_Quantiles=paste0(getwd(),"/Range_boundary.png") ){
+  # Default inputs
+  if( is.null(Year_Set)) Year_Set = 1:TmbData$n_t
+
+  # Plot range boundaries
+  calc_quantile = function(x, w, prob=0.5, doplot=FALSE){
+    DF = data.frame(x,w)[order(x),]
+    DF = cbind(DF, 'cumsum'=cumsum(DF$w)/sum(w))
+    Pred = rep(NA,length(prob))
+    for(i in 1:length(prob)){
+      Between = which( DF$cumsum>prob[i] )[1] - 1:0
+      Pred[i] = DF[Between[1],'x'] + (DF[Between[2],'x']-DF[Between[1],'x'])*(prob[i]-DF[Between[1],'cumsum'])/(DF[Between[2],'cumsum']-DF[Between[1],'cumsum'])
+    }
+    if(doplot==TRUE){
+      plot( x=DF$x, y=DF$cumsum, type="l" )
+      abline(v=Pred)
+    }
+    return( Pred )
+  }
+
+  # Extrapolation locations
+  Data_Extrap_Range = cbind( Data_Extrap[,c('Lat','Lon','N_km','E_km')], 'Include'=ifelse(Data_Extrap[,'Area_km2']>0, TRUE, FALSE) )
+  # Add and rescale density
+  for(t in 1:length(Year_Set)){
+    Data_Extrap_Range = cbind( Data_Extrap_Range, Report$Index_xtl[NN_Extrap$nn.idx,t,1] * (Data_Extrap[,'Area_km2'] / a_xl[NN_Extrap$nn.idx,1]) )
+    colnames( Data_Extrap_Range )[ncol(Data_Extrap_Range)] = paste0("Year_",Year_Set[t])
+    Data_Extrap_Range[,paste0("Year_",Year_Set[t])] = ifelse( Data_Extrap[,'Area_km2']==0 & a_xl[NN_Extrap$nn.idx,1]==0, 0, Data_Extrap_Range[,paste0("Year_",Year_Set[t])])
+  }
+
+  # Plot
+  png( file=FileName_Quantiles, width=6.5, height=6.5, res=200, units="in")
+    par( mfcol=c(2,2), mar=c(0,2,2,0), mgp=c(1.75,0.25,0), tck=-0.02, oma=c(4,0,0,0))
+    for(z in 1:4){
+      Range_Quantile = array( NA, dim=c(length(Year_Set),3), dimnames=list(NULL,c("min","mid","max")) )
+      for(t in 1:nrow(Range_Quantile)) Range_Quantile[t,] = calc_quantile( x=Data_Extrap_Range[,c('Lat','Lon','N_km','E_km')[z]], w=Data_Extrap_Range[,paste0('Year_',Year_Set[t])], prob=c(0.05,0.5,0.95), doplot=FALSE )
+      matplot( y=Range_Quantile, x=Year_Set, type="l", col="black", lty="solid", lwd=2, xlab="", ylab="", main=c('Latitude','Longitude','Nothings','Eastings')[z], xaxt="n" )
+      if(z %in% c(2,4)) axis(1)
+    }
+    mtext( side=1, text="Year", outer=TRUE, line=2)
+  dev.off()
+
+  # Return
+  Return = list( "Data_Extrap_Range"=Data_Extrap_Range )
+  return( invisible(Return) )
+}