diff --git a/DESCRIPTION b/DESCRIPTION
index c833dbf..cfe64ea 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,10 +1,15 @@
 Package: FishStatsUtils
 Type: Package
 Title: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox
-Version: 2.5.0
-Date: 2019-12-01
-Authors@R: person("James","Thorson", email="James.Thorson@noaa.gov", role=c("aut","cre"))
-Maintainer: James Thorson <James.Thorson@noaa.gov>
+Version: 2.6.0
+Date: 2020-04-11
+Authors@R: 
+    c(person(given = "James",
+             family = "Thorson",
+             role = c("aut", "cre"),
+             email = "James.Thorson@noaa.gov",
+             comment = c(ORCID = "0000-0001-7415-1010"))
+    )
 Description: FishStatsUtils contains utilities (shared code and data) used by multiple 
   packages (VAST, SpatialDeltaGLMM, MIST, Spatial_FA, SpatialDFA, surplus_production, EOFR) that are designed
   for spatio-temporal analysis of ecological data.  
@@ -29,17 +34,21 @@ Imports:
     MatrixModels,
     rgdal,
     ThorsonUtilities,
+    TMBhelper,
     abind,
     corpcor,
     pander,
     rnaturalearth,
-    formatR
+    rnaturalearthdata,
+    formatR,
+    splancs
 Depends:
     R (>= 3.1.0)
 Suggests:
     testthat
 Remotes: 
-    james-thorson/utilities
+    james-thorson/utilities,
+    kaskr/TMB_contrib_R/TMBhelper
 License: GPL-3
 LazyData: yes
 BuildVignettes: yes
diff --git a/NAMESPACE b/NAMESPACE
index 887e7f7..4fc0ef9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -36,6 +36,7 @@ export(calc_cov)
 export(calculate_proportion)
 export(combine_extrapolation_info)
 export(combine_lists)
+export(convert_shapefile)
 export(convert_version_name)
 export(fit_model)
 export(format_covariates)
@@ -70,4 +71,5 @@ export(plot_variable)
 export(project_coordinates)
 export(rotate_factors)
 export(sample_variable)
+export(simulate_data)
 export(summarize_covariance)
diff --git a/R/Geostat_Sim.R b/R/Geostat_Sim.R
deleted file mode 100644
index 4e02343..0000000
--- a/R/Geostat_Sim.R
+++ /dev/null
@@ -1,285 +0,0 @@
-
-#' Simulate data
-#'
-#' \code{Geostat_Sim} simulates data for use when testing \code{VAST}
-#'
-#' @param Sim_Settings an optional tagged list for specifying input parameters (see function code to determine settings)
-#' @inheritParams make_spatial_info
-#' @param Data_Geostat A data frame with column headers \code{c('Lon','Lat','Year','Vessel','AreaSwept_km2')} containing sample design to mimic
-#' @param standardize_fields Boolean, whether to ensure that random fields have sample mean and standard deviation equal to their inputted values
-
-#' @return Return Tagged list of output
-#' \describe{
-#'   \item{Data_Geostat}{Simulated data for analysis}
-#'   \item{B_tl}{True biomass for each year and stratum}
-#' }
-
-#' @examples
-#' ## Do not run (will be slow, due to simulating fine-scale spatial variation for many sites):
-#' ##
-#' ## # Prepare inputs
-#' ## Database = FishData::download_catch_rates( survey="Eastern_Bering_Sea", species_set="Gadus chalcogrammus", error_tol=0.01 )
-#' ## Data_Geostat = data.frame("Lat"=Database[,'Lat'], "Lon"=Database[,'Long'], "Year"=Database[,'Year'], "Vessel"="missing", "AreaSwept_km2"=1 )
-#' ## Extrapolation_List = FishStatsUtils::Prepare_Extrapolation_Data_Fn( Region="Eastern_Bering_Sea", strata.limits=data.frame( 'STRATA'="All_areas") )
-#' ##
-#' ## # Use function
-#' ## SimList = FishStatsUtils::Geostat_Sim(Sim_Settings=list(), Extrapolation_List, Data_Geostat=Data_Geostat )
-#' ## Data_Sim = SimList$Data_Geostat
-
-#' @export
-Geostat_Sim <-
-function(Sim_Settings, Extrapolation_List, Data_Geostat=NULL, MakePlot=FALSE, DateFile=paste0(getwd(),"/"), standardize_fields=FALSE ){
-  # Terminology
-  # O: Spatial component; E: spatiotemporal component
-  # O1: presence-absence; O2: positive catch rate
-  # P1 : linear predictor of presence/absence in link-space
-  # R1: predictor of presence/absence in natural-space (transformed out of link-space)
-  # v_i: vessel for sample i
-  # t_i: year for sample i
-
-  # Specify default values
-  Settings = list("beta1_mean"=0, "beta2_mean"=0, "beta1_slope"=0, "beta2_slope"=0, "beta1_sd"=0, "beta2_sd"=0, "Nyears"=10, "Nsamp_per_year"=600,
-    "Depth1_km"=0, "Depth1_km2"=0, "Dist1_sqrtkm"=0, "Depth2_km"=0, "Depth2_km2"=0, "Dist2_sqrtkm"=0,
-    "SigmaO1"=0.5, "SigmaO2"=0.5, "SigmaE1"=0.1, "SigmaE2"=0.1, "SigmaV1"=0, "SigmaV2"=0, "SigmaVY1"=0, "SigmaVY2"=0,
-    "Range1"=1000, "Range2"=500, "SigmaM"=1, "ObsModel"=c(2,0),
-    "Nages"=1, "M"=Inf, "K"=Inf, "Linf"=1, "W_alpha"=1, "W_beta"=3, "Selex_A50_mean"=0, "Selex_A50_sd"=0, "Selex_Aslope"=Inf )
-
-  # Replace defaults with provided values (if any)
-  for( i in seq_len(length(Sim_Settings)) ){
-    if(names(Sim_Settings)[i] %in% names(Settings)){
-      Settings[[match(names(Sim_Settings)[i],names(Settings))]] = Sim_Settings[[i]]
-    }
-  }
-
-  start_time = Sys.time()
-
-  # Local functions
-  RFsim = function(model, x, y, standardize=TRUE){
-    if( model_O1@par.general$var>0 ){
-      vec = RandomFields::RFsimulate(model=model, x=x, y=y)@data[,1]
-      if(standardize==TRUE) vec = (vec - mean(vec)) / sd(vec) * sqrt(model@par.general$var)
-    }else{
-      vec = rep(0,length(x))
-    }
-    return(vec)
-  }
-  rPoisGam = function( n=1, shape, scale, intensity ){
-    Num = rpois( n=n, lambda=intensity )
-    Biomass = rep(0, length=n)
-    for(i in which(Num>0) ) Biomass[i] = sum( rgamma(n=Num[i], shape=shape, scale=scale) )
-    return( Biomass )
-  }
-
-  # Initialize random-field objects
-  model_O1 = RandomFields::RMgauss(var=Settings[['SigmaO1']]^2, scale=Settings[['Range1']])
-  model_E1 = RandomFields::RMgauss(var=Settings[['SigmaE1']]^2, scale=Settings[['Range1']])
-  model_O2 = RandomFields::RMgauss(var=Settings[['SigmaO2']]^2, scale=Settings[['Range2']])
-  model_E2 = RandomFields::RMgauss(var=Settings[['SigmaE2']]^2, scale=Settings[['Range2']])
-
-  # Define sampling locations
-  if( is.null(Data_Geostat) ){
-    s_i = sample(1:nrow(Extrapolation_List$Data_Extrap), size=Settings[['Nyears']]*Settings[['Nsamp_per_year']]) #, nrow=Settings[['Nsamp_per_year']], ncol=Settings[['Nyears']])
-    t_i = rep( 1:Settings[['Nyears']], each=Settings[['Nsamp_per_year']])
-    v_i = rep(rep( 1:4, each=Settings[['Nsamp_per_year']]/4), Settings[['Nyears']])
-    w_i = rep(0.01, length(s_i))
-  }else{
-    NN_domain = RANN::nn2( data=Extrapolation_List$Data_Extrap[,c('Lon','Lat')], query=Data_Geostat[,c('Lon','Lat')], k=1)
-    s_i = NN_domain$nn.idx[,1]
-    t_i = Data_Geostat[,'Year'] - min(Data_Geostat[,'Year']) + 1
-    v_i = as.numeric(factor(Data_Geostat[,'Vessel']))
-    w_i = Data_Geostat[,'AreaSwept_km2']
-  }
-
-  # Duplicate locations if Nages>1
-  if( Settings[["Nages"]]>=2 ){
-    a_i = rep( 1:Settings[["Nages"]], times=length(s_i) )
-    s_i = rep( s_i, each=Settings[["Nages"]] )
-    t_i = rep( t_i, each=Settings[["Nages"]] )
-    v_i = rep( v_i, each=Settings[["Nages"]] )
-    w_i = rep( w_i, each=Settings[["Nages"]] )
-  }else{
-    a_i = rep( 1, times=length(s_i) )
-  }
-
-  # Generate locations, years, vessel for each sample i
-  loc_s = Extrapolation_List$Data_Extrap[,c('E_km','N_km','Lat','Lon')]
-  loc_i = loc_s[s_i,]
-
-  # Generate intercepts
-  exp_beta1_t = Settings[['beta1_mean']] + Settings[['beta1_slope']] * (1:max(t_i)-mean(1:max(t_i))/2)
-  beta1_t = rnorm( max(t_i), mean=exp_beta1_t, sd=Settings[['beta1_sd']] )
-  exp_beta2_t = Settings[['beta2_mean']] + Settings[['beta2_slope']] * (1:max(t_i)-mean(1:max(t_i))/2)
-  beta2_t = rnorm( max(t_i), mean=exp_beta2_t, sd=Settings[['beta2_sd']] )
-
-  # Simulate spatial and spatio-temporal components
-  O1_s = RFsim(model=model_O1, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
-  O2_s = RFsim(model=model_O2, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
-  message( "Finished variation that is constant over time" )
-  E1_st = E2_st = matrix(NA, nrow=nrow(loc_s), ncol=max(t_i) )
-  for(t in 1:max(t_i)){
-    E1_st[,t] = RFsim(model=model_E1, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
-    E2_st[,t] = RFsim(model=model_E2, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
-    message( "Finished variation for year ", t )
-  }
-
-  # Extract covariates
-  X_sj = array( 0, dim=c(nrow(Extrapolation_List$Data_Extrap),3), dimnames=list(NULL,c('Depth_km','Depth_km2','Rock_dist_')) )
-  for( colI in 1:ncol(X_sj)){
-    if(colnames(X_sj)[colI] %in% colnames(Extrapolation_List$Data_Extrap)) X_sj[,colI] = as.matrix(Extrapolation_List$Data_Extrap[,c('Depth_km','Depth_km2','Rock_dist_')[colI]])
-    X_sj = ifelse( is.na(X_sj), outer(rep(1,nrow(X_sj)),colMeans(X_sj,na.rm=TRUE)), X_sj)
-  }
-
-  # Simulate vessel effects
-  Vessel1_vy = array( rnorm( n=max(v_i)*max(t_i), mean=0, sd=Settings[['SigmaVY1']]), dim=c(max(v_i),max(t_i)) )
-  Vessel2_vy = array( rnorm( n=max(v_i)*max(t_i), mean=0, sd=Settings[['SigmaVY2']]), dim=c(max(v_i),max(t_i)) )
-  Vessel1_v = array( rnorm( n=max(v_i), mean=0, sd=Settings[['SigmaV1']]), dim=c(max(v_i)) )
-  Vessel2_v = array( rnorm( n=max(v_i), mean=0, sd=Settings[['SigmaV2']]), dim=c(max(v_i)) )
-
-  # Calculate covariate effect
-  eta1_s = as.matrix(X_sj) %*% unlist(Settings[c('Depth1_km','Depth1_km2','Dist1_sqrtkm')])
-  eta2_s = as.matrix(X_sj) %*% unlist(Settings[c('Depth2_km','Depth2_km2','Dist2_sqrtkm')])
-
-  # If Nages==1, then simulate as biomass-dynamic model
-  if( Settings[["Nages"]]==1 ){
-    # Calculate expected values
-    P1_st = P2_st = matrix(NA, nrow=nrow(loc_s), ncol=max(t_i) )
-    for(t in 1:max(t_i)){
-      P1_st[,t] = beta1_t[t] + O1_s + E1_st[,t] + eta1_s
-      P2_st[,t] = beta2_t[t] + O2_s + E2_st[,t] + eta2_s
-    }
-
-    # Calculate linear predictors
-    P1_i = P1_st[cbind(s_i,t_i)] + Vessel1_vy[cbind(v_i,t_i)] + Vessel1_v[v_i]
-    P2_i = P2_st[cbind(s_i,t_i)] + Vessel2_vy[cbind(v_i,t_i)] + Vessel2_v[v_i]
-
-    # Calculate true spatio-temporal biomass and annual abundance
-    if( Settings[['ObsModel']][2]==0 ){
-      B_ast = array( plogis(P1_st) * exp(P2_st) * outer(Extrapolation_List$Area_km2_x,rep(1,max(t_i))), dim=c(1,dim(P1_st)) )
-    }
-    if( Settings[['ObsModel']][2] %in% c(1,2) ){
-      B_ast = array( exp(P1_st) * exp(P2_st) * outer(Extrapolation_List$Area_km2_x,rep(1,max(t_i))), dim=c(1,dim(P2_st)) )
-    }
-    B_st = B_ast[1,,]
-
-    # Objects specific to Nages==1
-    Return = list( "P1_st"=P1_st, "P2_st"=P2_st )
-
-    # Calculate true spatio-temporal biomass and annual abundance
-    if( Settings[['ObsModel']][2]==0 ){
-      B_st = plogis(P1_st) * exp(P2_st) * outer(Extrapolation_List$Area_km2_x,rep(1,max(t_i)))
-    }
-    if( Settings[['ObsModel']][2] %in% c(1,2) ){
-      B_st = exp(P1_st) * exp(P2_st) * outer(Extrapolation_List$Area_km2_x,rep(1,max(t_i)))
-    }
-  }
-
-  # If Nages==1, then simulate as biomass-dynamic model
-  if( Settings[["Nages"]]>=2 ){
-    if( !(Settings[['ObsModel']][2] %in% c(1,2)) ) stop("If using age-structured simulation, please use 'Settings[['ObsModel']][2]=1'")
-
-    # Deviations in initial age-structure
-    E1_sa = matrix(NA, nrow=nrow(loc_s), ncol=Settings[["Nages"]] )
-    for(a in 2:Settings[["Nages"]]){
-      E1_sa[,a] = RFsim(model=model_E1, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
-      message( "Finished variation for age ", a, " in year 1" )
-    }
-
-    # Biomass at age
-    L_a = Settings[["Linf"]] * (1 - exp(-Settings[["K"]] * 1:Settings[["Nages"]]) )
-    W_a = Settings[["W_alpha"]] * L_a^Settings[["W_beta"]]
-
-    # Selectivity at age
-    Selex_A50_v = rnorm( max(v_i), mean=Settings[["Selex_A50_mean"]], sd=Settings[["Selex_A50_sd"]] )
-    Selex_av = NULL
-    for( v in 1:max(v_i)) Selex_av = cbind( Selex_av, plogis(1:Settings[["Nages"]], location=Selex_A50_v[v], scale=Settings[["Selex_Aslope"]]) )
-
-    # Numbers-density at age (N_ast) and Individual weight at age (W_ast) for age 1
-    N_ast = W_ast = array(NA, dim=c(Settings[["Nages"]],nrow(loc_s),max(t_i)) )
-    for(t in 1:max(t_i)){
-      # Individual weight at age
-      W_ast[,,t] = W_a %o% exp( beta2_t[t] + O2_s + E2_st[,t] + eta2_s )
-      # Numbers-density at age
-      N_ast[1,,t] = exp( beta1_t[t] + O1_s + E1_st[,t] + eta1_s )
-      for( a in 2:Settings[["Nages"]] ){
-        # Initial age-structure
-        if( t==1 ){
-          N_ast[a,,t] = exp( beta1_t[t] + O1_s + E1_sa[,a] + eta1_s - Settings[["M"]]*(a-1)  )
-        }
-        if( t>=2 ){
-          N_ast[a,,t] = exp(-Settings[["M"]]) * N_ast[a-1,,t-1]
-        }
-      }
-    }
-
-    # Calculate linear predictors
-    P1_i = log(Selex_av[cbind(a_i,v_i)] * N_ast[cbind(a_i,s_i,t_i)]) + Vessel1_vy[cbind(v_i,t_i)] + Vessel1_v[v_i]
-    P2_i = log(W_ast[cbind(a_i,s_i,t_i)]) + Vessel2_vy[cbind(v_i,t_i)] + Vessel2_v[v_i]
-
-    # Calculate true AVAILABLE spatio-temporal biomass and annual abundance
-    Selex_a = plogis(1:Settings[["Nages"]], location=Settings[["Selex_A50_mean"]], scale=Settings[["Selex_Aslope"]])
-    B_ast = N_ast * W_ast * (Selex_a %o% Extrapolation_List$Area_km2_x %o% rep(1,max(t_i)))
-    B_st = apply( B_ast, MARGIN=2:3, FUN=sum )
-    B_at = apply( B_ast, MARGIN=c(1,3), FUN=sum )
-
-    # Objects specific to Nages>=2
-    Return = list( "L_a"=L_a, "W_a"=W_a, "Selex_av"=Selex_av, "E1_sa"=E1_sa, "N_ast"=N_ast, "W_ast"=W_ast, "B_ast"=B_ast, "B_at"=B_at )
-  }
-
-  # Calculate linked-predictors
-  if( Settings[['ObsModel']][2]==0 ){
-    R1_i = plogis( P1_i )
-    R2_i = w_i * exp( P2_i )
-  }
-  if( Settings[['ObsModel']][2]==1 ){
-    R1_i = 1 - exp( -1*w_i*exp(P1_i) )
-    R2_i = w_i*exp(P1_i) / R1_i * exp(P2_i)
-  }
-  if( Settings[['ObsModel']][2]==2 ){
-    R1_i = w_i * exp(P1_i)
-    R2_i = exp(P2_i)
-  }
-
-  # Simulate catch and assemble data frame
-  if( Settings[['ObsModel']][1]==2 ){
-    b1_i = rbinom( n=length(R1_i), size=1, prob=R1_i )
-    b2_i = rlnorm( n=length(R2_i), meanlog=log(R2_i)-Settings[['SigmaM']]^2/2, sdlog=Settings[['SigmaM']])
-    b_i = b1_i * b2_i
-  }
-  if( Settings[['ObsModel']][1]==8 ){
-    if( Settings[['ObsModel']][2]!=2 ) stop("Must use 'ObsModel=c(8,2)'")
-    b_i = rep(NA,length(R1_i))
-    for(i in 1:length(b_i)) b_i[i] = rPoisGam( n=1, shape=Settings[['SigmaM']], scale=R2_i[i], intensity=R1_i[i] )
-  }
-  if( Settings[['ObsModel']][1]==10 ){
-    if( Settings[['ObsModel']][2]!=2 ) stop("Must use 'ObsModel=c(10,2)'")
-    b_i = rep(NA,length(R1_i))                            # R1_i(i)*R2_i(i), R1_i(i), invlogit(SigmaM(c_i(i),0))+1.0
-    for(i in 1:length(b_i)) b_i[i] = tweedie::rtweedie(n=1, mu=R1_i[i]*R2_i[i], phi=R1_i[i], power=plogis(Settings[['SigmaM']])+1.0)
-  }
-  Data_Geostat = cbind( "Catch_KG"=b_i, "Year"=t_i, "Vessel"=v_i, "Age"=a_i, "AreaSwept_km2"=w_i, "Lat"=loc_i[,'Lat'], "Lon"=loc_i[,'Lon'] )
-
-  # Calculate true center-of-gravity (COG)
-  COG_tm = array(NA, dim=c(max(t_i),4), dimnames=list(NULL,c("Lat","Lon","E_km","N_km")))
-  for( tI in 1:nrow(COG_tm)){
-  for( mI in 1:ncol(COG_tm)){
-    COG_tm[tI,mI] = weighted.mean( x=Extrapolation_List$Data_Extrap[,colnames(COG_tm)[mI]], w=B_st[,tI] )
-  }}
-
-  # Calculate stratified biomass
-  B_tl = array(NA, dim=c(max(t_i),ncol(Extrapolation_List$a_el)), dimnames=list(NULL,colnames(Extrapolation_List$a_el)))
-  for( l in 1:ncol(B_tl) ){
-    B_tl[,l] = colSums( B_st * outer(Extrapolation_List$a_el[,l]/Extrapolation_List$Area_km2_x,rep(1,max(t_i))), na.rm=TRUE )
-  }
-
-  # Timing
-  time_for_simulation = Sys.time()-start_time
-  message( "Total time: ", time_for_simulation )
-
-  # Objects shared by both
-  Return = c( Return, list("Data_Geostat"=data.frame(Data_Geostat), "B_tl"=B_tl, "B_st"=B_st, "COG_tm"=COG_tm, "beta1_t"=beta1_t,
-    "beta2_t"=beta2_t, "O1_s"=O1_s, "O2_s"=O2_s, "E1_st"=E1_st, "Vessel1_vy"=Vessel1_vy, "Vessel2_vy"=Vessel2_vy, "eta1_s"=eta1_s, "eta2_s"=eta2_s,
-    "Vessel1_v"=Vessel1_v, "Vessel2_v"=Vessel2_v, "E2_st"=E2_st, "s_i"=s_i, "t_i"=t_i, "a_i"=a_i, "w_i"=w_i,
-    "P1_i"=P1_i, "P2_i"=P2_i, "R1_i"=R1_i, "R2_i"=R2_i,
-    "time_for_simulation"=time_for_simulation, "Sim_Settings"=Settings) )
-  if( exists("b1_i")) Return = c(Return, list("b1_i"=b1_i, "b2_i"=b2_i))
-  return( Return )
-}
diff --git a/R/Heatmap_Legend.R b/R/Heatmap_Legend.R
deleted file mode 100644
index 26ad3aa..0000000
--- a/R/Heatmap_Legend.R
+++ /dev/null
@@ -1,13 +0,0 @@
-Heatmap_Legend <-
-function( colvec, heatrange, textmargin=NULL, labeltransform="uniform", dopar=TRUE, side=4 ){
-  if(dopar==TRUE) par( xaxs="i", yaxs="i", mar=c(1,0,1,2+ifelse(is.null(textmargin),0,1.5)), mgp=c(1.5,0.25,0), tck=-0.02 )
-  N = length(colvec)
-  Y = seq(heatrange[1], heatrange[2], length=N+1)
-  plot( 1, type="n", xlim=c(0,1), ylim=heatrange, xlab="", ylab="", main="", xaxt="n", yaxt="n", cex.main=1.5, xaxs="i", yaxs="i")
-  for( i in 1:N) polygon( x=c(0,1,1,0), y=Y[c(i,i,i+1,i+1)], col=colvec[i], border=NA)
-  if( labeltransform=="uniform" ) Labels = pretty(heatrange)
-  if( labeltransform=="inv_log10" ) Labels = 10^pretty(heatrange)
-  axis(side=4, at=pretty(heatrange), labels=Labels )
-  if(!is.null(textmargin)) mtext(side=side, text=textmargin, line=2, cex=1.5, las=0)
-  #mtext(side=1, outer=TRUE, line=1, "Legend")
-}
diff --git a/R/Prepare_XXXX_Extrapolation_Data_Fn.R b/R/Prepare_XXXX_Extrapolation_Data_Fn.R
index 5ef6f40..bd62074 100644
--- a/R/Prepare_XXXX_Extrapolation_Data_Fn.R
+++ b/R/Prepare_XXXX_Extrapolation_Data_Fn.R
@@ -355,7 +355,7 @@ function( strata.limits=NULL, epu_to_use = c('All', 'Georges_Bank','Mid_Atlantic
   }
   message("Using strata ", strata.limits)
 
-  if(tolower(epu_to_use) == "all") {
+  if(any(tolower(epu_to_use) %in% "all")) {
     epu_to_use <- c('Georges_Bank','Mid_Atlantic_Bight','Scotian_Shelf','Gulf_of_Maine','Other')
   }
 
@@ -368,10 +368,12 @@ function( strata.limits=NULL, epu_to_use = c('All', 'Georges_Bank','Mid_Atlantic
   if( length(strata.limits)==1 && strata.limits[1]=="EPU" ){
     # Specify epu by 'epu_to_use'
     Data_Extrap <- Data_Extrap[Data_Extrap$EPU %in% epu_to_use, ]
-    a_el = matrix(NA, nrow=nrow(Data_Extrap), ncol=length(unique(northwest_atlantic_grid[,'EPU'])), dimnames=list(NULL,unique(northwest_atlantic_grid[,'EPU'])) )
+    Data_Extrap$EPU <- droplevels(Data_Extrap$EPU)
+
+    a_el = matrix(NA, nrow=nrow(Data_Extrap), ncol=length(epu_to_use), dimnames=list(NULL, epu_to_use) )
     Area_km2_x = Data_Extrap[, "Area_in_survey_km2"]
     for(l in 1:ncol(a_el)){
-      a_el[,l] = ifelse( Data_Extrap[,'EPU']==unique(northwest_atlantic_grid[,'EPU'])[l], Area_km2_x, 0 )
+      a_el[,l] = ifelse( Data_Extrap[,'EPU'] %in% epu_to_use[[l]], Area_km2_x, 0 )
     }
   }else{
     # Specify strata by 'stratum_number'
diff --git a/R/calculate_proportion.R b/R/calculate_proportion.R
index 080ac54..b2df0db 100644
--- a/R/calculate_proportion.R
+++ b/R/calculate_proportion.R
@@ -21,7 +21,7 @@ calculate_proportion = function( TmbData, Index, Expansion_cz=NULL, Year_Set=NUL
   interval_width=1, width=6, height=6, xlab="Category", ylab="Proportion", ... ){
 
   # Warnings and errors
-  if( !all(TmbData[['FieldConfig']] %in% c(-2,-1)) ){
+  if( !all(TmbData[['FieldConfig']] %in% c(-3,-2,-1)) ){
     message("Derivation only included for independent categories")
     return( invisible("Not run") )
   }
@@ -39,7 +39,6 @@ calculate_proportion = function( TmbData, Index, Expansion_cz=NULL, Year_Set=NUL
   Index_tl = apply(Index_ctl,MARGIN=2:3,FUN=sum)
   SE_Index_tl = sqrt(apply(SE_Index_ctl^2,MARGIN=2:3,FUN=sum,na.rm=TRUE))
 
-
   # Approximate variance for proportions, and effective sample size
   Neff_ctl = var_Prop_ctl = array(NA,dim=dim(Prop_ctl))
   for( cI in 1:dim(var_Prop_ctl)[1]){
diff --git a/R/combine_lists.R b/R/combine_lists.R
index 36b8c8b..65e071d 100644
--- a/R/combine_lists.R
+++ b/R/combine_lists.R
@@ -9,7 +9,7 @@
 #' @return combined list.
 #'
 #' @export
-combine_lists = function( default, input, args_to_use=NULL ){
+combine_lists = function( default, input, args_to_use=NULL, use_partial_matching=FALSE ){
   output = default
   for( i in seq_along(input) ){
     if( names(input)[i] %in% names(default) ){
@@ -19,7 +19,14 @@ combine_lists = function( default, input, args_to_use=NULL ){
     }
   }
   if( !is.null(args_to_use) ){
-    output = output[intersect(names(output),args_to_use)]
+    # Exact matching
+    if( use_partial_matching==FALSE ){
+      output = output[ intersect(names(output),args_to_use) ]
+    }
+    # Partial matching
+    if( use_partial_matching==TRUE ){
+      stop( "`use_partial_matching=TRUE` is not implemented" )
+    }
   }
   return( output )
 }
diff --git a/R/convert_shapefile.R b/R/convert_shapefile.R
new file mode 100644
index 0000000..f453e9b
--- /dev/null
+++ b/R/convert_shapefile.R
@@ -0,0 +1,99 @@
+
+
+#' Convert shapefile to extrapolation-grid
+#'
+#' \code{convert_shapefile} reads in a shapefile and creates an extrapolation with a chosen resolution
+#'
+#' @inheritParams sp::CRS
+#' @inheritParams make_extraploation_info
+#'
+#' @param file_path path for shapefile on harddrive
+#' @param make_plots Boolean indicating whether to visualize inputs and outputs as maps
+#' @param projargs_for_shapefile projection-arguments (e.g., as parsed by \code{sp::CRS}), that are used when reading shapefile and overriding any projection arguments already saved there;  Default \code{projargs_for_shapefile=NULL} uses the projection arguments available in the shapefile
+#' @param projargs A character string of projection arguments used to project the shapefile prior to construction an extrapolation-grid; the arguments must be entered exactly as in the PROJ.4 documentation; Default \code{projargs=NULL} uses UTM and infers the UTM zone based on the centroid of the shapefile.
+
+#' @return extrapolation-grid
+
+#' @examples
+#' \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 )
+#' }
+
+#' @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 )
+  proj_orig = "+proj=longlat +ellps=WGS84 +no_defs"
+  shapefile_orig = sp::spTransform(shapefile_orig, CRSobj=sp::CRS(proj_orig) )
+  #shapefile_orig@proj4string = sp::CRS(proj_orig)
+
+  # Infer projargs if missing, and project
+  utm_zone = floor((mean(sp::bbox(shapefile_orig)[1,]) + 180) / 6) + 1
+  projargs_utm = paste0("+proj=utm +zone=",utm_zone," +ellps=WGS84 +datum=WGS84 +units=km +no_defs ")
+  if( is.null(projargs) || is.na(projargs) ){
+    projargs = projargs_utm
+  }
+  shapefile_proj = sp::spTransform(shapefile_orig, CRSobj=sp::CRS(projargs) )
+
+  # Determine bounds for box
+  bbox = shapefile_proj@bbox
+  bbox[,1] = floor(bbox[,1] / grid_dim_km[1]) * grid_dim_km[1]
+  bbox[,2] = ceiling(bbox[,2] / grid_dim_km[2]) * grid_dim_km[2]
+
+  # Make grid
+  grid_proj = expand.grid(X=seq(bbox[1,1],bbox[1,2],by=grid_dim_km[1]), Y=seq(bbox[2,1],bbox[2,2],by=grid_dim_km[2]))
+  grid_proj = sp::SpatialPointsDataFrame( coords=grid_proj[,c("X","Y")], data=grid_proj[,1:2], proj4string=sp::CRS(projargs) )
+
+  # Restrict points to those within a polygon
+  grid_proj@data = sp::over(grid_proj, shapefile_proj)
+  grid_proj = subset( grid_proj, !apply(grid_proj@data,MARGIN=1,FUN=function(vec){all(is.na(vec))}) )
+
+  # Return to original coordinates
+  grid_orig = sp::spTransform( grid_proj, CRSobj=sp::CRS(proj_orig) )
+  grid_orig = as.data.frame(grid_orig)
+
+  # Combine
+  extrapolation_grid = data.frame( "Lat"=grid_orig[,'Y'], "Lon"=grid_orig[,'X'], "Area_km2"=prod(grid_dim_km), # "Stratum"=grid_orig[,'AreaName'],
+    "Include"=1, "E_km"=grid_proj@coords[,'X'], "N_km"=grid_proj@coords[,'Y'] )
+
+  # make plots
+  if( make_plots==TRUE ){
+    par( mfrow=c(2,2), mar=c(3,3,1,1), mgp=c(2,0.5,0), tck=-0.02 )
+
+    # Plot projected shapefile
+    sp::plot( shapefile_proj, main="shapefile in original projection" )
+    axis(1); axis(2); box()
+
+    # Plot projected grid
+    sp::plot( grid_proj, main="Grid in new projection" )
+    axis(1); axis(2); box()
+
+    # plot original shapefile
+    sp::plot( shapefile_orig, main="shapefile in new projection" )
+    axis(1); axis(2); box()
+
+    # plot original shapefile
+    sp::plot( grid_orig[,c('X','Y')], main="Grid in original coordinates" )
+    axis(1); axis(2); box()
+  }
+
+  # Compare areas
+  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'] )
+
+  # Messages
+  Area_ratio = area_grid_proj / area_shapefile_orig
+  if( quiet==FALSE ){
+    message( "Area of projected extrapolation-grid is ", formatC(Area_ratio*100,format="f",digits=2), "% of the original shapefile area" )
+  }
+  if( Area_ratio>(1+area_tolerance) | Area_ratio<(1-area_tolerance) ){
+    stop( "Please use a different projection to decrease this conversion error; perhaps `projargs=", projargs_utm, "`" )
+  }
+
+  # Return output
+  Return = list( "extrapolation_grid"=extrapolation_grid, "projargs"=projargs )
+  return( Return )
+}
+
diff --git a/R/deprecated.R b/R/deprecated.R
index 9e5d1f7..3dc1dd2 100644
--- a/R/deprecated.R
+++ b/R/deprecated.R
@@ -1004,6 +1004,19 @@ function(MappingDetails, Mat, PlotDF, MapSizeRatio=c('Width(in)'=4,'Height(in)'=
   return( invisible(list("Par"=Par)) )
 }
 
+Heatmap_Legend <-
+function( colvec, heatrange, textmargin=NULL, labeltransform="uniform", dopar=TRUE, side=4 ){
+  if(dopar==TRUE) par( xaxs="i", yaxs="i", mar=c(1,0,1,2+ifelse(is.null(textmargin),0,1.5)), mgp=c(1.5,0.25,0), tck=-0.02 )
+  N = length(colvec)
+  Y = seq(heatrange[1], heatrange[2], length=N+1)
+  plot( 1, type="n", xlim=c(0,1), ylim=heatrange, xlab="", ylab="", main="", xaxt="n", yaxt="n", cex.main=1.5, xaxs="i", yaxs="i")
+  for( i in 1:N) polygon( x=c(0,1,1,0), y=Y[c(i,i,i+1,i+1)], col=colvec[i], border=NA)
+  if( labeltransform=="uniform" ) Labels = pretty(heatrange)
+  if( labeltransform=="inv_log10" ) Labels = 10^pretty(heatrange)
+  axis(side=4, at=pretty(heatrange), labels=Labels )
+  if(!is.null(textmargin)) mtext(side=side, text=textmargin, line=2, cex=1.5, las=0)
+  #mtext(side=1, outer=TRUE, line=1, "Legend")
+}
 
 #' @export
 plot_lines = function( x, y, ybounds, fn=lines, col_bounds="black", bounds_type="whiskers", border=NA,
@@ -1180,3 +1193,291 @@ format_covariates = function( Lat_e, Lon_e, t_e, Cov_ep, Extrapolation_List, Spa
   return( Return )
 }
 
+
+#' Simulate data
+#'
+#' \code{Geostat_Sim} simulates data for use when testing \code{VAST}
+#'
+#' @param Sim_Settings an optional tagged list for specifying input parameters (see function code to determine settings)
+#' @inheritParams make_spatial_info
+#' @param Data_Geostat A data frame with column headers \code{c('Lon','Lat','Year','Vessel','AreaSwept_km2')} containing sample design to mimic
+#' @param standardize_fields Boolean, whether to ensure that random fields have sample mean and standard deviation equal to their inputted values
+
+#' @return Return Tagged list of output
+#' \describe{
+#'   \item{Data_Geostat}{Simulated data for analysis}
+#'   \item{B_tl}{True biomass for each year and stratum}
+#' }
+
+#' @examples
+#' ## Do not run (will be slow, due to simulating fine-scale spatial variation for many sites):
+#' ##
+#' ## # Prepare inputs
+#' ## Database = FishData::download_catch_rates( survey="Eastern_Bering_Sea", species_set="Gadus chalcogrammus", error_tol=0.01 )
+#' ## Data_Geostat = data.frame("Lat"=Database[,'Lat'], "Lon"=Database[,'Long'], "Year"=Database[,'Year'], "Vessel"="missing", "AreaSwept_km2"=1 )
+#' ## Extrapolation_List = FishStatsUtils::Prepare_Extrapolation_Data_Fn( Region="Eastern_Bering_Sea", strata.limits=data.frame( 'STRATA'="All_areas") )
+#' ##
+#' ## # Use function
+#' ## SimList = FishStatsUtils::Geostat_Sim(Sim_Settings=list(), Extrapolation_List, Data_Geostat=Data_Geostat )
+#' ## Data_Sim = SimList$Data_Geostat
+
+#' @export
+Geostat_Sim <-
+function(Sim_Settings, Extrapolation_List, Data_Geostat=NULL, MakePlot=FALSE, DateFile=paste0(getwd(),"/"), standardize_fields=FALSE ){
+  # Terminology
+  # O: Spatial component; E: spatiotemporal component
+  # O1: presence-absence; O2: positive catch rate
+  # P1 : linear predictor of presence/absence in link-space
+  # R1: predictor of presence/absence in natural-space (transformed out of link-space)
+  # v_i: vessel for sample i
+  # t_i: year for sample i
+
+  # Warning
+  warning( "`Geostat_Sim` is deprecated, please use `obj$env$simulate(.)` instead for supported simulation method")
+
+  # Specify default values
+  Settings = list("beta1_mean"=0, "beta2_mean"=0, "beta1_slope"=0, "beta2_slope"=0, "beta1_sd"=0, "beta2_sd"=0, "Nyears"=10, "Nsamp_per_year"=600,
+    "Depth1_km"=0, "Depth1_km2"=0, "Dist1_sqrtkm"=0, "Depth2_km"=0, "Depth2_km2"=0, "Dist2_sqrtkm"=0,
+    "SigmaO1"=0.5, "SigmaO2"=0.5, "SigmaE1"=0.1, "SigmaE2"=0.1, "SigmaV1"=0, "SigmaV2"=0, "SigmaVY1"=0, "SigmaVY2"=0,
+    "Range1"=1000, "Range2"=500, "SigmaM"=1, "ObsModel"=c(2,0),
+    "Nages"=1, "M"=Inf, "K"=Inf, "Linf"=1, "W_alpha"=1, "W_beta"=3, "Selex_A50_mean"=0, "Selex_A50_sd"=0, "Selex_Aslope"=Inf )
+
+  # Replace defaults with provided values (if any)
+  for( i in seq_len(length(Sim_Settings)) ){
+    if(names(Sim_Settings)[i] %in% names(Settings)){
+      Settings[[match(names(Sim_Settings)[i],names(Settings))]] = Sim_Settings[[i]]
+    }
+  }
+
+  start_time = Sys.time()
+
+  # Local functions
+  RFsim = function(model, x, y, standardize=TRUE){
+    if( model_O1@par.general$var>0 ){
+      vec = RandomFields::RFsimulate(model=model, x=x, y=y)@data[,1]
+      if(standardize==TRUE) vec = (vec - mean(vec)) / sd(vec) * sqrt(model@par.general$var)
+    }else{
+      vec = rep(0,length(x))
+    }
+    return(vec)
+  }
+  rPoisGam = function( n=1, shape, scale, intensity ){
+    Num = rpois( n=n, lambda=intensity )
+    Biomass = rep(0, length=n)
+    for(i in which(Num>0) ) Biomass[i] = sum( rgamma(n=Num[i], shape=shape, scale=scale) )
+    return( Biomass )
+  }
+
+  # Initialize random-field objects
+  model_O1 = RandomFields::RMgauss(var=Settings[['SigmaO1']]^2, scale=Settings[['Range1']])
+  model_E1 = RandomFields::RMgauss(var=Settings[['SigmaE1']]^2, scale=Settings[['Range1']])
+  model_O2 = RandomFields::RMgauss(var=Settings[['SigmaO2']]^2, scale=Settings[['Range2']])
+  model_E2 = RandomFields::RMgauss(var=Settings[['SigmaE2']]^2, scale=Settings[['Range2']])
+
+  # Define sampling locations
+  if( is.null(Data_Geostat) ){
+    s_i = sample(1:nrow(Extrapolation_List$Data_Extrap), size=Settings[['Nyears']]*Settings[['Nsamp_per_year']]) #, nrow=Settings[['Nsamp_per_year']], ncol=Settings[['Nyears']])
+    t_i = rep( 1:Settings[['Nyears']], each=Settings[['Nsamp_per_year']])
+    v_i = rep(rep( 1:4, each=Settings[['Nsamp_per_year']]/4), Settings[['Nyears']])
+    w_i = rep(0.01, length(s_i))
+  }else{
+    NN_domain = RANN::nn2( data=Extrapolation_List$Data_Extrap[,c('Lon','Lat')], query=Data_Geostat[,c('Lon','Lat')], k=1)
+    s_i = NN_domain$nn.idx[,1]
+    t_i = Data_Geostat[,'Year'] - min(Data_Geostat[,'Year']) + 1
+    v_i = as.numeric(factor(Data_Geostat[,'Vessel']))
+    w_i = Data_Geostat[,'AreaSwept_km2']
+  }
+
+  # Duplicate locations if Nages>1
+  if( Settings[["Nages"]]>=2 ){
+    a_i = rep( 1:Settings[["Nages"]], times=length(s_i) )
+    s_i = rep( s_i, each=Settings[["Nages"]] )
+    t_i = rep( t_i, each=Settings[["Nages"]] )
+    v_i = rep( v_i, each=Settings[["Nages"]] )
+    w_i = rep( w_i, each=Settings[["Nages"]] )
+  }else{
+    a_i = rep( 1, times=length(s_i) )
+  }
+
+  # Generate locations, years, vessel for each sample i
+  loc_s = Extrapolation_List$Data_Extrap[,c('E_km','N_km','Lat','Lon')]
+  loc_i = loc_s[s_i,]
+
+  # Generate intercepts
+  exp_beta1_t = Settings[['beta1_mean']] + Settings[['beta1_slope']] * (1:max(t_i)-mean(1:max(t_i))/2)
+  beta1_t = rnorm( max(t_i), mean=exp_beta1_t, sd=Settings[['beta1_sd']] )
+  exp_beta2_t = Settings[['beta2_mean']] + Settings[['beta2_slope']] * (1:max(t_i)-mean(1:max(t_i))/2)
+  beta2_t = rnorm( max(t_i), mean=exp_beta2_t, sd=Settings[['beta2_sd']] )
+
+  # Simulate spatial and spatio-temporal components
+  O1_s = RFsim(model=model_O1, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
+  O2_s = RFsim(model=model_O2, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
+  message( "Finished variation that is constant over time" )
+  E1_st = E2_st = matrix(NA, nrow=nrow(loc_s), ncol=max(t_i) )
+  for(t in 1:max(t_i)){
+    E1_st[,t] = RFsim(model=model_E1, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
+    E2_st[,t] = RFsim(model=model_E2, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
+    message( "Finished variation for year ", t )
+  }
+
+  # Extract covariates
+  X_sj = array( 0, dim=c(nrow(Extrapolation_List$Data_Extrap),3), dimnames=list(NULL,c('Depth_km','Depth_km2','Rock_dist_')) )
+  for( colI in 1:ncol(X_sj)){
+    if(colnames(X_sj)[colI] %in% colnames(Extrapolation_List$Data_Extrap)) X_sj[,colI] = as.matrix(Extrapolation_List$Data_Extrap[,c('Depth_km','Depth_km2','Rock_dist_')[colI]])
+    X_sj = ifelse( is.na(X_sj), outer(rep(1,nrow(X_sj)),colMeans(X_sj,na.rm=TRUE)), X_sj)
+  }
+
+  # Simulate vessel effects
+  Vessel1_vy = array( rnorm( n=max(v_i)*max(t_i), mean=0, sd=Settings[['SigmaVY1']]), dim=c(max(v_i),max(t_i)) )
+  Vessel2_vy = array( rnorm( n=max(v_i)*max(t_i), mean=0, sd=Settings[['SigmaVY2']]), dim=c(max(v_i),max(t_i)) )
+  Vessel1_v = array( rnorm( n=max(v_i), mean=0, sd=Settings[['SigmaV1']]), dim=c(max(v_i)) )
+  Vessel2_v = array( rnorm( n=max(v_i), mean=0, sd=Settings[['SigmaV2']]), dim=c(max(v_i)) )
+
+  # Calculate covariate effect
+  eta1_s = as.matrix(X_sj) %*% unlist(Settings[c('Depth1_km','Depth1_km2','Dist1_sqrtkm')])
+  eta2_s = as.matrix(X_sj) %*% unlist(Settings[c('Depth2_km','Depth2_km2','Dist2_sqrtkm')])
+
+  # If Nages==1, then simulate as biomass-dynamic model
+  if( Settings[["Nages"]]==1 ){
+    # Calculate expected values
+    P1_st = P2_st = matrix(NA, nrow=nrow(loc_s), ncol=max(t_i) )
+    for(t in 1:max(t_i)){
+      P1_st[,t] = beta1_t[t] + O1_s + E1_st[,t] + eta1_s
+      P2_st[,t] = beta2_t[t] + O2_s + E2_st[,t] + eta2_s
+    }
+
+    # Calculate linear predictors
+    P1_i = P1_st[cbind(s_i,t_i)] + Vessel1_vy[cbind(v_i,t_i)] + Vessel1_v[v_i]
+    P2_i = P2_st[cbind(s_i,t_i)] + Vessel2_vy[cbind(v_i,t_i)] + Vessel2_v[v_i]
+
+    # Calculate true spatio-temporal biomass and annual abundance
+    if( Settings[['ObsModel']][2]==0 ){
+      B_ast = array( plogis(P1_st) * exp(P2_st) * outer(Extrapolation_List$Area_km2_x,rep(1,max(t_i))), dim=c(1,dim(P1_st)) )
+    }
+    if( Settings[['ObsModel']][2] %in% c(1,2) ){
+      B_ast = array( exp(P1_st) * exp(P2_st) * outer(Extrapolation_List$Area_km2_x,rep(1,max(t_i))), dim=c(1,dim(P2_st)) )
+    }
+    B_st = B_ast[1,,]
+
+    # Objects specific to Nages==1
+    Return = list( "P1_st"=P1_st, "P2_st"=P2_st )
+
+    # Calculate true spatio-temporal biomass and annual abundance
+    if( Settings[['ObsModel']][2]==0 ){
+      B_st = plogis(P1_st) * exp(P2_st) * outer(Extrapolation_List$Area_km2_x,rep(1,max(t_i)))
+    }
+    if( Settings[['ObsModel']][2] %in% c(1,2) ){
+      B_st = exp(P1_st) * exp(P2_st) * outer(Extrapolation_List$Area_km2_x,rep(1,max(t_i)))
+    }
+  }
+
+  # If Nages==1, then simulate as biomass-dynamic model
+  if( Settings[["Nages"]]>=2 ){
+    if( !(Settings[['ObsModel']][2] %in% c(1,2)) ) stop("If using age-structured simulation, please use 'Settings[['ObsModel']][2]=1'")
+
+    # Deviations in initial age-structure
+    E1_sa = matrix(NA, nrow=nrow(loc_s), ncol=Settings[["Nages"]] )
+    for(a in 2:Settings[["Nages"]]){
+      E1_sa[,a] = RFsim(model=model_E1, x=loc_s[,'E_km'], y=loc_s[,'N_km'], standardize=standardize_fields)
+      message( "Finished variation for age ", a, " in year 1" )
+    }
+
+    # Biomass at age
+    L_a = Settings[["Linf"]] * (1 - exp(-Settings[["K"]] * 1:Settings[["Nages"]]) )
+    W_a = Settings[["W_alpha"]] * L_a^Settings[["W_beta"]]
+
+    # Selectivity at age
+    Selex_A50_v = rnorm( max(v_i), mean=Settings[["Selex_A50_mean"]], sd=Settings[["Selex_A50_sd"]] )
+    Selex_av = NULL
+    for( v in 1:max(v_i)) Selex_av = cbind( Selex_av, plogis(1:Settings[["Nages"]], location=Selex_A50_v[v], scale=Settings[["Selex_Aslope"]]) )
+
+    # Numbers-density at age (N_ast) and Individual weight at age (W_ast) for age 1
+    N_ast = W_ast = array(NA, dim=c(Settings[["Nages"]],nrow(loc_s),max(t_i)) )
+    for(t in 1:max(t_i)){
+      # Individual weight at age
+      W_ast[,,t] = W_a %o% exp( beta2_t[t] + O2_s + E2_st[,t] + eta2_s )
+      # Numbers-density at age
+      N_ast[1,,t] = exp( beta1_t[t] + O1_s + E1_st[,t] + eta1_s )
+      for( a in 2:Settings[["Nages"]] ){
+        # Initial age-structure
+        if( t==1 ){
+          N_ast[a,,t] = exp( beta1_t[t] + O1_s + E1_sa[,a] + eta1_s - Settings[["M"]]*(a-1)  )
+        }
+        if( t>=2 ){
+          N_ast[a,,t] = exp(-Settings[["M"]]) * N_ast[a-1,,t-1]
+        }
+      }
+    }
+
+    # Calculate linear predictors
+    P1_i = log(Selex_av[cbind(a_i,v_i)] * N_ast[cbind(a_i,s_i,t_i)]) + Vessel1_vy[cbind(v_i,t_i)] + Vessel1_v[v_i]
+    P2_i = log(W_ast[cbind(a_i,s_i,t_i)]) + Vessel2_vy[cbind(v_i,t_i)] + Vessel2_v[v_i]
+
+    # Calculate true AVAILABLE spatio-temporal biomass and annual abundance
+    Selex_a = plogis(1:Settings[["Nages"]], location=Settings[["Selex_A50_mean"]], scale=Settings[["Selex_Aslope"]])
+    B_ast = N_ast * W_ast * (Selex_a %o% Extrapolation_List$Area_km2_x %o% rep(1,max(t_i)))
+    B_st = apply( B_ast, MARGIN=2:3, FUN=sum )
+    B_at = apply( B_ast, MARGIN=c(1,3), FUN=sum )
+
+    # Objects specific to Nages>=2
+    Return = list( "L_a"=L_a, "W_a"=W_a, "Selex_av"=Selex_av, "E1_sa"=E1_sa, "N_ast"=N_ast, "W_ast"=W_ast, "B_ast"=B_ast, "B_at"=B_at )
+  }
+
+  # Calculate linked-predictors
+  if( Settings[['ObsModel']][2]==0 ){
+    R1_i = plogis( P1_i )
+    R2_i = w_i * exp( P2_i )
+  }
+  if( Settings[['ObsModel']][2]==1 ){
+    R1_i = 1 - exp( -1*w_i*exp(P1_i) )
+    R2_i = w_i*exp(P1_i) / R1_i * exp(P2_i)
+  }
+  if( Settings[['ObsModel']][2]==2 ){
+    R1_i = w_i * exp(P1_i)
+    R2_i = exp(P2_i)
+  }
+
+  # Simulate catch and assemble data frame
+  if( Settings[['ObsModel']][1]==2 ){
+    b1_i = rbinom( n=length(R1_i), size=1, prob=R1_i )
+    b2_i = rlnorm( n=length(R2_i), meanlog=log(R2_i)-Settings[['SigmaM']]^2/2, sdlog=Settings[['SigmaM']])
+    b_i = b1_i * b2_i
+  }
+  if( Settings[['ObsModel']][1]==8 ){
+    if( Settings[['ObsModel']][2]!=2 ) stop("Must use 'ObsModel=c(8,2)'")
+    b_i = rep(NA,length(R1_i))
+    for(i in 1:length(b_i)) b_i[i] = rPoisGam( n=1, shape=Settings[['SigmaM']], scale=R2_i[i], intensity=R1_i[i] )
+  }
+  if( Settings[['ObsModel']][1]==10 ){
+    if( Settings[['ObsModel']][2]!=2 ) stop("Must use 'ObsModel=c(10,2)'")
+    b_i = rep(NA,length(R1_i))                            # R1_i(i)*R2_i(i), R1_i(i), invlogit(SigmaM(c_i(i),0))+1.0
+    for(i in 1:length(b_i)) b_i[i] = tweedie::rtweedie(n=1, mu=R1_i[i]*R2_i[i], phi=R1_i[i], power=plogis(Settings[['SigmaM']])+1.0)
+  }
+  Data_Geostat = cbind( "Catch_KG"=b_i, "Year"=t_i, "Vessel"=v_i, "Age"=a_i, "AreaSwept_km2"=w_i, "Lat"=loc_i[,'Lat'], "Lon"=loc_i[,'Lon'] )
+
+  # Calculate true center-of-gravity (COG)
+  COG_tm = array(NA, dim=c(max(t_i),4), dimnames=list(NULL,c("Lat","Lon","E_km","N_km")))
+  for( tI in 1:nrow(COG_tm)){
+  for( mI in 1:ncol(COG_tm)){
+    COG_tm[tI,mI] = weighted.mean( x=Extrapolation_List$Data_Extrap[,colnames(COG_tm)[mI]], w=B_st[,tI] )
+  }}
+
+  # Calculate stratified biomass
+  B_tl = array(NA, dim=c(max(t_i),ncol(Extrapolation_List$a_el)), dimnames=list(NULL,colnames(Extrapolation_List$a_el)))
+  for( l in 1:ncol(B_tl) ){
+    B_tl[,l] = colSums( B_st * outer(Extrapolation_List$a_el[,l]/Extrapolation_List$Area_km2_x,rep(1,max(t_i))), na.rm=TRUE )
+  }
+
+  # Timing
+  time_for_simulation = Sys.time()-start_time
+  message( "Total time: ", time_for_simulation )
+
+  # Objects shared by both
+  Return = c( Return, list("Data_Geostat"=data.frame(Data_Geostat), "B_tl"=B_tl, "B_st"=B_st, "COG_tm"=COG_tm, "beta1_t"=beta1_t,
+    "beta2_t"=beta2_t, "O1_s"=O1_s, "O2_s"=O2_s, "E1_st"=E1_st, "Vessel1_vy"=Vessel1_vy, "Vessel2_vy"=Vessel2_vy, "eta1_s"=eta1_s, "eta2_s"=eta2_s,
+    "Vessel1_v"=Vessel1_v, "Vessel2_v"=Vessel2_v, "E2_st"=E2_st, "s_i"=s_i, "t_i"=t_i, "a_i"=a_i, "w_i"=w_i,
+    "P1_i"=P1_i, "P2_i"=P2_i, "R1_i"=R1_i, "R2_i"=R2_i,
+    "time_for_simulation"=time_for_simulation, "Sim_Settings"=Settings) )
+  if( exists("b1_i")) Return = c(Return, list("b1_i"=b1_i, "b2_i"=b2_i))
+  return( Return )
+}
diff --git a/R/fit_model.R b/R/fit_model.R
index 4d6389a..00cc6bd 100644
--- a/R/fit_model.R
+++ b/R/fit_model.R
@@ -3,7 +3,32 @@
 #'
 #' \code{fit_model} fits a spatio-temporal model to data
 #'
-#' This function is the user-interface for the functions that determine the extrapolation-grid, define spatial objects, build covariates from a formula interface, assemble data, build model, estimate parameters, and check for obvious problems with the estimates.
+#' This function is the user-interface for the multiple mid-level functions that
+#' perform separate components of a spatio-temporal analysis:
+#' \itemize{
+#' \item determine the extrapolation-grid \code{\link{make_extrapolation_info}},
+#' \item define spatial objects \code{\link{make_spatial_info}},
+#' \item build covariates from a formula interface \code{\link{make_covariates}},
+#' \item assemble data \code{\link{make_data}},
+#' \item build model \code{\link{make_model}},
+#' \item estimate parameters \code{\link[TMBhelper]{fit_tmb}}, and
+#' \item check for obvious problems with the estimates \code{\link{check_fit}}.
+#' }
+#' Please see reference documetation for each of those functions (e.g., \code{?make_extrapolation_info}) to see a list of arguments used by each mid-level function.
+#'
+#' Specifically, the mid-level functions called by \code{fit_model(.)} look for arguments in the following order of precedence (from highest to lowest precedence):
+#' \enumerate{
+#' \item \code{fit_model(.)} prioritizes using named arguments passed directly to \code{fit_model(.)}. If arguments are passed this way, they are used instead of other options below.
+#' \item If an argument is not passed supplied directly to \code{fit_model(.)}, then \code{fit_model(.)} looks for elements in input \code{settings}, as typically created by \code{\link{make_settings}}.
+#' \item If an argument is not supplied via (1) or (2) above, then each mid-level function uses default values defined in those function arguments, e.g., see \code{args(make_extrapolation_info)} for defaults for function \code{make_extrapolation_info(.)}
+#' }
+#' Collectively, this order of precedence allows users to specify inputs for a specific project via input method (1), the package author to change defaults through changes in the settings
+#' defined for a given purpose in \code{make_settings(.)} via input method (2), while still defaulting to package defaults via option (3).
+#'
+#' Variables are indexed internally for locations \code{g}, categories \code{c}, and times \code{y}.
+#' Location index \code{g} represents Longitude-Latitude \code{fit$extrapolation_list$Data_Extrap[which(fit$spatial_list$g_e==g),c('Lon','Lat')]};
+#' Time index \code{y} represents time \code{fit$year_labels}; and
+#' Category \code{g} corresponds to values in \code{fit$data_list$g_i}.
 #'
 #' @inheritParams make_extrapolation_info
 #' @inheritParams make_spatial_info
@@ -11,15 +36,25 @@
 #' @inheritParams VAST::make_data
 #' @inheritParams VAST::make_model
 #' @inheritParams TMBhelper::fit_tmb
-#' @param settings Output from \code{make_settings}
+#' @param settings Output from \code{\link{make_settings}}
 #' @param run_model Boolean indicating whether to run the model or simply return the inputs and built TMB object
 #' @param test_fit Boolean indicating whether to apply \code{VAST::check_fit} before calculating standard errors, to test for parameters hitting bounds etc; defaults to TRUE
-#' @param ... additional arguments to pass to \code{FishStatsUtils::make_extrapolation_info}, \code{FishStatsUtils::make_spatial_info}, \code{VAST::make_data}, \code{VAST::make_model}, or \code{TMBhelper::fit_tmb}, where arguments are matched by name against each function.  If an argument doesn't match, it is still passed to \code{VAST::make_data}
+#' @param ... additional arguments to pass to \code{\link{make_extrapolation_info}}, \code{\link{make_spatial_info}}, \code{\link[VAST]{make_data}}, \code{\link[VAST]{make_model}}, or \code{\link[TMBhelper]{fit_tmb}}, where arguments are matched by name against each function.  If an argument doesn't match, it is still passed to \code{\link[VAST]{make_data}}
 #'
-#' @return Returns a tagged list of internal objects, the TMB object, and slot \code{parameter_estimates} containing the MLE estimates
+#' @return Object of class \code{fit_model}, containing formatted inputs and outputs from VAST
+#' \describe{
+#'   \item{parameter_estimates}{Output from \code{\link[TMBhelper]{fit_tmb}}; see that documentation for definition of contents}
+#'   \item{extrapolation_list}{Output from \code{\link{make_extrapolation_info}}; see that documentation for definition of contents}
+#'   \item{spatial_list}{Output from \code{\link{make_spatial_info}}; see that documentation for definition of contents}
+#'   \item{data_list}{Output from \code{\link[VAST]{make_data}}; see that documentation for definition of contents}
+#'   \item{tmb_list}{Output from \code{\link[VAST]{make_model}}; see that documentation for definition of contents}
+#'   \item{ParHat}{Tagged list of maximum likelihood estimatesion of fixed effects and empirical Bayes estimates of random effects, following format of initial values generated by \code{\link[VAST]{make_parameters}}; see that documentation for definition of contents}
+#'   \item{Report}{Tagged list of VAST outputs. For example, estimated density for grid \code{g}, category \code{c}, and time \code{y} is available as \code{fit$Report$D_gcy[g,c,y]}; see Details section for description of indexing}
+#' }
 #'
 #' @family wrapper functions
 #' @seealso \code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots
+#' @seealso \code{\link{summary.fit_model}} for methods to summarize output, including obtain a dataframe of estimated densities
 #'
 #' @examples
 #' \dontrun{
@@ -47,10 +82,12 @@
 #' }
 #'
 #' @export
+#' @md
+# Using https://cran.r-project.org/web/packages/roxygen2/vignettes/rd-formatting.html for guidance on markdown-enabled documentation
 fit_model = function( settings, Lat_i, Lon_i, t_iz, b_i, a_i, c_iz=rep(0,length(b_i)),
   v_i=rep(0,length(b_i)), working_dir=paste0(getwd(),"/"),
   Xconfig_zcp=NULL, covariate_data, formula=~0, Q_ik=NULL, newtonsteps=1,
-  silent=TRUE, run_model=TRUE, test_fit=TRUE, ... ){
+  silent=TRUE, build_model=TRUE, run_model=TRUE, test_fit=TRUE, ... ){
 
   # Capture extra arguments to function
   extra_args = list(...)
@@ -70,58 +107,73 @@ fit_model = function( settings, Lat_i, Lon_i, t_iz, b_i, a_i, c_iz=rep(0,length(
 
   # Build extrapolation grid
   message("\n### Making extrapolation-grid")
-  extrapolation_args_default = list(Region=settings$Region, strata.limits=settings$strata.limits, zone=settings$zone)
-  extrapolation_args_input = extra_args[intersect(names(extra_args),formalArgs(make_extrapolation_info))]
-  extrapolation_args_input = combine_lists( input=extrapolation_args_input, default=extrapolation_args_default )
+  extrapolation_args_default = list(Region=settings$Region, strata.limits=settings$strata.limits, zone=settings$zone,
+    max_cells=settings$max_cells)
+  extrapolation_args_input = combine_lists( input=extra_args, default=extrapolation_args_default, args_to_use=formalArgs(make_extrapolation_info) )
   extrapolation_list = do.call( what=make_extrapolation_info, args=extrapolation_args_input )
 
   # Build information regarding spatial location and correlation
   message("\n### Making spatial information")
   spatial_args_default = list(grid_size_km=settings$grid_size_km, n_x=settings$n_x, Method=settings$Method, Lon_i=Lon_i, Lat_i=Lat_i,
-    Extrapolation_List=extrapolation_list, DirPath=working_dir, Save_Results=TRUE, fine_scale=settings$fine_scale)
-  spatial_args_input = extra_args[intersect(names(extra_args),formalArgs(make_spatial_info))]
-  spatial_args_input = combine_lists( input=spatial_args_input, default=spatial_args_default )
+    Extrapolation_List=extrapolation_list, DirPath=working_dir, Save_Results=TRUE, fine_scale=settings$fine_scale, knot_method=settings$knot_method)
+  spatial_args_input = combine_lists( input=extra_args, default=spatial_args_default, args_to_use=formalArgs(make_spatial_info) )
   spatial_list = do.call( what=make_spatial_info, args=spatial_args_input )
 
   # Build data
+  # Do *not* restrict inputs to formalArgs(make_data) because other potential inputs are still parsed by make_data for backwards compatibility
   message("\n### Making data object") # VAST::
   if(missing(covariate_data)) covariate_data = NULL
   data_args_default = list("Version"=settings$Version, "FieldConfig"=settings$FieldConfig, "OverdispersionConfig"=settings$OverdispersionConfig,
     "RhoConfig"=settings$RhoConfig, "VamConfig"=settings$VamConfig, "ObsModel"=settings$ObsModel, "c_iz"=c_iz, "b_i"=b_i, "a_i"=a_i, "v_i"=v_i,
     "s_i"=spatial_list$knot_i-1, "t_iz"=t_iz, "spatial_list"=spatial_list, "Options"=settings$Options, "Aniso"=settings$use_anisotropy,
     Xconfig_zcp=Xconfig_zcp, covariate_data=covariate_data, formula=formula, Q_ik=Q_ik)
-  data_args_input = combine_lists( input=extra_args, default=data_args_default )
+  data_args_input = combine_lists( input=extra_args, default=data_args_default )  # Do *not* use args_to_use
   data_list = do.call( what=make_data, args=data_args_input )
 
   # Build object
   message("\n### Making TMB object")
   model_args_default = list("TmbData"=data_list, "RunDir"=working_dir, "Version"=settings$Version,
-    "RhoConfig"=settings$RhoConfig, "loc_x"=spatial_list$loc_x, "Method"=spatial_list$Method)
-  model_args_input = extra_args[intersect(names(extra_args),formalArgs(make_model))]
-  model_args_input = combine_lists( input=model_args_input, default=model_args_default )
+    "RhoConfig"=settings$RhoConfig, "loc_x"=spatial_list$loc_x, "Method"=spatial_list$Method, "build_model"=build_model)
+  model_args_input = combine_lists( input=extra_args, default=model_args_default, args_to_use=formalArgs(make_model) )
   tmb_list = do.call( what=make_model, args=model_args_input )
-  if(silent==TRUE) tmb_list$Obj$env$beSilent()
 
   # Run the model or optionally don't
-  if( run_model==FALSE ){
+  if( run_model==FALSE | build_model==FALSE ){
     # Build and output
+    input_args = list( "extra_args"=extra_args, "extrapolation_args_input"=extrapolation_args_input,
+      "model_args_input"=model_args_input, "spatial_args_input"=spatial_args_input,
+      "data_args_input"=data_args_input )
     Return = list("data_frame"=data_frame, "extrapolation_list"=extrapolation_list, "spatial_list"=spatial_list,
       "data_list"=data_list, "tmb_list"=tmb_list, "year_labels"=year_labels, "years_to_plot"=years_to_plot,
-      "settings"=settings)
+      "settings"=settings, "input_args"=input_args)
     class(Return) = "fit_model"
     return(Return)
   }
+  if(silent==TRUE) tmb_list$Obj$env$beSilent()
+
+  # Check for obvious problems with model
+  if( test_fit==TRUE ){
+    message("\n### Testing model at initial values")
+    LogLike0 = tmb_list$Obj$fn( tmb_list$Obj$par )
+    Gradient0 = tmb_list$Obj$gr( tmb_list$Obj$par )
+    if( any( Gradient0==0 ) ){
+      message("\n")
+      stop("Please check model structure; some parameter has a gradient of zero at starting values\n", call.=FALSE)
+    }else{
+      message("Looks good: All fixed effects have a nonzero gradient")
+    }
+  }
 
   # Optimize object
   message("\n### Estimating parameters")
   # have user override upper, lower, and loopnum
   optimize_args_default1 = combine_lists( default=list(lower=tmb_list$Lower, upper=tmb_list$Upper, loopnum=2),
-    input=extra_args[intersect(names(extra_args),formalArgs(TMBhelper::fit_tmb))] )
+    input=extra_args, args_to_use=formalArgs(TMBhelper::fit_tmb) )
   # auto-override user inputs for optimizer-related inputs for first test run
   optimize_args_input1 = list(obj=tmb_list$Obj, savedir=NULL, newtonsteps=0, bias.correct=FALSE,
     control=list(eval.max=10000,iter.max=10000,trace=1), quiet=TRUE, getsd=FALSE )
   # combine
-  optimize_args_input1 = combine_lists( default=optimize_args_default1, input=optimize_args_input1 )
+  optimize_args_input1 = combine_lists( default=optimize_args_default1, input=optimize_args_input1, args_to_use=formalArgs(TMBhelper::fit_tmb) )
   parameter_estimates = do.call( what=TMBhelper::fit_tmb, args=optimize_args_input1 )
 
   # Check fit of model (i.e., evidence of non-convergence based on bounds, approaching zero, etc)
@@ -129,7 +181,7 @@ fit_model = function( settings, Lat_i, Lon_i, t_iz, b_i, a_i, c_iz=rep(0,length(
     problem_found = VAST::check_fit( parameter_estimates )
     if( problem_found==TRUE ){
       message("\n")
-      stop("Please change model structure to avoid problems with parameter estimates and then re-try\n", call.=FALSE)
+      stop("Please change model structure to avoid problems with parameter estimates and then re-try; see details in `?check_fit`\n", call.=FALSE)
     }
   }
 
@@ -138,11 +190,9 @@ fit_model = function( settings, Lat_i, Lon_i, t_iz, b_i, a_i, c_iz=rep(0,length(
     savedir=working_dir, bias.correct=settings$bias.correct, newtonsteps=newtonsteps,
     bias.correct.control=list(sd=FALSE, split=NULL, nsplit=1, vars_to_correct=settings$vars_to_correct),
     control=list(eval.max=10000,iter.max=10000,trace=1), loopnum=1)
-  # user over-rides all default inputs
-  optimize_args_input2 = extra_args[intersect(names(extra_args),formalArgs(TMBhelper::fit_tmb))]
-  # combine
-  optimize_args_input2 = combine_lists( input=optimize_args_input2, default=optimize_args_default2 )
-  # start from MLE
+  # combine while over-riding defaults using user inputs
+  optimize_args_input2 = combine_lists( input=extra_args, default=optimize_args_default2, args_to_use=formalArgs(TMBhelper::fit_tmb) )
+  # over-ride inputs to start from previous MLE
   optimize_args_input2 = combine_lists( input=list(startpar=parameter_estimates$par), default=optimize_args_input2 )
   parameter_estimates = do.call( what=TMBhelper::fit_tmb, args=optimize_args_input2 )
 
@@ -153,7 +203,8 @@ fit_model = function( settings, Lat_i, Lon_i, t_iz, b_i, a_i, c_iz=rep(0,length(
   # Build and output
   input_args = list( "extra_args"=extra_args, "extrapolation_args_input"=extrapolation_args_input,
     "model_args_input"=model_args_input, "spatial_args_input"=spatial_args_input,
-    "optimize_args_input1"=optimize_args_input1, "optimize_args_input2"=optimize_args_input2)
+    "optimize_args_input1"=optimize_args_input1, "optimize_args_input2"=optimize_args_input2,
+    "data_args_input"=data_args_input )
   Return = list("data_frame"=data_frame, "extrapolation_list"=extrapolation_list, "spatial_list"=spatial_list,
     "data_list"=data_list, "tmb_list"=tmb_list, "parameter_estimates"=parameter_estimates, "Report"=Report,
     "ParHat"=ParHat, "year_labels"=year_labels, "years_to_plot"=years_to_plot, "settings"=settings,
@@ -222,9 +273,14 @@ plot.fit_model <- function(x, what="results", ...)
 #' Extract summary of spatial estimates
 #'
 #' @title Extract spatial estimates
-#' @param fit Output from \code{\link{fit_model}}
+#' @param x Output from \code{\link{fit_model}}
 #' @param what Boolean indicating what to summarize; only option is `density`
 #' @param ... Not used
+#'
+#' \code{what="density"} returns a tagged list containing element \code{Density_dataframe},
+#' which lists the estimated density for every Latitude-Longitude-Year-Category combination
+#' for every modelled location in the extrapolation-grid.
+#'
 #' @return NULL
 #' @method summary fit_model
 #' @export
diff --git a/R/make_extrapolation_info.R b/R/make_extrapolation_info.R
index 9aec469..f018248 100644
--- a/R/make_extrapolation_info.R
+++ b/R/make_extrapolation_info.R
@@ -3,23 +3,38 @@
 #'
 #' \code{make_extrapolation_data} builds an object used to determine areas to extrapolation densities to when calculating indices
 #'
-#' To do area-weighted extrapolation of estimated density for use in calculating abundance indices, it is necessary to have a precise measurement of the footprint for a given survey design. Using VAST, analysts do this by including an "extrapolation grid" where densities are predicted at the location of each grid cell and where each grid cell is associated with a known area within a given survey design. Collaborators have worked with the package author to include the extrapolation-grid for several regions automatically in FishStatsUtils, but for new regions an analyst must either detect the grid automatically using \code{Region="Other"} or input an extrapolation-grid manually using \code{Region="User"}.  The extrapolation is also used to determine where to drawn pixels when plotting predictions of density.
+#' To do area-weighted extrapolation of estimated density for use in calculating abundance indices,
+#' it is necessary to have a precise measurement of the footprint for a given survey design.
+#' Using VAST, analysts do this by including an "extrapolation grid" where densities are predicted
+#' at the location of each grid cell and where each grid cell is associated with a known area within a given survey design.
+#' Collaborators have worked with the package author to include the extrapolation-grid for several
+#' regions automatically in FishStatsUtils. For new regions an analyst can either (1) detect
+#' the grid automatically using \code{Region="Other"}, or (2) input an extrapolation-grid manually
+#' using \code{Region="User"}, or supply a GIS shapefile \code{Region="[directory_path/file_name].shp"}.
+#' The extrapolation is also used to determine where to drawn pixels when plotting predictions of density.
+#' If a user supplies a character-vector with more than one of these, then they are combined to
+#' assemble a combined extrapolation-grid.
+#'
+#' When supplying a shapefile, I recommend using a UTM projection for projargs, which appears to have lower
+#' projection errors regarding total area than rnaturalearth.
 #'
 #' @inheritParams sp::CRS
 #' @inheritParams Calc_Kmeans
+#' @inheritParams convert_shapefile
 #'
-#' @param Region a character vector, where each element that is matched against potential values to determine the region for the extrapolation grid. Current options are "california_current", "west_coast_hook_and_line", "british_columbia", "eastern_bering_sea", "northern_bering_sea", "bering_sea_slope", "st_matthews_island", "aleutian_islands", "gulf_of_alaska", "northwest_atlantic", "south_africa", "gulf_of_st_lawrence", "new_zealand", "habcam", "gulf_of_mexico", "stream_network", "user", or "other"
+#' @param Region a character vector, where each element that is matched against potential values to determine the region for the extrapolation grid. Current options are "california_current", "west_coast_hook_and_line", "british_columbia", "eastern_bering_sea", "northern_bering_sea", "bering_sea_slope", "st_matthews_island", "aleutian_islands", "gulf_of_alaska", "northwest_atlantic", "south_africa", "gulf_of_st_lawrence", "new_zealand", "habcam", "gulf_of_mexico", "ATL-IBTS-Q1", "ATL-IBTS-Q4", "BITS", "BTS", "BTS-VIIA", "EVHOE", "IE-IGFS", "NIGFS", "NS_IBTS", "PT-IBTS", "SP-ARSA", "SP-NORTH", "SP-PORC", "stream_network", "user", "other", or the absolute path and file name for a GIS shapefile
 #' @param strata.limits an input for determining stratification of indices (see example script)
 #' @param zone UTM zone used for projecting Lat-Lon to km distances; use \code{zone=NA} by default to automatically detect UTM zone from the location of extrapolation-grid samples
 #' @param flip_around_dateline used applies when using UTM projection, where {flip_around_dateline=TRUE} causes code to convert given latitude on other side of globe (as helpful when data straddle dateline); default value depends upon \code{Region} used
 #' @param create_strata_per_region Boolean indicating whether to create a single stratum for all regions listed in \code{Region} (the default), or a combined stratum in addition to a stratum for each individual Region
-#' @param observations_LL a matrix with two columns (labeled 'Lat' and 'Lon') giving latitude and longitude for each observation; only used when \code{Region="user"}
-#' @param input_grid a matrix with three columns (labeled 'Lat', 'Lon', and 'Area_km2') giving latitude, longitude, and area for each cell of a user-supplied grid; only used when \code{Region="other"}
+#' @param observations_LL a matrix with two columns (labeled 'Lat' and 'Lon') giving latitude and longitude for each observation; only used when \code{Region="other"}
+#' @param input_grid a matrix with three columns (labeled 'Lat', 'Lon', and 'Area_km2') giving latitude, longitude, and area for each cell of a user-supplied grid; only used when \code{Region="user"}
 #' @param grid_dim_km numeric-vector with length two, giving the distance in km between cells in the automatically generated extrapolation grid; only used if \code{Region="other"}
 #' @param maximum_distance_from_sample maximum distance that an extrapolation grid cell can be from the nearest sample and still be included in area-weighted extrapolation of density; only used if \code{Region="other"}
 #' @param grid_dim_LL same as \code{grid_dim_km} except measured in latitude-longitude coordinates; only used if \code{Region="other"}
 #' @param grid_in_UTM Boolean stating whether to automatically generate an extrapolation grid based on sampling locations in km within the UTM projection of within Lat-Lon coordinates; only used if \code{Region="other"}
 #' @param strata_to_use strata to include by default for the BC coast extrapolation grid; only used if \code{Region="british_columbia"}
+#' @param epu_to_use EPU to include for the Northwest Atlantic (NWA) extrapolation grid, default is "All"; only used if \code{Region="northwest_atlantic"}
 #' @param survey survey to use for New Zealand extrapolation grid; only used if \code{Region="new_zealand"}
 #' @param region which coast to use for South Africa extrapolation grid; only used if \code{Region="south_africa"}
 #' @param surveyname area of West Coast to include in area-weighted extrapolation for California Current; only used if \code{Region="california_current"}
@@ -37,14 +52,19 @@
 
 #' @export
 make_extrapolation_info = function( Region, projargs=NA, zone=NA, strata.limits=data.frame('STRATA'="All_areas"),
-  create_strata_per_region=FALSE, max_cells=Inf, input_grid=NULL, observations_LL=NULL, grid_dim_km=c(2,2),
+  create_strata_per_region=FALSE, max_cells=NULL, input_grid=NULL, observations_LL=NULL, grid_dim_km=c(2,2),
   maximum_distance_from_sample=NULL, grid_in_UTM=TRUE, grid_dim_LL=c(0.1,0.1),
   region=c("south_coast","west_coast"), strata_to_use=c('SOG','WCVI','QCS','HS','WCHG'),
-  survey="Chatham_rise", surveyname='propInWCGBTS', flip_around_dateline, nstart=100, ... ){
+  epu_to_use=c('All','Georges_Bank','Mid_Atlantic_Bight','Scotian_Shelf','Gulf_of_Maine','Other')[1],
+  survey="Chatham_rise", surveyname='propInWCGBTS', flip_around_dateline, nstart=100,
+  area_tolerance=0.05, ... ){
 
   # Note: flip_around_dateline must appear in arguments for argument-matching in fit_model
   # However, it requires a different default value for different regions; hence the input format being used.
 
+  # Backwards compatibility for when settings didn't include max_cells, such that settings$max_cells=NULL
+  if(is.null(max_cells)) max_cells = Inf
+
   for( rI in seq_along(Region) ){
     Extrapolation_List = NULL
     if( tolower(Region[rI]) == "california_current" ){
@@ -85,7 +105,7 @@ make_extrapolation_info = function( Region, projargs=NA, zone=NA, strata.limits=
     }
     if( tolower(Region[rI]) == "northwest_atlantic" ){
       if(missing(flip_around_dateline)) flip_around_dateline = FALSE
-      Extrapolation_List = Prepare_NWA_Extrapolation_Data_Fn( strata.limits=strata.limits, projargs=projargs, zone=zone, flip_around_dateline=flip_around_dateline, ... )
+      Extrapolation_List = Prepare_NWA_Extrapolation_Data_Fn( strata.limits=strata.limits, epu_to_use=epu_to_use, projargs=projargs, zone=zone, flip_around_dateline=flip_around_dateline, ... )
     }
     if( tolower(Region[rI]) == "south_africa" ){
       if(missing(flip_around_dateline)) flip_around_dateline = FALSE
@@ -107,6 +127,18 @@ make_extrapolation_info = function( Region, projargs=NA, zone=NA, strata.limits=
       if(missing(flip_around_dateline)) flip_around_dateline = FALSE
       Extrapolation_List = Prepare_GOM_Extrapolation_Data_Fn( strata.limits=strata.limits, projargs=projargs, zone=zone, flip_around_dateline=flip_around_dateline, ... )
     }
+    if( toupper(Region[rI]) %in% c("ATL-IBTS-Q1","ATL-IBTS-Q4","BITS","BTS","BTS-VIIA","EVHOE","IE-IGFS","NIGFS","NS_IBTS","PT-IBTS","SP-ARSA","SP-NORTH","SP-PORC") ){
+      if( Region[rI]=="SP-ARSA" ) stop("There's some problem with `SP-ARSA` which precludes it's use")
+      Conversion = convert_shapefile( file_path=paste0(system.file("region_shapefiles",package="FishStatsUtils"),"/",toupper(Region[rI]),"/Shapefile.shp"),
+        projargs_for_shapefile="+proj=longlat +ellps=WGS84 +no_defs", projargs=projargs, grid_dim_km=grid_dim_km, area_tolerance=area_tolerance, ... )
+      Extrapolation_List = list( "a_el"=matrix(Conversion$extrapolation_grid[,'Area_km2'],ncol=1), "Data_Extrap"=Conversion$extrapolation_grid,
+        "zone"=NA, "projargs"=Conversion$projargs, "flip_around_dateline"=FALSE, "Area_km2_x"=Conversion$extrapolation_grid[,'Area_km2'])
+    }
+    if( file.exists(Region[rI]) ){
+      Conversion = convert_shapefile( file_path=Region[rI], projargs=projargs, grid_dim_km=grid_dim_km, area_tolerance=area_tolerance, ... )
+      Extrapolation_List = list( "a_el"=matrix(Conversion$extrapolation_grid[,'Area_km2'],ncol=1), "Data_Extrap"=Conversion$extrapolation_grid,
+        "zone"=NA, "projargs"=Conversion$projargs, "flip_around_dateline"=FALSE, "Area_km2_x"=Conversion$extrapolation_grid[,'Area_km2'])
+    }
     if( tolower(Region[rI]) == "stream_network" ){
       if( is.null(input_grid)){
         stop("Because you're using a stream network, please provide 'input_grid' input")
@@ -129,7 +161,7 @@ make_extrapolation_info = function( Region, projargs=NA, zone=NA, strata.limits=
     }
     if( is.null(Extrapolation_List) ){
       if( is.null(observations_LL)){
-        stop("Because you're using a new Region[rI], please provide 'observations_LL' input")
+        stop("Because you're using a new Region[rI], please provide 'observations_LL' input with columns named `Lat` and `Lon`")
       }
       if(missing(flip_around_dateline)) flip_around_dateline = FALSE
       Extrapolation_List = Prepare_Other_Extrapolation_Data_Fn( strata.limits=strata.limits, observations_LL=observations_LL,
diff --git a/R/make_map_info.R b/R/make_map_info.R
index 9d69317..728946f 100644
--- a/R/make_map_info.R
+++ b/R/make_map_info.R
@@ -148,16 +148,17 @@ make_map_info = function( Region, Extrapolation_List, spatial_list=NULL, NN_Extr
   }
 
   #
-  if( fine_scale==TRUE | spatial_list$Method=="Stream_network" ){
-    PlotDF[,'x2i'] = NA
-    PlotDF[ which(Extrapolation_List[["Area_km2_x"]]>0),'x2i'] = 1:length(which(Extrapolation_List[["Area_km2_x"]]>0))
-  }else{
-    PlotDF[,'x2i'] = NN_Extrap$nn.idx[,1]
-  }
-  # Exceptions
-  if( tolower(Region) == "eastern_bering_sea" ){
-    PlotDF = PlotDF[which(PlotDF[,'Lon']<0),]
-  }
+  #if( fine_scale==TRUE | spatial_list$Method=="Stream_network" ){
+  #  PlotDF[,'x2i'] = NA
+  #  PlotDF[ which(Extrapolation_List[["Area_km2_x"]]>0),'x2i'] = 1:length(which(Extrapolation_List[["Area_km2_x"]]>0))
+  #}else{
+  #  PlotDF[,'x2i'] = NN_Extrap$nn.idx[,1]
+  #}
+  ## Exceptions
+  #if( tolower(Region) == "eastern_bering_sea" ){
+  #  PlotDF = PlotDF[which(PlotDF[,'Lon']<0),]
+  #}
+  PlotDF[,'x2i'] = spatial_list$g_e
 
   # Plotting zone has to match zone used for Extrapolation_List, for COG estimates (in UTM via Z_xm) to match map UTM
   if( is.numeric(Extrapolation_List$zone) ){
diff --git a/R/make_settings.R b/R/make_settings.R
index 9940b86..3164482 100644
--- a/R/make_settings.R
+++ b/R/make_settings.R
@@ -25,13 +25,21 @@
 make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE,
   strata.limits=data.frame('STRATA'="All_areas"), zone=NA, FieldConfig, RhoConfig,
   OverdispersionConfig, ObsModel, bias.correct, Options, use_anisotropy,
-  vars_to_correct, Version, treat_nonencounter_as_zero, n_categories, VamConfig ){
+  vars_to_correct, Version, treat_nonencounter_as_zero, n_categories, VamConfig,
+  max_cells, knot_method ){
 
   # Get version
   if(missing(Version)) Version = FishStatsUtils::get_latest_version()
+  purpose_found = FALSE
 
+  ###################
   # Index standardization
+  ###################
+
+  # Deprecated
   if( tolower(purpose) == "index" ){
+    purpose_found = TRUE
+    warning( "The package author recommends using purpose=`index2` for updated defaults; purpose=`index` is retained for backwards compatibility but not recommended" )
     if( convert_version_name(Version) >= convert_version_name("VAST_v7_0_0") ){
       if(missing(FieldConfig)) FieldConfig = matrix( "IID", ncol=2, nrow=3, dimnames=list(c("Omega","Epsilon","Beta"),c("Component_1","Component_2")) )
     }else{
@@ -45,10 +53,35 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE,
     if(missing(treat_nonencounter_as_zero)) treat_nonencounter_as_zero = TRUE
     if(missing(Options)) Options =  c("SD_site_logdensity"=FALSE, "Calculate_Range"=TRUE, "Calculate_effective_area"=TRUE, "treat_nonencounter_as_zero"=treat_nonencounter_as_zero )
     if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" )
+    if(missing(knot_method)) knot_method = "samples"
+    if(missing(max_cells)) max_cells = Inf
+  }
+
+  # Current
+  if( tolower(purpose) == "index2" ){
+    purpose_found = TRUE
+    if( convert_version_name(Version) >= convert_version_name("VAST_v7_0_0") ){
+      if(missing(FieldConfig)) FieldConfig = matrix( "IID", ncol=2, nrow=3, dimnames=list(c("Omega","Epsilon","Beta"),c("Component_1","Component_2")) )
+    }else{
+      if(missing(FieldConfig)) FieldConfig = c("Omega1"="IID", "Epsilon1"="IID", "Omega2"="IID", "Epsilon2"="IID")
+    }
+    if(missing(RhoConfig)) RhoConfig = c("Beta1"=0, "Beta2"=0, "Epsilon1"=0, "Epsilon2"=0)
+    if(missing(VamConfig)) VamConfig = c("Method"=0, "Rank"=0, "Timing"=0)
+    if(missing(OverdispersionConfig)) OverdispersionConfig = c("Eta1"=0, "Eta2"=0)
+    if(missing(ObsModel)) ObsModel = c(2,1)
+    if(missing(bias.correct)) bias.correct = TRUE
+    if(missing(treat_nonencounter_as_zero)) treat_nonencounter_as_zero = TRUE
+    if(missing(Options)) Options =  c("SD_site_logdensity"=FALSE, "Calculate_Range"=TRUE, "Calculate_effective_area"=TRUE, "treat_nonencounter_as_zero"=treat_nonencounter_as_zero )
+    if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" )
+    if(missing(knot_method)) knot_method = "grid"
+    if(missing(max_cells)) max_cells = max( 2000, n_x*10 )
   }
 
+  ###################
   # Condition and density
+  ###################
   if( tolower(purpose) == "condition_and_density" ){
+    purpose_found = TRUE
     if( convert_version_name(Version) >= convert_version_name("VAST_v7_0_0") ){
       if(missing(FieldConfig)) FieldConfig = matrix( c(2,2,"IID",0,0,"IID"), ncol=2, nrow=3, dimnames=list(c("Omega","Epsilon","Beta"),c("Component_1","Component_2")) )
     }else{
@@ -62,10 +95,15 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE,
     if(missing(treat_nonencounter_as_zero)) treat_nonencounter_as_zero = FALSE
     if(missing(Options)) Options =  c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=TRUE, "treat_nonencounter_as_zero"=treat_nonencounter_as_zero )
     if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" )
+    if(missing(knot_method)) knot_method = "samples"
+    if(missing(max_cells)) max_cells = Inf
   }
 
+  ###################
   # Spatial model of intermediate complexity for ecosystems (MICE-in-space)
+  ###################
   if( tolower(purpose) %in% c("mice","interactions") ){
+    purpose_found = TRUE
     if(missing(n_categories)) stop("Must supply `n_categories` when using purpose==`mice`")
     if( convert_version_name(Version) >= convert_version_name("VAST_v7_0_0") ){
       if(missing(FieldConfig)) FieldConfig = matrix( c(n_categories,n_categories,"IID",n_categories,n_categories,"IID"), ncol=2, nrow=3, dimnames=list(c("Omega","Epsilon","Beta"),c("Component_1","Component_2")) )
@@ -79,10 +117,15 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE,
     if(missing(bias.correct)) bias.correct = TRUE
     if(missing(Options)) Options =  c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=FALSE, "Calculate_Fratio"=TRUE, "Estimate_B0"=TRUE )
     if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl", "Bratio_cyl" )
+    if(missing(knot_method)) knot_method = "samples"
+    if(missing(max_cells)) max_cells = Inf
   }
 
+  ###################
   # Spatial model for ordinating species
+  ###################
   if( tolower(purpose) %in% c("ordination") ){
+    purpose_found = TRUE
     if(missing(n_categories)) stop("Must supply `n_categories` when using purpose==`ordination`")
     if( convert_version_name(Version) >= convert_version_name("VAST_v7_0_0") ){
       if(missing(FieldConfig)) FieldConfig = matrix( c(n_categories,n_categories,n_categories,n_categories,n_categories,n_categories), ncol=2, nrow=3, dimnames=list(c("Omega","Epsilon","Beta"),c("Component_1","Component_2")) )
@@ -96,10 +139,19 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE,
     if(missing(bias.correct)) bias.correct = FALSE
     if(missing(Options)) Options =  c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=TRUE, "Project_factors"=TRUE )
     if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" )
+    if(missing(knot_method)) knot_method = "samples"
+    if(missing(max_cells)) max_cells = Inf
   }
 
+  ###################
   # Spatial EOF analysis
+  ###################
+
   if( tolower(purpose) %in% c("eof") ){
+    purpose_found = TRUE
+    if( convert_version_name(Version) >= convert_version_name("VAST_v9_2_0") ){
+      warning( "The package author recommends using purpose=`eof2` for improved model specification; purpose=`eof` is retained for backwards compatibility but not recommended" )
+    }
     if(missing(n_categories)) stop("Must supply `n_categories` when using purpose==`eof`")
     if( convert_version_name(Version) >= convert_version_name("VAST_v7_0_0") ){
       if(missing(FieldConfig)) FieldConfig = matrix( c(0,n_categories,"IID", 0,0,"IID"), ncol=2, nrow=3, dimnames=list(c("Omega","Epsilon","Beta"),c("Component_1","Component_2")) )
@@ -113,10 +165,33 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE,
     if(missing(bias.correct)) bias.correct = FALSE
     if(missing(Options)) Options =  c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=TRUE, "Project_factors"=TRUE )
     if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" )
+    if(missing(knot_method)) knot_method = "samples"
+    if(missing(max_cells)) max_cells = Inf
   }
 
+  if( tolower(purpose) %in% c("eof2") ){
+    purpose_found = TRUE
+    if( convert_version_name(Version) < convert_version_name("VAST_v9_2_0") ){
+      stop("Cannot use purpose=`eof2` with VAST versions < 9.2.0")
+    }
+    if(missing(FieldConfig)) FieldConfig = matrix( c("IID","Identity","IID",2, 0,0,"IID","Identity"), ncol=2, nrow=4, dimnames=list(c("Omega","Epsilon","Beta","Epsilon_year"),c("Component_1","Component_2")) )
+    if(missing(RhoConfig)) RhoConfig = c("Beta1"=0, "Beta2"=3, "Epsilon1"=0, "Epsilon2"=0)
+    if(missing(VamConfig)) VamConfig = c("Method"=0, "Rank"=0, "Timing"=0)
+    if(missing(OverdispersionConfig)) OverdispersionConfig = c("Eta1"=0, "Eta2"=0)
+    if(missing(ObsModel)) ObsModel = c(2,1)
+    if(missing(bias.correct)) bias.correct = FALSE
+    if(missing(Options)) Options =  c("SD_site_logdensity"=FALSE, "Calculate_Range"=FALSE, "Calculate_effective_area"=FALSE, "Calculate_Cov_SE"=TRUE, "Project_factors"=TRUE )
+    if(missing(vars_to_correct)) vars_to_correct = c( "Index_cyl" )
+    if(missing(knot_method)) knot_method = "grid"
+    if(missing(max_cells)) max_cells = max( 2000, n_x*10 )
+  }
+
+  ###################
+  # Other settings and formatting
+  ###################
+
   # Check for bad input
-  if( !( tolower(purpose) %in% c("index","condition_and_density","mice","ordination","eof")) ){
+  if( purpose_found==FALSE ){
     stop("'purpose' is currently set up only for index-standardization models, correlations between condition and density, MICE-in-space models, ordination, or empirical-orthogonal-function analysis")
   }
 
@@ -131,7 +206,7 @@ make_settings = function( n_x, Region, purpose="index", fine_scale=TRUE,
   # Bundle and export
   settings = list("Version"=Version, "n_x"=n_x, "Region"=Region, "strata.limits"=strata.limits, "zone"=zone, "FieldConfig"=FieldConfig, "RhoConfig"=RhoConfig,
     "VamConfig"=VamConfig, "OverdispersionConfig"=OverdispersionConfig, "ObsModel"=ObsModel, "vars_to_correct"=vars_to_correct,
-    "Options"=Options, "grid_size_km"=grid_size_km,
+    "Options"=Options, "grid_size_km"=grid_size_km, "max_cells"=max_cells, "knot_method"=knot_method,
     "Method"=Method, "use_anisotropy"=use_anisotropy, "fine_scale"=fine_scale, "bias.correct"=bias.correct )
   return(settings)
 }
diff --git a/R/make_spatial_info.R b/R/make_spatial_info.R
index 6193f18..29f8a72 100644
--- a/R/make_spatial_info.R
+++ b/R/make_spatial_info.R
@@ -10,7 +10,7 @@
 #' @param n_x the number of vertices in the SPDE mesh (determines the spatial resolution when Method="Mesh")
 #' @param Lon_i Longitude for each sample
 #' @param Lat_i Latitude for each sample
-#' @param knot_method whether to determine location of GMRF vertices based on the location of samples \code{knot_method=`samples`} or extrapolation-grid cells within the specified strata \code{knot_method='grid'}
+#' @param knot_method whether to determine location of GMRF vertices based on the location of samples \code{knot_method=`samples`} or extrapolation-grid cells within the specified strata \code{knot_method='grid'}; default \code{knot_method=NULL} is coerced to \code{knot_method=`samples`}
 #' @param Method a character of either "Grid" or "Mesh" where "Grid" is a 2D AR1 process, and "Mesh" is the SPDE method with geometric anisotropy
 #' @param fine_scale a Boolean indicating whether to ignore (\code{fine_scale=FALSE}) or account for (\code{fine_scale=TRUE}) fine-scale spatial heterogeneity;  See details for more informatino
 #' @param Extrapolation_List the output from \code{Prepare_Extrapolation_Data_Fn}
@@ -34,7 +34,7 @@
 #'
 
 #' @export
-make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method="samples", Method="Mesh",
+make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method=NULL, Method="Mesh",
   grid_size_km=50, grid_size_LL=1, fine_scale=FALSE, Network_sz_LL=NULL,
   iter.max=1000, randomseed=1, nstart=100, DirPath=paste0(getwd(),"/"), Save_Results=FALSE,
   LON_intensity, LAT_intensity, ... ){
@@ -47,6 +47,9 @@ make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method
     stop("Please use fine_scale=FALSE with stream network spatial model; feature fine_scale=TRUE not yet supported for stream network spatial model.")
   }
 
+  # Backwards compatibility for when settings didn't include knot_method, such that settings$knot_method=NULL
+  if(is.null(knot_method)) knot_method = "samples"
+
   # Backwards compatible option for different extrapolation grid
   if( missing(LON_intensity) & missing(LAT_intensity) ){
     if( knot_method=="samples" ){
@@ -178,10 +181,19 @@ make_spatial_info = function( n_x, Lon_i, Lat_i, Extrapolation_List, knot_method
     a_gl = matrix(dist_inp, nrow=n_x)
   }
 
+  # Moving
+  if( fine_scale==TRUE | Method=="Stream_network" ){
+    g_e = rep(NA, length(Extrapolation_List[["Area_km2_x"]]))
+    g_e[ which(Extrapolation_List[["Area_km2_x"]]>0) ] = 1:length(which(Extrapolation_List[["Area_km2_x"]]>0))
+  }else{
+    g_e = PolygonList$NN_Extrap$nn.idx[,1]
+    g_e[ which(Extrapolation_List[["Area_km2_x"]]==0) ] = NA
+  }
+
   # Return
   Return = list( "fine_scale"=fine_scale, "A_is"=A_is, "A_gs"=A_gs, "n_x"=n_x, "n_s"=n_s, "n_g"=nrow(a_gl), "n_i"=nrow(loc_i),
     "MeshList"=MeshList, "GridList"=GridList, "a_gl"=a_gl, "a_xl"=a_gl, "Kmeans"=Kmeans, "knot_i"=knot_i,
-    "loc_i"=as.matrix(loc_i), "loc_x"=as.matrix(loc_x), "loc_g"=as.matrix(loc_g),
+    "loc_i"=as.matrix(loc_i), "loc_x"=as.matrix(loc_x), "loc_g"=as.matrix(loc_g), "g_e"=g_e,
     "Method"=Method, "PolygonList"=PolygonList, "NN_Extrap"=PolygonList$NN_Extrap, "knot_method"=knot_method,
     "latlon_x"=latlon_x, "latlon_g"=latlon_g, "latlon_i"=latlon_i )
   class(Return) = "make_spatial_info"
diff --git a/R/plot_data.R b/R/plot_data.R
index 7a2abc1..873dc28 100644
--- a/R/plot_data.R
+++ b/R/plot_data.R
@@ -18,7 +18,7 @@
 plot_data = function( Extrapolation_List, Spatial_List, Data_Geostat=NULL,
   Lat_i=Data_Geostat[,'Lat'], Lon_i=Data_Geostat[,'Lon'], Year_i=Data_Geostat[,'Year'], PlotDir=paste0(getwd(),"/"),
   Plot1_name="Data_and_knots.png", Plot2_name="Data_by_year.png", col="red", cex=0.01, pch=19,
-  Year_Set, projargs='+proj=longlat', map_resolution="medium", land_color="grey", ...){
+  Year_Set, projargs='+proj=longlat', map_resolution="medium", land_color="grey", country=NULL, ...){
 
   # Check for issues
   if( is.null(Lat_i) | is.null(Lon_i) | is.null(Year_i) ){
@@ -46,35 +46,41 @@ plot_data = function( Extrapolation_List, Spatial_List, Data_Geostat=NULL,
   CRS_proj = sp::CRS( projargs )
 
   # Data for mapping
-  map_data = rnaturalearth::ne_countries(scale=switch(map_resolution, "low"=110, "medium"=50, "high"=10, 50 ))
+  map_data = rnaturalearth::ne_countries(scale=switch(map_resolution, "low"=110, "medium"=50, "high"=10, 50), country=country)
   map_data = sp::spTransform(map_data, CRSobj=CRS_proj)
 
   # Plot data and grid
-  png( file=paste0(PlotDir,Plot1_name), width=6, height=6, res=200, units="in")
-    par( mfrow=c(2,2), mar=c(3,3,2,0), mgp=c(1.75,0.25,0) )
-    which_rows = which( Extrapolation_List[["Area_km2_x"]]>0 & rowSums(Extrapolation_List[["a_el"]])>0 )
-    plot( Extrapolation_List$Data_Extrap[which_rows,c('Lon','Lat')], cex=0.01, main="Extrapolation (Lat-Lon)" )
-    sp::plot( map_data, col=land_color, add=TRUE )
-    if( !any(is.na(Extrapolation_List$Data_Extrap[,c('E_km','N_km')])) ){
-      plot( Extrapolation_List$Data_Extrap[which_rows,c('E_km','N_km')], cex=0.01, main="Extrapolation (North-East)" )
-    }
-    plot( Spatial_List$loc_x, col="red", pch=20, main="Knots (North-East)")
-  dev.off()
+  if( !missing(Extrapolation_List) & !missing(Spatial_List) ){
+    png( file=paste0(PlotDir,Plot1_name), width=6, height=6, res=200, units="in")
+      par( mfrow=c(2,2), mar=c(3,3,2,0), mgp=c(1.75,0.25,0) )
+      which_rows = which( Extrapolation_List[["Area_km2_x"]]>0 & rowSums(Extrapolation_List[["a_el"]])>0 )
+      plot( Extrapolation_List$Data_Extrap[which_rows,c('Lon','Lat')], cex=0.01, main="Extrapolation (Lat-Lon)" )
+      sp::plot( map_data, col=land_color, add=TRUE )
+      if( !any(is.na(Extrapolation_List$Data_Extrap[,c('E_km','N_km')])) ){
+        plot( Extrapolation_List$Data_Extrap[which_rows,c('E_km','N_km')], cex=0.01, main="Extrapolation (North-East)" )
+      }
+      plot( Spatial_List$loc_x, col="red", pch=20, main="Knots (North-East)")
+    dev.off()
+  }
 
   # Plot data by year
   # Use Data_Geostat, instead of TmbData, because I want raw locations, not locations at knots
   if(missing(Year_Set)) Year_Set = min(Year_i):max(Year_i)
+  if( !any(unique(Year_i) %in% Year_Set) ) Year_Set = sort(unique(Year_i))
     Nrow = ceiling( sqrt(length(Year_Set)) )
     Ncol = ceiling( length(Year_Set)/Nrow )
-  png( file=paste0(PlotDir,Plot2_name), width=Ncol*2, height=Nrow*2, res=200, units="in")
+  if(!is.null(Plot2_name)) png( file=paste0(PlotDir,Plot2_name), width=Ncol*2, height=Nrow*2, res=200, units="in")
     par( mfrow=c(Nrow,Ncol), mar=c(0,0,2,0), mgp=c(1.75,0.25,0), oma=c(4,4,0,0) )
     for( t in 1:length(Year_Set) ){
-      Which = which( Year_i == Year_Set[t] )
-      plot( x=Lon_i[Which], y=Lat_i[Which], cex=cex[Which], main=Year_Set[t], xlim=range(Lon_i), ylim=range(Lat_i), xaxt="n", yaxt="n", col=col[Which], pch=pch[Which], ... )
+      plot( 1, type="n", xlim=range(Lon_i), ylim=range(Lat_i), main=Year_Set[t], xaxt="n", yaxt="n" )
       sp::plot( map_data, col=land_color, add=TRUE )
+      Which = which( Year_i == Year_Set[t] )
+      if( length(Which)>0 ){
+        points( x=Lon_i[Which], y=Lat_i[Which], cex=cex[Which], col=col[Which], pch=pch[Which], ... )
+      }
       if( t>(length(Year_Set)-Ncol) ) axis(1)
       if( t%%Ncol == 1 ) axis(2)
       mtext( side=c(1,2), text=c("Longitude","Latitude"), outer=TRUE, line=1)
     }
-  dev.off()
+  if(!is.null(Plot2_name))dev.off()
 }
diff --git a/R/plot_factors.R b/R/plot_factors.R
index 304e4d3..dd7312b 100644
--- a/R/plot_factors.R
+++ b/R/plot_factors.R
@@ -5,17 +5,23 @@
 #'
 #' @inheritParams plot_overdispersion
 #' @inheritParams summarize_covariance
+#' @inheritParams rotate_factors
+#' @inheritParams plot_maps
 #' @param Year_Set plotting-names for time dimension
 #' @param mapdetails_list output from \code{FishStatsUtils::MapDetails_Fn}
 #' @param Dim_year Plotting dimension (row,column) for plot of years (default: square with sufficient size for number of years)
 #' @param Dim_species Plotting dimension (row,column) for plot of categories (default: square with sufficient size for number of categories)
 #' @param plotdir directory for saving plots
 #' @param land_color color for filling in land (use \code{land_color=rgb(0,0,0,alpha=0)} for transparent land)
-#' @param ... additional arguments passed to \code{plot_maps(.)} when plotting spatio-temporal terms Epsilon
+#' @param ... additional arguments passed to \code{plot_maps(.)} and/or \code{plot_variable(.)} when plotting factor-values on a map
 
 #' @export
 plot_factors = function( Report, ParHat, Data, SD=NULL, Year_Set=NULL, category_names=NULL, RotationMethod="PCA",
-  mapdetails_list=NULL, Dim_year=NULL, Dim_species=NULL, plotdir=paste0(getwd(),"/"), land_color="grey", ... ){
+  mapdetails_list=NULL, Dim_year=NULL, Dim_species=NULL, plotdir=paste0(getwd(),"/"), land_color="grey", zlim=NA,
+  testcutoff=1e-4, ... ){
+
+  #
+  if(is.null(mapdetails_list)) message( "`plot_factors(.) skipping plots because argument `mapdetails_list` is missing")
 
   # Extract Options and Options_vec (depends upon version)
   if( all(c("Options","Options_vec") %in% names(Data)) ){
@@ -27,14 +33,21 @@ plot_factors = function( Report, ParHat, Data, SD=NULL, Year_Set=NULL, category_
     Options = Data$Options_list$Options
   }
 
-  # Adds intercept defaults to FieldConfig if missing
-  if( is.vector(Data[["FieldConfig"]]) && length(Data[["FieldConfig"]])==4 ){
-    Data[["FieldConfig"]] = rbind( matrix(Data[["FieldConfig"]],ncol=2,dimnames=list(c("Omega","Epsilon"),c("Component_1","Component_2"))), "Beta"=c("Beta1"=-2,"Beta2"=-2) )
-  }else{
-    if( !is.matrix(Data[["FieldConfig"]]) || !all(dim(Data[["FieldConfig"]])==c(3,2)) ){
-      stop("`FieldConfig` has the wrong dimensions in `Summarize_Covariance`")
-    }
+  #### Deals with backwards compatibility for FieldConfig
+  # Converts from 4-vector to 3-by-2 matrix
+  if( is.vector(Data$FieldConfig) && length(Data$FieldConfig)==4 ){
+    Data$FieldConfig = rbind( matrix(Data$FieldConfig,ncol=2,dimnames=list(c("Omega","Epsilon"),c("Component_1","Component_2"))), "Beta"=c(-2,-2) )
+  }
+  # Converts from 3-by-2 matrix to 4-by-2 matrix
+  if( is.matrix(Data$FieldConfig) & all(dim(Data$FieldConfig)==c(3,2)) ){
+    Data$FieldConfig = rbind( Data$FieldConfig, "Epsilon_time"=c(-3,-3) )
+  }
+  # Checks for errors
+  if( !is.matrix(Data$FieldConfig) || !all(dim(Data$FieldConfig)==c(4,2)) ){
+    stop("`FieldConfig` has the wrong dimensions in `plot_factors`")
   }
+  # Renames
+  dimnames(Data$FieldConfig) = list( c("Omega","Epsilon","Beta","Epsilon_time"), c("Component_1","Component_2") )
 
   # Fill in missing inputs
   if( "D_xct" %in% names(Report) ){
@@ -56,52 +69,72 @@ plot_factors = function( Report, ParHat, Data, SD=NULL, Year_Set=NULL, category_
   Dim_species = Dim(length(category_names))
 
   # Extract loadings matrices (more numerically stable than extracting covariances, and then re-creating Cholesky)
-  Psi2prime_list = Psiprime_list = Lprime_SE_list = Hinv_list = L_SE_list = Lprime_list = L_list = vector("list", length=6)    # Add names at end so that NULL doesn't interfere
+  Psi2prime_list = Psiprime_list = Lprime_SE_list = Hinv_list = L_SE_list = Lprime_list = L_list = vector("list", length=8)    # Add names at end so that NULL doesn't interfere
 
   # Loop through
-  for(i in 1:6){
+  for(i in 1:8){
+
     # Variable names
-    Par_name = c("Omega1", "Epsilon1", "Beta1", "Omega2", "Epsilon2", "Beta2")[i]
-    if(Par_name == "Omega1"){ Var_name = "Omegainput1_sf"; Var2_name = "Omegainput1_gf" }
-    if(Par_name == "Epsilon1"){ Var_name = "Epsiloninput1_sft"; Var2_name = "Epsiloninput1_gft" }
-    if(Par_name == "Beta1"){ Var_name = "beta1_ft"; Var2_name = "missing" }
-    if(Par_name == "Omega2"){ Var_name = "Omegainput2_sf"; Var2_name = "Omegainput2_gf" }
-    if(Par_name == "Epsilon2"){ Var_name = "Epsiloninput2_sft"; Var2_name = "Epsiloninput2_gft" }
-    if(Par_name == "Beta2"){ Var_name = "beta2_ft"; Var2_name = "missing" }
+    Par_name = c("Omega1", "Epsilon1", "Beta1", "EpsilonTime1", "Omega2", "Epsilon2", "Beta2", "EpsilonTime2")[i]
+    Lpar_name = c("L_omega1_z", "L_epsilon1_z", "L_beta1_z", "Ltime_epsilon1_z", "L_omega2_z", "L_epsilon2_z", "L_beta2_z", "Ltime_epsilon2_z")[i]
+
+    # Backwards compatible loading of variables and names
+    if(Par_name == "Omega1"){ Var_name = "Omegainput1_sf"; Var2_name = "Omegainput1_gf"; L_name = "L_omega1_cf" }
+    if(Par_name == "Epsilon1"){ Var_name = "Epsiloninput1_sft"; Var2_name = "Epsiloninput1_gft"; L_name = "L_epsilon1_cf" }
+    if(Par_name == "Beta1"){ Var_name = "beta1_ft"; Var2_name = "missing"; L_name = "L_beta1_cf" }
+    if(Par_name == "EpsilonTime1"){ Var_name = "Epsiloninput1_sff"; Var2_name = "Epsiloninput1_gff"; L_name = "Ltime_epsilon1_tf" }
+    if(Par_name == "Omega2"){ Var_name = "Omegainput2_sf"; Var2_name = "Omegainput2_gf"; L_name = "L_omega2_cf" }
+    if(Par_name == "Epsilon2"){ Var_name = "Epsiloninput2_sft"; Var2_name = "Epsiloninput2_gft"; L_name = "L_epsilon2_cf" }
+    if(Par_name == "Beta2"){ Var_name = "beta2_ft"; Var2_name = "missing"; L_name = "L_beta2_cf" }
+    if(Par_name == "EpsilonTime2"){ Var_name = "Epsiloninput2_sff"; Var2_name = "Epsiloninput2_gff"; L_name = "Ltime_epsilon2_tf" }
 
     # Continue if component is included
     if( as.vector(Data[["FieldConfig"]])[i] > 0 ){
+
       # Get loadings matrix
-      L_list[[i]] = calc_cov( L_z=ParHat[[paste0("L_",tolower(Par_name),"_z")]], n_f=as.vector(Data[["FieldConfig"]])[i], n_c=Data$n_c, returntype="loadings_matrix" )
-      rownames(L_list[[i]]) = category_names
+      if( "L_beta1_cf" %in% names(Report) ){
+        L_list[[i]] = Report[[L_name]]
+      }else{
+        L_list[[i]] = calc_cov( L_z=ParHat[[Lpar_name]], n_f=as.vector(Data[["FieldConfig"]])[i], n_c=Data$n_c, returntype="loadings_matrix" )
+      }
+      if( !Par_name %in% c("EpsilonTime1","EpsilonTime2") ){
+        rownames(L_list[[i]]) = category_names
+      }
 
       # Load SE
       if( class(SD)=="sdreport" ){
-        #L_SE_list[[i]] = calc_cov( L_z=ParHat_SE[[paste0("L_",tolower(Par_name),"_z")]], n_f=as.vector(Data[["FieldConfig"]])[i], n_c=Data$n_c, returntype="loadings_matrix" )
-        #rownames(L_SE_list[[i]]) = category_names
-        rowindex = grep( paste0("L_",tolower(Par_name),"_z"), rownames(SD$cov.fixed) )
-        L_rz = mvtnorm::rmvnorm( n=1e3, mean=ParHat[[paste0("L_",tolower(Par_name),"_z")]], sigma=SD$cov.fixed[rowindex,rowindex] )
-        L_rcf = array(NA, dim=c(nrow(L_rz),dim(L_list[[i]])) )
-        for( rI in 1:nrow(L_rz) ){
-          L_rcf[rI,,] = calc_cov( L_z=L_rz[rI,], n_f=as.vector(Data[["FieldConfig"]])[i], n_c=Data$n_c, returntype="loadings_matrix" )
+        rowindex = grep( Lpar_name, rownames(SD$cov.fixed) )
+        if( length(rowindex)>0 ){
+          L_rz = mvtnorm::rmvnorm( n=1e3, mean=ParHat[[Lpar_name]], sigma=SD$cov.fixed[rowindex,rowindex] )
+          L_rcf = array(NA, dim=c(nrow(L_rz),dim(L_list[[i]])) )
+          for( rI in 1:nrow(L_rz) ){
+            L_rcf[rI,,] = calc_cov( L_z=L_rz[rI,], n_f=as.vector(Data[["FieldConfig"]])[i], n_c=nrow(L_list[[i]]), returntype="loadings_matrix" )
+          }
+          Lmean_cf = apply(L_rcf, MARGIN=2:3, FUN=mean)
+          Lsd_cf = apply(L_rcf, MARGIN=2:3, FUN=sd)
+          L_SE_list[[i]] = Lsd_cf
+          if( nrow(L_SE_list[[i]])==length(category_names) ){
+            rownames(L_SE_list[[i]]) = category_names
+          }
         }
-        Lmean_cf = apply(L_rcf, MARGIN=2:3, FUN=mean)
-        Lsd_cf = apply(L_rcf, MARGIN=2:3, FUN=sd)
-        L_SE_list[[i]] = Lsd_cf
-        rownames(L_SE_list[[i]]) = category_names
       }
 
       # Get covariance
-      Psi_sjt = ParHat[[Var_name]]
+      if(Var_name%in%names(ParHat)) Psi_sjt = ParHat[[Var_name]]
+      if(Var_name%in%names(Report)) Psi_sjt = Report[[Var_name]]
       Psi_gjt = Report[[Var2_name]]
-      ## the betas are transposed compared to others so fix that here
-      if(Var_name %in% c("beta1_ft", "beta2_ft")){
+      ## the betas and EpsilonTimeare transposed compared to others so fix that here
+      if( Par_name %in% c("Beta1","Beta2") ){
         Psi_sjt <- t(Psi_sjt)
       }
+      if( Par_name %in% c("EpsilonTime1","EpsilonTime2") ){
+        Psi_sjt = aperm( Psi_sjt, c(1,3,2) )
+        Psi_gjt = aperm( Psi_gjt, c(1,3,2) )
+      }
       if(is.null(Psi_sjt)){
         stop(paste("Covariance is empty for parameter", Var_name))
       }
-      logkappa = unlist(ParHat[c('logkappa1','logkappa2')])[c(1,1,1,2,2,2)[i]]
+      logkappa = unlist(ParHat[c('logkappa1','logkappa2')])[c(1,1,1,1,2,2,2,2)[i]]
       if(Options_vec[8]==0){
         tau = 1 / (exp(logkappa) * sqrt(4*pi));
       }else if(Options_vec[8]==1){
@@ -109,65 +142,97 @@ plot_factors = function( Report, ParHat, Data, SD=NULL, Year_Set=NULL, category_
       }else stop("Check 'Options_vec[8]' for allowable entries")
 
       # Rotate stuff
-      Var_rot = rotate_factors( L_pj=L_list[[i]], Psi=Psi_sjt/tau, RotationMethod=RotationMethod, testcutoff=1e-4 )
+      Var_rot = rotate_factors( L_pj=L_list[[i]], Psi=Psi_sjt/tau, RotationMethod=RotationMethod, testcutoff=testcutoff, quiet=TRUE )
+      if( Par_name %in% c("EpsilonTime1","EpsilonTime2") ){
+        Var_rot$Psi_rot = aperm( Var_rot$Psi_rot, c(1,3,2) )
+      }
       Report_tmp = list("D_xct"=Var_rot$Psi_rot, "Epsilon1_sct"=Var_rot$Psi_rot, "Epsilon2_sct"=Var_rot$Psi_rot)
       Lprime_list[[i]] = Var_rot$L_pj_rot
-      rownames(Lprime_list[[i]]) = category_names
+      if( !Par_name %in% c("EpsilonTime1","EpsilonTime2") ){
+        rownames(Lprime_list[[i]]) = category_names
+      }
       Psiprime_list[[i]] = Var_rot$Psi_rot
       Hinv_list[[i]] = Var_rot$Hinv
 
       # Extract SEs if available
       if( class(SD)=="sdreport" ){
-        rowindex = grep( paste0("L_",tolower(Par_name),"_z"), rownames(SD$cov.fixed) )
-        L_rz = mvtnorm::rmvnorm( n=1e3, mean=ParHat[[paste0("L_",tolower(Par_name),"_z")]], sigma=SD$cov.fixed[rowindex,rowindex] )
-        Lprime_rcf = array(NA, dim=c(nrow(L_rz),dim(L_list[[i]])) )
-        for( rI in 1:nrow(L_rz) ){
-          tmpmat = calc_cov( L_z=L_rz[rI,], n_f=as.vector(Data[["FieldConfig"]])[i], n_c=Data$n_c, returntype="loadings_matrix" )
-          Lprime_rcf[rI,,] = rotate_factors( L_pj=tmpmat, RotationMethod="PCA", testcutoff=1e-4, quiet=TRUE )$L_pj_rot
+        rowindex = grep( Lpar_name, rownames(SD$cov.fixed) )
+        if( length(rowindex)>0 ){
+          L_rz = mvtnorm::rmvnorm( n=1e3, mean=ParHat[[Lpar_name]], sigma=SD$cov.fixed[rowindex,rowindex] )
+          Lprime_rcf = array(NA, dim=c(nrow(L_rz),dim(L_list[[i]])) )
+          for( rI in 1:nrow(L_rz) ){
+            tmpmat = calc_cov( L_z=L_rz[rI,], n_f=as.vector(Data[["FieldConfig"]])[i], n_c=nrow(L_list[[i]]), returntype="loadings_matrix" )
+            Lprime_rcf[rI,,] = rotate_factors( L_pj=tmpmat, RotationMethod="PCA", testcutoff=testcutoff, quiet=TRUE )$L_pj_rot
+            # Rotate to maximize correlation with original factors, to prevent effect of label switching on SEs
+            for( fI in 1:ncol(Lprime_list[[i]]) ){
+              Lprime_rcf[rI,,fI] = Lprime_rcf[rI,,fI] * sign(cor(Lprime_rcf[rI,,fI],Lprime_list[[i]][,fI]))
+            }
+          }
+          Lmean_cf = apply(Lprime_rcf, MARGIN=2:3, FUN=mean)
+          Lsd_cf = apply(Lprime_rcf, MARGIN=2:3, FUN=sd)
+          Lprime_SE_list[[i]] = Lsd_cf
+          if( nrow(L_SE_list[[i]])==length(category_names) ){
+            rownames(Lprime_SE_list[[i]]) = category_names
+          }
         }
-        Lmean_cf = apply(Lprime_rcf, MARGIN=2:3, FUN=mean)
-        Lsd_cf = apply(Lprime_rcf, MARGIN=2:3, FUN=sd)
-        Lprime_SE_list[[i]] = Lsd_cf
-        rownames(Lprime_SE_list[[i]]) = category_names
       }
 
       # Extract projected factors is available
       if( !is.null(Psi_gjt) ){
-        Var2_rot = rotate_factors( L_pj=L_list[[i]], Psi=Psi_gjt/tau, RotationMethod=RotationMethod, testcutoff=1e-4 )
-        Report_tmp = list("D_xct"=Var2_rot$Psi_rot, "Epsilon1_sct"=Var2_rot$Psi_rot, "Epsilon2_sct"=Var2_rot$Psi_rot)
+        Var2_rot = rotate_factors( L_pj=L_list[[i]], Psi=Psi_gjt/tau, RotationMethod=RotationMethod, testcutoff=testcutoff, quiet=TRUE )
+        if( Par_name %in% c("EpsilonTime1","EpsilonTime2") ){
+          Var2_rot$Psi_rot = aperm( Var2_rot$Psi_rot, c(1,3,2) )
+        }
+        Report2_tmp = list("D_xct"=Var2_rot$Psi_rot, "Epsilon1_sct"=Var2_rot$Psi_rot, "Epsilon2_sct"=Var2_rot$Psi_rot)
         Psi2prime_list[[i]] = Var2_rot$Psi_rot
+      }else{
+        Report2_tmp = NULL
       }
 
       # Plot loadings
       Dim_factor = Dim( as.vector(Data[["FieldConfig"]])[i] )
       png( file=paste0(plotdir,"Factor_loadings--",Par_name,".png"), width=Dim_factor[2]*4, height=Dim_factor[1]*4, units="in", res=200 )
-        par( mfrow=Dim_factor, mar=c(0,2,2,0) )
-        for( cI in 1:as.vector(Data[["FieldConfig"]])[i] ) FishStatsUtils::plot_loadings( L_pj=Var_rot$L_pj_rot, whichfactor=cI )
+        par( mfrow=Dim_factor, mar=c(2,2,1,0), oma=c(0,0,0,0), mgp=c(2,0.5,0), tck=-0.02 )
+        for( cI in 1:as.vector(Data[["FieldConfig"]])[i] ){
+          if( Par_name %in% c("EpsilonTime1","EpsilonTime2") ){
+            plot_loadings( L_pj=Var_rot$L_pj_rot, whichfactor=cI, At=Year_Set, LabelPosition="Side" )
+          }else{
+            plot_loadings( L_pj=Var_rot$L_pj_rot, whichfactor=cI, At=1:nrow(Var_rot$L_pj_rot) )
+          }
+        }
       dev.off()
 
       # Plot factors
-      if( !is.null(mapdetails_list) ){
+      if( !is.null(mapdetails_list) & !is.null(Report2_tmp) ){
 
         # Plot Epsilon
         # Use plot_maps to automatically make one figure per factor
         if( Par_name %in% c("Epsilon1","Epsilon2") ){
-          plot_maps(plot_set=c(6,6,NA,7,7,NA)[i], Report=Report_tmp, PlotDF=mapdetails_list[["PlotDF"]], MapSizeRatio=mapdetails_list[["MapSizeRatio"]],
+          plot_maps(plot_set=c(6,6,NA,6,7,7,NA,7)[i], Report=Report2_tmp, PlotDF=mapdetails_list[["PlotDF"]], MapSizeRatio=mapdetails_list[["MapSizeRatio"]],
             working_dir=plotdir, Year_Set=Year_Set, category_names=paste0("Factor_",1:dim(Var_rot$Psi_rot)[2]),
-            legend_x=mapdetails_list[["Legend"]]$x/100, legend_y=mapdetails_list[["Legend"]]$y/100, ...)
+            legend_x=mapdetails_list[["Legend"]]$x/100, legend_y=mapdetails_list[["Legend"]]$y/100, zlim=zlim, ...)
         }  #
 
         # Plot Omega
         # Use plot_variable to plot all factors on single figure
         if( Par_name %in% c("Omega1", "Omega2")){
-        #  Mat_sf = apply(Report_tmp$D_xct, MARGIN=1:2, FUN=mean)
-        #  plot_maps(plot_set=c(6,6,NA,7,7,NA)[i], Report=Report_tmp, PlotDF=mapdetails_list[["PlotDF"]], MapSizeRatio=mapdetails_list[["MapSizeRatio"]],
-        #    working_dir=plotdir, category_names=paste0("Factor_",1:dim(Var_rot$Psi_rot)[2]),
-        #    legend_x=mapdetails_list[["Legend"]]$x/100, legend_y=mapdetails_list[["Legend"]]$y/100)
-          plot_variable( Y_gt=array(Report_tmp$D_xct[,,1],dim=dim(Report_tmp$D_xct)[1:2]), map_list=mapdetails_list, working_dir=plotdir,
-            panel_labels=paste0("Factor_",1:dim(Var_rot$Psi_rot)[2]), file_name=paste0("Factor_maps--",Par_name) )
+          plot_variable( Y_gt=array(Report_tmp$D_xct[,,1],dim=dim(Report2_tmp$D_xct)[1:2]), map_list=mapdetails_list, working_dir=plotdir,
+            panel_labels=paste0("Factor_",1:dim(Var_rot$Psi_rot)[2]), file_name=paste0("Factor_maps--",Par_name), ... )
         }
 
         ## Doesn't make sense to make maps of beta factors since they aren't spatial
+
+        # Plot EpsilonTime
+        if( Par_name %in% c("EpsilonTime1","EpsilonTime2") ){
+          if( dim(Var_rot$Psi_rot)[2] == length(category_names) ){
+            factor_names = category_names
+          }else{
+            factor_names = paste0("Factor_",1:dim(Var_rot$Psi_rot)[2])
+          }
+          plot_maps(plot_set=c(6,6,NA,6,7,7,NA,7)[i], Report=Report2_tmp, PlotDF=mapdetails_list[["PlotDF"]], MapSizeRatio=mapdetails_list[["MapSizeRatio"]],
+            working_dir=plotdir, category_names=factor_names, Panel="Year",
+            legend_x=mapdetails_list[["Legend"]]$x/100, legend_y=mapdetails_list[["Legend"]]$y/100, zlim=zlim, ...)
+        }  #
       }
     }else{
       Lprime_SE_list[[i]] = L_SE_list[[i]] = L_SE_list[[i]] = Psi2prime_list[[i]] = Psiprime_list[[i]] = Lprime_list[[i]] = L_list[[i]] = "Element not estimated, and therefore empty"
@@ -175,7 +240,7 @@ plot_factors = function( Report, ParHat, Data, SD=NULL, Year_Set=NULL, category_
   }
 
   # Return stuff invisibly
-  names(Hinv_list) = names(Psi2prime_list) = names(Psiprime_list) = names(Lprime_SE_list) = names(L_SE_list) = names(Lprime_list) = names(L_list) = c("Omega1", "Epsilon1", "Beta1", "Omega2", "Epsilon2", "Beta2")
+  names(Hinv_list) = names(Psi2prime_list) = names(Psiprime_list) = names(Lprime_SE_list) = names(L_SE_list) = names(Lprime_list) = names(L_list) = c("Omega1", "Epsilon1", "Beta1", "Epsilon1Time1", "Omega2", "Epsilon2", "Beta2", "Epsilon1Time2")
   Return = list("Loadings"=L_list, "Rotated_loadings"=Lprime_list, "Rotated_factors"=Psiprime_list, "Rotated_projected_factors"=Psi2prime_list, "Rotation_matrices"=Hinv_list)
   if( class(SD)=="sdreport" ){
     Return[["Loadings_SE"]] = L_SE_list
diff --git a/R/plot_loadings.R b/R/plot_loadings.R
index 55b7359..c88ee1f 100644
--- a/R/plot_loadings.R
+++ b/R/plot_loadings.R
@@ -18,25 +18,34 @@
 #' }
 
 #' @export
-plot_loadings = function( L_pj, whichfactor=1, addtitle=TRUE, LabelPosition="Right", Buffer=c(0,0.1), Labels=rownames(L_pj), Cex=1.2,
-  legend_text="Proportion of explained variance", ... ){
+plot_loadings = function( L_pj, At=1:nrow(L_pj), whichfactor=1, addtitle=TRUE, LabelPosition="Right",
+  Buffer=c(0,0.1), Labels=rownames(L_pj), Cex=1.2, legend_text="Proportion of explained variance", ... ){
+
+  # Check for bugs
+  if(length(At) != nrow(L_pj)) stop("Check argument `At` in `plot_loadings(.)`")
 
   # Plotting window
   Ylim = range(L_pj)
   Ylim[1] = ifelse( Ylim[1]>0, 0, Ylim[1] )
   Ylim[2] = ifelse( Ylim[2]<0, 0, Ylim[2] )
   Ylim = Ylim + diff(Ylim)*Buffer
-  plot(1, type="n", xlim=c(0.5,nrow(L_pj)+0.5), ylim=Ylim, xlab="", ylab="", xaxt="n", xaxs="i", ... )
+  Xlim = c(-0.5,0.5) + range(At)
+
+  # Start plot
+  plot(1, type="n", xlim=Xlim, ylim=Ylim, xlab="", ylab="", xaxt="n", xaxs="i", ... )
   if(LabelPosition=="Xaxis") axis( side=1, at=1:nrow(L_pj), labels=TRUE, las=3)
   if(addtitle==TRUE) mtext( text=paste("Factor",whichfactor), side=3, line=0.1, adj=0)
   abline(h=0)
 
   # Loop through categories and plot each
   for(p in 1:nrow(L_pj)){
-    lines(y=c(0,L_pj[p,whichfactor]), x=rep(p,2), lwd=5)
-    if(LabelPosition=="Right") text(x=(p+0.2), y=0+(0.02*max(L_pj[p,whichfactor])), labels=Labels[p], srt=90, pos=4, cex=Cex)
-    if(LabelPosition=="Above") text(x=p, y=max(L_pj), labels=Labels[p], srt=90, pos=3, cex=Cex)
-    if(LabelPosition=="Below") text(x=p, y=min(L_pj), labels=Labels[p], srt=90, pos=1, cex=Cex)
+    lines(y=c(0,L_pj[p,whichfactor]), x=rep(At[p],2), lwd=5)
+    if(LabelPosition=="Right") text(x=(At[p]+0.2), y=0+(0.02*max(L_pj[p,whichfactor])), labels=Labels[p], srt=90, pos=4, cex=Cex)
+    if(LabelPosition=="Above") text(x=At[p], y=max(L_pj), labels=Labels[p], srt=90, pos=3, cex=Cex)
+    if(LabelPosition=="Below") text(x=At[p], y=min(L_pj), labels=Labels[p], srt=90, pos=1, cex=Cex)
+    if(LabelPosition=="Side") axis(1)
+  }
+  if(!is.null(legend_text)){
+    legend( "top", legend=paste0(legend_text," ",round(100*sum(L_pj[,whichfactor]^2)/sum(L_pj^2),1),"%"), bty="n")
   }
-  if(!is.null(legend_text)) legend( "top", legend=paste0(legend_text," ",round(100*sum(L_pj[,whichfactor]^2)/sum(L_pj^2),1),"%"), bty="n")
 }
diff --git a/R/plot_maps.r b/R/plot_maps.r
index 2f383f9..5e6885f 100644
--- a/R/plot_maps.r
+++ b/R/plot_maps.r
@@ -6,6 +6,8 @@
 #'
 #' @inheritParams plot_variable
 #' @inheritParams sample_variable
+#' @inheritParams rnaturalearth::ne_countries
+
 #' @param plot_set integer-vector defining plots to create
 #' \describe{
 #'   \item{plot_set=1}{Probability of encounter/non-encounter}
@@ -43,7 +45,7 @@ plot_maps <-
 function(plot_set=3, Obj=NULL, PlotDF, Sdreport=NULL, projargs='+proj=longlat',
          Panel="Category", Year_Set=NULL, Years2Include=NULL, category_names=NULL, quiet=FALSE,
          working_dir=paste0(getwd(),"/"), MapSizeRatio, n_cells, plot_value="estimate", n_samples=100,
-         Report, TmbData, ...){
+         Report, TmbData, zlim=NULL, country=NULL, ...){
 
   # Local functions
   extract_value = function( Sdreport, Report, Obj, variable_name, plot_value="estimate", n_samples ){
@@ -138,7 +140,9 @@ function(plot_set=3, Obj=NULL, PlotDF, Sdreport=NULL, projargs='+proj=longlat',
 
     # Extract elements
     Array_xct = NULL
-    plot_code <- c("encounter_prob", "pos_catch", "ln_density", "", "", "epsilon_1", "epsilon_2", "linear_predictor_1", "linear_predictor_2", "density_CV", "covariates", "total_density", "covariate_effects_1", "covariate_effects_2", "omega_1", "omega_2")[plot_num]
+    plot_code <- c("encounter_prob", "pos_catch", "ln_density", "", "", "epsilon_1", "epsilon_2",
+      "linear_predictor_1", "linear_predictor_2", "density_CV", "covariates", "total_density",
+      "covariate_effects_1", "covariate_effects_2", "omega_1", "omega_2")[plot_num]
 
     # Extract matrix to plot
     if(plot_num==1){
@@ -275,27 +279,33 @@ function(plot_set=3, Obj=NULL, PlotDF, Sdreport=NULL, projargs='+proj=longlat',
       if("dhat_ktp" %in% names(Report)) stop()
       if("dpred_ktp" %in% names(Report)) stop()
     }
-    #if(plot_num==15){
-    #  # Spatial effects for probability of encounter
-    #  if( quiet==FALSE ) message(" # Plotting spatial effects (Omega) for 1st linear predictor")
-    #  if("D_xt"%in%names(Report)) stop()
-    #  if("D_xct"%in%names(Report)) stop()
-    #  if("D_xcy"%in%names(Report)) Array_xct = Report$Omega1_sc %o% 1
-    #  if("D_gcy"%in%names(Report)) Array_xct = Report$Omega1_gc %o% 1
-    #  if("dhat_ktp" %in% names(Report)) stop()
-    #  if("dpred_ktp" %in% names(Report)) stop()
-    #}
-    #if(plot_num==16){
-    #  # Spatial effects for positive catch rates
-    #  if( quiet==FALSE ) message(" # Plotting spatial effects (Omega) for 2nd linear predictor")
-    #  if("D_xt"%in%names(Report)) stop()
-    #  if("D_xct"%in%names(Report)) stop()
-    #  if("D_xcy"%in%names(Report)) Array_xct = Report$Omega2_sc %o% 1
-    #  if("D_gcy"%in%names(Report)) Array_xct = Report$Omega2_gc %o% 1
-    #  if("dhat_ktp" %in% names(Report)) stop()
-    #  if("dpred_ktp" %in% names(Report)) stop()
-    #}
+    if(plot_num==15){
+      # Spatial effects for probability of encounter
+      if( quiet==FALSE ) message(" # Plotting spatial effects (Omega) for 1st linear predictor")
+      if("D_xt"%in%names(Report)) stop()
+      if("D_xct"%in%names(Report)) stop()
+      if("D_xcy"%in%names(Report)) Array_xct = Report$Omega1_sc %o% 1
+      if("D_gcy"%in%names(Report)) Array_xct = Report$Omega1_gc %o% 1
+      if("dhat_ktp" %in% names(Report)) stop()
+      if("dpred_ktp" %in% names(Report)) stop()
+    }
+    if(plot_num==16){
+      # Spatial effects for positive catch rates
+      if( quiet==FALSE ) message(" # Plotting spatial effects (Omega) for 2nd linear predictor")
+      if("D_xt"%in%names(Report)) stop()
+      if("D_xct"%in%names(Report)) stop()
+      if("D_xcy"%in%names(Report)) Array_xct = Report$Omega2_sc %o% 1
+      if("D_gcy"%in%names(Report)) Array_xct = Report$Omega2_gc %o% 1
+      if("dhat_ktp" %in% names(Report)) stop()
+      if("dpred_ktp" %in% names(Report)) stop()
+    }
     if( is.null(Array_xct)) stop("Problem with `plot_num` in `plot_maps(.)")
+    if( any(abs(Array_xct)==Inf) ) stop("plot_maps(.) has some element of output that is Inf or -Inf, please check results")
+    if( all(Years2Include %in% 1:dim(Array_xct)[3]) ){
+      years_to_include = Years2Include
+    }else{
+      years_to_include = 1:dim(Array_xct)[3]
+    }
 
     # Plot for each category
     if( tolower(Panel)=="category" ){
@@ -305,23 +315,25 @@ function(plot_set=3, Obj=NULL, PlotDF, Sdreport=NULL, projargs='+proj=longlat',
         if(length(dim(Array_xct))==2) Return = Mat_xt = Array_xct
         if(length(dim(Array_xct))==3) Return = Mat_xt = array(as.vector(Array_xct[,cI,]),dim=dim(Array_xct)[c(1,3)])
 
-        file_name = paste0(plot_code, ifelse(Nplot>1, paste0("--",category_names[cI]), "") )
-        plot_args = plot_variable( Y_gt=Mat_xt[,Years2Include,drop=FALSE], map_list=list("PlotDF"=PlotDF, "MapSizeRatio"=MapSizeRatio), projargs=projargs, working_dir=working_dir,
-          panel_labels=Year_Set[Years2Include], file_name=file_name, n_cells=n_cells, ... )
+        file_name = paste0(plot_code, ifelse(Nplot>1, paste0("--",category_names[cI]), ""), ifelse(is.function(plot_value),"-transformed","-predicted") )
+        plot_args = plot_variable( Y_gt=Mat_xt[,years_to_include,drop=FALSE],
+          map_list=list("PlotDF"=PlotDF, "MapSizeRatio"=MapSizeRatio), projargs=projargs, working_dir=working_dir,
+          panel_labels=Year_Set[years_to_include], file_name=file_name, n_cells=n_cells, zlim=zlim, country=country, ... )
       }
     }
     # Plot for each year
     if( tolower(Panel)=="year" ){
-      Nplot = length(Years2Include)
+      Nplot = length(years_to_include)
       for( tI in 1:Nplot){
-        if(length(dim(Array_xct))==2) Mat_xc = Array_xct[,Years2Include[tI],drop=TRUE]
-        if(length(dim(Array_xct))==3) Mat_xc = Array_xct[,,Years2Include[tI],drop=TRUE]
+        if(length(dim(Array_xct))==2) Mat_xc = Array_xct[,years_to_include[tI],drop=TRUE]
+        if(length(dim(Array_xct))==3) Mat_xc = Array_xct[,,years_to_include[tI],drop=TRUE]
         Return = Mat_xc = array( as.vector(Mat_xc), dim=c(dim(Array_xct)[1],Ncategories)) # Reformat to make sure it has same format for everything
 
         # Do plot
-        file_name = paste0(plot_code, ifelse(Nplot>1, paste0("--",Year_Set[Years2Include][tI]), "") )
-        plot_args = plot_variable( Y_gt=Mat_xc, map_list=list("PlotDF"=PlotDF, "MapSizeRatio"=MapSizeRatio), projargs=projargs, working_dir=working_dir,
-          panel_labels=category_names, file_name=file_name, n_cells=n_cells, ... )
+        file_name = paste0(plot_code, ifelse(Nplot>1, paste0("--",Year_Set[years_to_include][tI]), ""), ifelse(is.function(plot_value),"-transformed","-predicted") )
+        plot_args = plot_variable( Y_gt=Mat_xc, map_list=list("PlotDF"=PlotDF, "MapSizeRatio"=MapSizeRatio),
+          projargs=projargs, working_dir=working_dir,
+          panel_labels=category_names, file_name=file_name, n_cells=n_cells, zlim=zlim, country=country, ... )
       }
     }
   }
diff --git a/R/plot_range_edge.R b/R/plot_range_edge.R
index 91f92a9..bbf4642 100644
--- a/R/plot_range_edge.R
+++ b/R/plot_range_edge.R
@@ -11,12 +11,13 @@
 #' @inheritParams sample_variable
 #' @param working_dir Directory for plots
 #' @param quantiles vector
+#' @param calculate_relative_to_average Boolean, whether to calculate edge in UTM coordinates (default), or instead calculate relative to median across all years. The latter reduces standard errors, and is appropriate when checking significance for comparison across years for a single species.  The former (default) is appropriate for checking significance for comparison across species.
 #'
 
 #' @export
 plot_range_edge = function( Sdreport, Obj, Year_Set=NULL, Years2Include=NULL, strata_names=NULL,
   category_names=NULL, working_dir=paste0(getwd(),"/"), quantiles=c(0.05,0.95), n_samples=100,
-  interval_width=1, width=NULL, height=NULL, ...){
+  interval_width=1, width=NULL, height=NULL, calculate_relative_to_average=FALSE, seed=123456, ...){
 
   # Unpack
   Report = Obj$report()
@@ -55,55 +56,41 @@ plot_range_edge = function( Sdreport, Obj, Year_Set=NULL, Years2Include=NULL, st
   }
 
   ##### Local function
-  ## Sample from GMRF using sparse precision
-  #rmvnorm_prec <- function(mu, prec, n.sims) {
-  #  z <- matrix(rnorm(length(mu) * n.sims), ncol=n.sims)
-  #  L <- Matrix::Cholesky(prec, super=TRUE)
-  #  z <- Matrix::solve(L, z, system = "Lt") ## z = Lt^-1 %*% z
-  #  z <- Matrix::solve(L, z, system = "Pt") ## z = Pt    %*% z
-  #  z <- as.matrix(z)
-  #  mu + z
-  #}
-  #
-  ## Sample from densities
-  #u_zr = rmvnorm_prec( mu=Obj$env$last.par.best, prec=Sdreport$jointPrecision, n.sims=n_samples)
-  #D_gcyr = array( NA, dim=c(dim(Report$D_gcy),n_samples) )
-  #for( rI in 1:n_samples ){
-  #  if( rI%%max(1,floor(n_samples/10)) == 1 ) message( "Obtaining sample ", rI, " from predictive distribution for density" )
-  #  D_gcyr[,,,rI] = Obj$report( par=u_zr[,rI] )$D_gcy
-  #}
-  D_gcyr = sample_variable( Sdreport=Sdreport, Obj=Obj, variable_name="D_gcy", n_samples=n_samples )
+  D_gcyr = sample_variable( Sdreport=Sdreport, Obj=Obj, variable_name="D_gcy", n_samples=n_samples, seed=seed )
 
   # Calculate quantiles from observed and sampled densities D_gcy
-  Z_zm = TmbData$Z_gm
-  E_zctm = array(NA, dim=c(length(quantiles),dim(Report$D_gcy)[2:3],ncol(Z_zm)) )
-  E_zctmr = array(NA, dim=c(length(quantiles),dim(Report$D_gcy)[2:3],ncol(Z_zm),n_samples) )
-  prop_zctm = array(NA, dim=c(dim(Report$D_gcy)[1:3],ncol(Z_zm)) )
-  prop_zctmr = array(NA, dim=c(dim(Report$D_gcy)[1:3],ncol(Z_zm),n_samples) )
+  E_zctm = array(NA, dim=c(length(quantiles),dim(Report$D_gcy)[2:3],ncol(TmbData$Z_gm)) )
+  E_zctmr = array(NA, dim=c(length(quantiles),dim(Report$D_gcy)[2:3],ncol(TmbData$Z_gm),n_samples) )
+  Mean_cmr = array(NA, dim=c(dim(Report$D_gcy)[2],ncol(TmbData$Z_gm),n_samples) )
+  prop_zctm = array(NA, dim=c(dim(Report$D_gcy)[1:3],ncol(TmbData$Z_gm)) )
+  prop_zctmr = array(NA, dim=c(dim(Report$D_gcy)[1:3],ncol(TmbData$Z_gm),n_samples) )
   for( rI in 0:n_samples ){
   for( mI in 1:ncol(TmbData$Z_gm) ){
     order_g = order(TmbData$Z_gm[,mI], decreasing=FALSE)
-    Z_zm[,mI] = Z_zm[order_g,mI]
     if(rI==0) prop_zctm[,,,mI] = apply( Report$D_gcy, MARGIN=2:3, FUN=function(vec){cumsum(vec[order_g])/sum(vec)} )
     if(rI>=0) prop_zctmr[,,,mI,rI] = apply( D_gcyr[,,,rI,drop=FALSE], MARGIN=2:3, FUN=function(vec){cumsum(vec[order_g])/sum(vec)} )
 
     # Calculate edge
-    for( zI in 1:dim(E_zctm)[1] ){
     for( cI in 1:dim(E_zctm)[2] ){
-    for( tI in 1:dim(E_zctm)[3] ){
-      if(rI==0){
-        index_tmp = which.min( (prop_zctm[,cI,tI,mI]-quantiles[zI])^2 )
-        #return( TmbData$Z_gm[order_g[index_tmp],mI] )
-        E_zctm[zI,cI,tI,mI] = TmbData$Z_gm[order_g[index_tmp],mI]
-        #return( E_zctm[zI,cI,tI,mI] )
-      }
       if(rI>=1){
-        index_tmp = which.min( (prop_zctmr[,cI,tI,mI,rI]-quantiles[zI])^2 )
-        #return( TmbData$Z_gm[order_g[index_tmp],mI] )
-        E_zctmr[zI,cI,tI,mI,rI] = TmbData$Z_gm[order_g[index_tmp],mI]
-        #return( TmbData$Z_gm[order_g[index_tmp],mI] )
+        if( calculate_relative_to_average==TRUE ){
+          Mean_cmr[cI,mI,rI] = weighted.mean( as.vector(TmbData$Z_gm[,mI]%o%rep(1,dim(Report$D_gcy)[3])), w=as.vector(D_gcyr[,cI,,rI]) )
+        }else{
+          Mean_cmr[cI,mI,rI] = 0
+        }
       }
-    }}}
+      for( zI in 1:dim(E_zctm)[1] ){
+      for( tI in 1:dim(E_zctm)[3] ){
+        if(rI==0){
+          index_tmp = which.min( (prop_zctm[,cI,tI,mI]-quantiles[zI])^2 )
+          E_zctm[zI,cI,tI,mI] = TmbData$Z_gm[order_g[index_tmp],mI]
+        }
+        if(rI>=1){
+          index_tmp = which.min( (prop_zctmr[,cI,tI,mI,rI]-quantiles[zI])^2 )
+          E_zctmr[zI,cI,tI,mI,rI] = TmbData$Z_gm[order_g[index_tmp],mI] - Mean_cmr[cI,mI,rI]
+        }
+      }}
+    }
   }}
   SE_zctm = apply( E_zctmr, MARGIN=1:4, FUN=sd )
   Edge_zctm = abind::abind( "Estimate"=E_zctm, "Std. Error"=SE_zctm, along=5 )
diff --git a/R/plot_results.R b/R/plot_results.R
index de167b8..e307397 100644
--- a/R/plot_results.R
+++ b/R/plot_results.R
@@ -4,6 +4,7 @@
 #' \code{plot_results} plots diagnostics, results, and indices for a given fitted model
 #'
 #' This function takes a fitted VAST model and generates a standard set of diagnostic and visualization plots.
+#' It does this by calling a series of mid-level plotting functions; see list of functions in Value section of documentation.
 #'
 #' @param fit Output from \code{fit_model}
 #' @inheritParams fit_model
@@ -14,6 +15,15 @@
 #' @param ... additional settings to pass to \code{FishStatsUtils::plot_maps}
 #'
 #' @return Invisibly returns a tagged list of outputs generated by standard plots.
+#' \describe{
+#'   \item{Q}{Output from \code{\link{plot_quantile_diagnostic}}; see that documentation for definition of optional arguments and contents}
+#'   \item{Enc_prob}{Output from \code{\link{plot_encounter_diagnostic}}; see that documentation for definition of optional arguments and contents}
+#'   \item{Index}{Output from \code{\link{plot_biomass_index}}; see that documentation for definition of optional arguments and contents}
+#'   \item{Proportions}{Output from \code{\link{calculate_proportion}}; see that documentation for definition of optional arguments and contents}
+#'   \item{Range}{Output from \code{\link{plot_range_index}}; see that documentation for definition of optional arguments and contents}
+#'   \item{Dens_xt}{Output from \code{\link{plot_maps}}; see that documentation for definition of optional arguments and contents}
+#'   \item{Edge}{Output from \code{\link{plot_range_edge}}; see that documentation for definition of optional arguments and contents}
+#' }
 #'
 #' @family wrapper functions
 #' @seealso \code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots
@@ -21,7 +31,8 @@
 #' @export
 plot_results = function( fit, settings=fit$settings, plot_set=3, working_dir=paste0(getwd(),"/"),
   year_labels=fit$year_labels, years_to_plot=fit$years_to_plot, use_biascorr=TRUE, map_list,
-  category_names, check_residuals=TRUE, projargs='+proj=longlat', zrange, n_samples=100, ... ){
+  category_names, check_residuals=TRUE, projargs='+proj=longlat', zrange, n_samples=100,
+  calculate_relative_to_average=FALSE, ... ){
 
   # Check for known issues
   if( is.null(fit$Report)) stop("`fit$Report` is missing, please check inputs")
@@ -36,7 +47,7 @@ plot_results = function( fit, settings=fit$settings, plot_set=3, working_dir=pas
   plot_data_args = combine_lists( "input"=plot_data_args, "default"=list(Extrapolation_List=fit$extrapolation_list,
     Spatial_List=fit$spatial_list, Lat_i=fit$data_frame[,'Lat_i'], Lon_i=fit$data_frame[,'Lon_i'],
     Year_i=fit$data_frame[,'t_i'], PlotDir=working_dir, Year_Set=year_labels), "args_to_use"=formalArgs(plot_data) )
-  Dens_xt = do.call( what=plot_data, args=plot_data_args )
+  do.call( what=plot_data, args=plot_data_args )
 
   # PLot settings
   if( missing(map_list) ){
@@ -53,11 +64,14 @@ plot_results = function( fit, settings=fit$settings, plot_set=3, working_dir=pas
   plot_anisotropy( FileName=paste0(working_dir,"Aniso.png"), Report=fit$Report, TmbData=fit$data_list )
 
   # Plot index
+  plot_biomass_index_args = list(...)
   if( !is.null(fit$parameter_estimates$SD) ){
     message("\n### Making plot of abundance index")
     #if( !all(is.numeric(year_labels)) ) stop("`plot_biomass_index` isn't built to handle non-numeric `year_labels`")
-    Index = plot_biomass_index( DirName=working_dir, TmbData=fit$data_list, Sdreport=fit$parameter_estimates$SD, Year_Set=year_labels,
-      Years2Include=years_to_plot, use_biascorr=use_biascorr, category_names=category_names )
+    plot_biomass_index_args = combine_lists( "input"=plot_biomass_index_args, "default"=list(DirName=working_dir,
+      TmbData=fit$data_list, Sdreport=fit$parameter_estimates$SD, Year_Set=year_labels,
+      Years2Include=years_to_plot, use_biascorr=use_biascorr, category_names=category_names), "args_to_use"=formalArgs(plot_biomass_index) )
+    Index = do.call( what=plot_biomass_index, args=plot_biomass_index_args )
   }else{
     Index = "Not run"
     message("\n### Skipping plot of abundance index; must re-run with standard errors to plot")
@@ -88,14 +102,15 @@ plot_results = function( fit, settings=fit$settings, plot_set=3, working_dir=pas
   }
 
   # Plot range edges
-  if( "jointPrecision"%in%names(fit$parameter_estimates$SD) & n_samples>0 ){
+  if( "jointPrecision"%in%names(fit$parameter_estimates$SD) & n_samples>0 & fit$data_list$Options_list$Options['Calculate_Range']==TRUE ){
     message("\n### Making plot of range edges")
     Edge = plot_range_edge( Obj=fit$tmb_list$Obj, Sdreport=fit$parameter_estimates$SD,
       working_dir=working_dir, Year_Set=year_labels, Years2Include=years_to_plot,
-      category_names=category_names, n_samples=n_samples, quantiles=c(0.05,0.5,0.95) )
+      category_names=category_names, n_samples=n_samples, quantiles=c(0.05,0.5,0.95),
+      calculate_relative_to_average=calculate_relative_to_average )
   }else{
     Edge = "Not run"
-    message("\n### Skipping plot of range edge; only possible if `getJointPrecision=TRUE` and `n_samples`>0")
+    message("\n### Skipping plot of range edge; only possible if `getJointPrecision=TRUE`, `Options['Calculate_Range']=TRUE`, and `n_samples`>0")
   }
 
   # Plot densities
@@ -105,7 +120,7 @@ plot_results = function( fit, settings=fit$settings, plot_set=3, working_dir=pas
   plot_maps_args = combine_lists( "input"=plot_maps_args, "default"=list(plot_set=plot_set, category_names=category_names, TmbData=fit$data_list,
     Report=fit$Report, Sdreport=fit$parameter_estimates$SD, PlotDF=map_list[["PlotDF"]], MapSizeRatio=map_list[["MapSizeRatio"]],
     working_dir=working_dir, Year_Set=year_labels, Years2Include=years_to_plot, legend_x=map_list[["Legend"]]$x/100, legend_y=map_list[["Legend"]]$y/100,
-    Obj=fit$tmb_list$Obj), "args_to_use"=formalArgs(plot_maps) )
+    Obj=fit$tmb_list$Obj, projargs=projargs), "args_to_use"=formalArgs(plot_maps) )
   Dens_xt = do.call( what=plot_maps, args=plot_maps_args )
 
   # Plot quantile-quantile plot
@@ -127,8 +142,8 @@ plot_results = function( fit, settings=fit$settings, plot_set=3, working_dir=pas
   }
 
   # return
-  Return = list( "Q"=Q, "Index"=Index, "Proportions"=Proportions, "Range"=Range, "Dens_xt"=Dens_xt, "Edge"=Edge,
-    "map_list"=map_list, "plot_maps_args"=plot_maps_args )
+  Return = list( "Q"=Q, "Enc_prob"=Enc_prob, "Index"=Index, "Proportions"=Proportions, "Range"=Range, "Dens_xt"=Dens_xt, "Edge"=Edge,
+    "map_list"=map_list, "plot_maps_args"=plot_maps_args, "plot_biomass_index_args"=plot_biomass_index_args )
   return( invisible(Return) )
 }
 
diff --git a/R/plot_variable.R b/R/plot_variable.R
index 43b09e2..da27a47 100644
--- a/R/plot_variable.R
+++ b/R/plot_variable.R
@@ -4,7 +4,7 @@
 #' @description
 #' \code{plot_variable} plots a map and fills in regions with colors to represent intensity in an areal-interpretion of model results
 #'
-#' See \url{https://proj.org/operations/projections/index.html} for a list of projections to pass via \code{projargs}. I often prefer \code{projargs='+proj=natearth +lat_0=0 +units=km'} where argument \code{+lat_0} allows the user to center eastings on a specified latitude. If maps are generating visual artefacts, please try using argument \code{country} to simplify the polygons used to represent land features.
+#' See \url{https://proj.org/operations/projections/index.html} for a list of projections to pass via \code{projargs}. I often prefer \code{projargs='+proj=natearth +lon_0=0 +units=km'} where argument \code{+lon_0} allows the user to center eastings on a specified longitude. If maps are generating visual artefacts, please try using argument \code{country} to simplify the polygons used to represent land features.
 #'
 #' @inheritParams sp::CRS
 #' @inheritParams rnaturalearth::ne_countries
@@ -14,7 +14,7 @@
 #' @param legend_x two numeric values (generally between 0 and 1, but slightly lower/higher values generate colorbars that are just outside the plotting window) giving left and right-hand location of color legend
 #' @param legend_y two numeric values (see legend_y) giving bottom and top location of color legend
 #' @param map_list output from \code{FishStatsUtils::make_map_info}
-#' @param zlim range for defining bounds of color scale
+#' @param zlim range for defining bounds of color scale.  If \code{zlim=NULL}, then a constant scale is inferred from the range of \code{Y_gt} and a color-legend is plotted in the last panel.  If \code{zlim=NA} then a different range is used in each panel from the range of \code{Y_gt[,t]} and a color-legend is plotted in every panel.
 #' @param add boolean indicating whether to add plot to an existing panel figure, or to define a new panel figure
 #' @param outermargintext vector defining text to plot in outer margins of panel figure
 #' @param panel_labels vector defining titles to use for each panel; defaults to blank
@@ -26,7 +26,7 @@
 plot_variable <-
 function( Y_gt, map_list, panel_labels, projargs='+proj=longlat', map_resolution="medium",
          file_name="density", working_dir=paste0(getwd(),"/"), Format="png", Res=200, add=FALSE,
-         outermargintext=c("Eastings","Northings"), zlim, col, mar=c(0,0,2,0), oma=c(4,4,0,0),
+         outermargintext=c("Eastings","Northings"), zlim=NULL, col, mar=c(0,0,2,0), oma=c(4,4,0,0),
          legend_x=c(0,0.05), legend_y=c(0.05,0.45), cex.legend=1, mfrow, land_color="grey",
          n_cells, xlim, ylim, country=NULL, ...){
 
@@ -38,7 +38,7 @@ function( Y_gt, map_list, panel_labels, projargs='+proj=longlat', map_resolution
   if( is.vector(Y_gt)){
     Y_gt = matrix(Y_gt, ncol=1)
   }
-  if( missing(zlim)){
+  if( is.null(zlim)){
     zlim = range(Y_gt, na.rm=TRUE)
   }
   if( missing(map_list) || is.null(map_list$MapSizeRatio) ){
@@ -66,7 +66,7 @@ function( Y_gt, map_list, panel_labels, projargs='+proj=longlat', map_resolution
   if( is.function(col)){
     col = col(1000)
   }
-  if( all(is.numeric(c(legend_x,legend_y))) ){
+  if( !any(is.na(c(legend_x,legend_y))) ){
     if( any(c(legend_x,legend_y) > 1.2) | any(c(legend_x,legend_y) < -0.2) ){
       stop("Check values for `legend_x` and `legend_y`")
     }
@@ -80,7 +80,7 @@ function( Y_gt, map_list, panel_labels, projargs='+proj=longlat', map_resolution
 
   # Data for mapping
   #map_data = rnaturalearth::ne_coastline(scale=switch(map_resolution, "low"=110, "medium"=50, "high"=10, 50 ))# , continent="america")
-  map_data = rnaturalearth::ne_countries(scale=switch(map_resolution, "low"=110, "medium"=50, "high"=10, 50 ), country=country)
+  map_data = rnaturalearth::ne_countries(scale=switch(map_resolution, "low"=110, "medium"=50, "high"=10, 50), country=country)
   map_data = sp::spTransform(map_data, CRSobj=CRS_proj)
 
   ###################
@@ -126,7 +126,9 @@ function( Y_gt, map_list, panel_labels, projargs='+proj=longlat', map_resolution
     Raster_proj = plotKML::vect2rast( Points_proj, cell.size=cell.size )
     if(missing(xlim)) xlim = Raster_proj@bbox[1,]
     if(missing(ylim)) ylim = Raster_proj@bbox[2,]
-    image( Raster_proj, col=col, zlim=zlim, xlim=xlim, ylim=ylim )
+    Zlim = zlim
+    if(is.na(Zlim[1])) Zlim = range(Y_gt[,tI],na.rm=TRUE)
+    image( Raster_proj, col=col, zlim=Zlim, xlim=xlim, ylim=ylim )
 
     # Plot maps using rnaturalearth
     sp::plot( map_data, col=land_color, add=TRUE )
@@ -134,22 +136,22 @@ function( Y_gt, map_list, panel_labels, projargs='+proj=longlat', map_resolution
     # Title and box
     title( panel_labels[tI], line=0.1, cex.main=ifelse(is.null(Par$cex.main), 1.5, Par$cex.main), cex=ifelse(is.null(Par$cex.main), 1.5, Par$cex.main) )
     box()
-  }
 
-  # Include legend
-  if( all(is.numeric(c(legend_x,legend_y))) ){
-    xl = (1-legend_x[1])*par('usr')[1] + (legend_x[1])*par('usr')[2]
-    xr = (1-legend_x[2])*par('usr')[1] + (legend_x[2])*par('usr')[2]
-    yb = (1-legend_y[1])*par('usr')[3] + (legend_y[1])*par('usr')[4]
-    yt = (1-legend_y[2])*par('usr')[3] + (legend_y[2])*par('usr')[4]
-    if( diff(legend_y) > diff(legend_x) ){
-      align = c("lt","rb")[2]
-      gradient = c("x","y")[2]
-    }else{
-      align = c("lt","rb")[1]
-      gradient = c("x","y")[1]
+    # Include legend
+    if( !any(is.na(c(legend_x,legend_y))) & (tI==ncol(Y_gt) | is.na(zlim[1])) ){
+      xl = (1-legend_x[1])*par('usr')[1] + (legend_x[1])*par('usr')[2]
+      xr = (1-legend_x[2])*par('usr')[1] + (legend_x[2])*par('usr')[2]
+      yb = (1-legend_y[1])*par('usr')[3] + (legend_y[1])*par('usr')[4]
+      yt = (1-legend_y[2])*par('usr')[3] + (legend_y[2])*par('usr')[4]
+      if( diff(legend_y) > diff(legend_x) ){
+        align = c("lt","rb")[2]
+        gradient = c("x","y")[2]
+      }else{
+        align = c("lt","rb")[1]
+        gradient = c("x","y")[1]
+      }
+      plotrix::color.legend(xl=xl, yb=yb, xr=xr, yt=yt, legend=round(seq(Zlim[1],Zlim[2],length=4),1), rect.col=col, cex=cex.legend, align=align, gradient=gradient)
     }
-    plotrix::color.legend(xl=xl, yb=yb, xr=xr, yt=yt, legend=round(seq(zlim[1],zlim[2],length=4),1), rect.col=col, cex=cex.legend, align=align, gradient=gradient)
   }
 
   # Margin text
diff --git a/R/rotate_factors.R b/R/rotate_factors.R
index 764aeeb..04f05c8 100644
--- a/R/rotate_factors.R
+++ b/R/rotate_factors.R
@@ -5,7 +5,7 @@
 #'
 #' @param Cov_jj Covariance calculated from loadings matrix
 #' @param L_pj Loadings matrix for `p` categories and `j` factors (calculated from \code{Cov_jj} if it is provided)
-#' @param Psi Array of factors (1st dimension: spatial knots;  2nd dimension: factors;  3rd dimension:  time)
+#' @param Psi_sjt Array of factors (1st dimension: spatial knots;  2nd dimension: factors;  3rd dimension:  time)
 #' @param RotationMethod Method used for rotation, Options: "PCA" (recommended) or "Varimax"
 #' @param testcutoff tolerance for numerical rounding when confirming that rotation doesn't effect results
 
@@ -18,18 +18,18 @@
 #' }
 
 #' @export
-rotate_factors = function( Cov_jj=NULL, L_pj=NULL, Psi=NULL, RotationMethod="PCA", testcutoff=1e-10,
+rotate_factors = function( Cov_jj=NULL, L_pj=NULL, Psi_sjt=NULL, RotationMethod="PCA", testcutoff=1e-10,
   quiet=FALSE ){
 
   # If missing time, add a third dimension
-  if( length(dim(Psi))==2 ){
-    Psi = array( Psi, dim=c(dim(Psi),1) )
+  if( length(dim(Psi_sjt))==2 ){
+    Psi_sjt = array( Psi_sjt, dim=c(dim(Psi_sjt),1) )
   }
 
   # Local functions
   approx_equal = function(m1,m2,denominator=mean(m1+m2),d=1e-10) (2*abs(m1-m2)/denominator) < d
   trunc_machineprec = function(n) ifelse(n<1e-10,0,n)
-  Nknots = dim(Psi)[1]
+  Nknots = dim(Psi_sjt)[1]
   Nfactors = ncol(L_pj)
   Nyears = nrow(L_pj)
 
@@ -53,11 +53,11 @@ rotate_factors = function( Cov_jj=NULL, L_pj=NULL, Psi=NULL, RotationMethod="PCA
     Hinv = varimax( L_pj, normalize=FALSE )
     H = solve(Hinv$rotmat)
     L_pj_rot = L_pj %*% Hinv$rotmat
-    if( !is.null(Psi) ){
-      Psi_rot = array(NA, dim=dim(Psi))
+    if( !is.null(Psi_sjt) ){
+      Psi_rot = array(NA, dim=dim(Psi_sjt))
       # My new factors
       for( n in 1:Nknots ){
-        Psi_rot[n,,] = H %*% Psi[n,,]
+        Psi_rot[n,,] = H %*% Psi_sjt[n,,]
       }
     }
   }
@@ -76,18 +76,18 @@ rotate_factors = function( Cov_jj=NULL, L_pj=NULL, Psi=NULL, RotationMethod="PCA
     rownames(L_pj_rot) = rownames(L_pj)
     # My new factors
     H = corpcor::pseudoinverse(L_pj_rot) %*% L_pj
-    Hinv = list("rotmat"=solve(H))
-    if( !is.null(Psi) ){
-      Psi_rot = array(NA, dim=dim(Psi))
+    Hinv = list("rotmat"=corpcor::pseudoinverse(H))
+    if( !is.null(Psi_sjt) ){
+      Psi_rot = array(NA, dim=dim(Psi_sjt))
       for( n in 1:Nknots ){
-        Psi_rot[n,,] = H %*% Psi[n,,]
+        Psi_rot[n,,] = H %*% Psi_sjt[n,,]
       }
     }
   }
 
   # Flip around
   for( j in 1:dim(L_pj_rot)[2] ){
-    if( !is.null(Psi) ){
+    if( !is.null(Psi_sjt) ){
       Psi_rot[,j,] = Psi_rot[,j,] * sign(sum(L_pj_rot[,j]))
     }
     L_pj_rot[,j] = L_pj_rot[,j] * sign(sum(L_pj_rot[,j]))
@@ -97,14 +97,18 @@ rotate_factors = function( Cov_jj=NULL, L_pj=NULL, Psi=NULL, RotationMethod="PCA
   # Check covariance matrix
     # Should be identical for rotated and unrotated
   if( !is.na(testcutoff) ){
-    if( !all(approx_equal(L_pj%*%t(L_pj),L_pj_rot%*%t(L_pj_rot), d=testcutoff)) ) stop("Covariance matrix is changed by rotation")
+    if( !all(approx_equal(L_pj%*%t(L_pj),L_pj_rot%*%t(L_pj_rot), d=testcutoff)) ){
+      stop("Covariance matrix is changed by rotation")
+    }
     # Check linear predictor
       # Should give identical predictions as unrotated
-    if( !is.null(Psi) ){
-      for(i in 1:dim(Psi)[[1]]){
-      for(j in 1:dim(Psi)[[3]]){
-        MaxDiff = max(L_pj%*%Psi[i,,j] - L_pj_rot%*%Psi_rot[i,,j])
-        if( !all(approx_equal(L_pj%*%Psi[i,,j],L_pj_rot%*%Psi_rot[i,,j], d=testcutoff, denominator=1)) ) stop(paste0("Linear predictor is wrong for site ",i," and time ",j," with difference ",MaxDiff))
+    if( !is.null(Psi_sjt) ){
+      for(i in 1:dim(Psi_sjt)[[1]]){
+      for(j in 1:dim(Psi_sjt)[[3]]){
+        MaxDiff = max(L_pj%*%Psi_sjt[i,,j] - L_pj_rot%*%Psi_rot[i,,j])
+        if( !all(approx_equal(L_pj%*%Psi_sjt[i,,j],L_pj_rot%*%Psi_rot[i,,j], d=testcutoff, denominator=1)) ){
+          stop(paste0("Linear predictor is wrong for site ",i," and time ",j," with difference ",MaxDiff))
+        }
       }}
     }
     # Check rotation matrix
@@ -112,13 +116,15 @@ rotate_factors = function( Cov_jj=NULL, L_pj=NULL, Psi=NULL, RotationMethod="PCA
       # Doesn't have det(R) = 1; determinant(Hinv$rotmat)!=1 ||
     Diag = Hinv$rotmat %*% t(Hinv$rotmat)
     diag(Diag) = ifelse( diag(Diag)==0,1,diag(Diag) )
-    if( !all(approx_equal(Diag,diag(Nfactors), d=testcutoff)) ) stop("Rotation matrix is not a rotation")
+    if( !all(approx_equal(Diag,diag(Nfactors), d=testcutoff)) ){
+      stop("Rotation matrix is not a rotation")
+    }
   }
 
   # Return stuff
   Return = list( "L_pj_rot"=L_pj_rot, "Hinv"=Hinv, "L_pj"=L_pj )
   if(RotationMethod=="PCA") Return[["Eigen"]] = Eigen
-  if( !is.null(Psi) ){
+  if( !is.null(Psi_sjt) ){
     Return[["Psi_rot"]] = Psi_rot
   }
   return( Return )
diff --git a/R/sample_variable.R b/R/sample_variable.R
index ab3f69e..5e90a31 100644
--- a/R/sample_variable.R
+++ b/R/sample_variable.R
@@ -11,10 +11,11 @@
 #' @param Obj Fitted TMB object from package `VAST`, i.e., output from `\code{fit_model(...)$tmb_list$Obj}`
 #' @param variable_name name of variable available in report using \code{Obj$report()}
 #' @param n_samples number of samples from the joint predictive distribution for fixed and random effects.  Default is 100, which is slow.
+#' @param seed integer used to set random-number seed when sampling variables, as passed to \code{set.seed(.)}
 #'
 
 #' @export
-sample_variable = function( Sdreport, Obj, variable_name, n_samples=100 ){
+sample_variable = function( Sdreport, Obj, variable_name, n_samples=100, seed=123456 ){
 
   # Informative error messages
   if( !("jointPrecision" %in% names(Sdreport)) ){
@@ -26,17 +27,18 @@ sample_variable = function( Sdreport, Obj, variable_name, n_samples=100 ){
 
   #### Local function
   # Sample from GMRF using sparse precision
-  rmvnorm_prec <- function(mu, prec, n.sims) {
+  rmvnorm_prec <- function(mu, prec, n.sims, seed) {
+    set.seed(seed)
     z <- matrix(rnorm(length(mu) * n.sims), ncol=n.sims)
     L <- Matrix::Cholesky(prec, super=TRUE)
     z <- Matrix::solve(L, z, system = "Lt") ## z = Lt^-1 %*% z
     z <- Matrix::solve(L, z, system = "Pt") ## z = Pt    %*% z
     z <- as.matrix(z)
-    mu + z
+    return(mu + z)
   }
 
   # Sample from joint distribution
-  u_zr = rmvnorm_prec( mu=Obj$env$last.par.best, prec=Sdreport$jointPrecision, n.sims=n_samples)
+  u_zr = rmvnorm_prec( mu=Obj$env$last.par.best, prec=Sdreport$jointPrecision, n.sims=n_samples, seed=seed)
 
   # Extract variable for each sample
   message( "# Obtaining samples from predictive distribution for variable ", variable_name )
diff --git a/R/simulate_data.R b/R/simulate_data.R
new file mode 100644
index 0000000..f7cc047
--- /dev/null
+++ b/R/simulate_data.R
@@ -0,0 +1,71 @@
+
+#' @title
+#' Simulate new dynamics and sampling data
+#'
+#' @description
+#' \code{simulate_data} conducts a parametric bootstrap to simulate new data and potentially simulate new population dynamics and associated variables
+#'
+#' @param fit output form \code{fit_model(.)}
+#' @param type integer stating what type of simulation to use.  \code{type=1} simulates new data conditional upon estimated fixed and random effects.  \code{type=2} simulates new random effects conditional upon fixed effects, and new data conditional upon both.  \code{type=3} simulates new fixed and random effects from the joint precision matrix, and new data conditional upon these values.
+#' @param random_seed integer passed to \code{\link[base]{set.seed}} whenever \code{type=3}, where the default value \code{random_seed=NULL} resets the random-number seed. Argument no effect when \code{type!=3} because TMB has no interface for setting the random-number seed in C++
+#'
+
+#' @return Report object containing new data and population variables including
+#' \describe{
+#'   \item{b_i}{New simulated data}
+#'   \item{D_gcy}{Density for each grid cell g, category c, and year y}
+#'   \item{Index_cyl}{Index of abundance for each category c, year y, and stratum l}
+#' }
+
+#' @export
+simulate_data = function( fit, type=1, random_seed=NULL ){
+
+  # Sample from GMRF using sparse precision
+  rmvnorm_prec <- function(mu, prec, n.sims, random_seed ) {
+    set.seed( random_seed )
+    z = matrix(rnorm(length(mu) * n.sims), ncol=n.sims)
+    L = Matrix::Cholesky(prec, super=TRUE)
+    z = Matrix::solve(L, z, system = "Lt") ## z = Lt^-1 %*% z
+    z = Matrix::solve(L, z, system = "Pt") ## z = Pt    %*% z
+    z = as.matrix(z)
+    return(mu + z)
+  }
+
+  # Extract stuff
+  Obj = fit$tmb_list$Obj
+  if( !is.null(random_seed) & type!=3 ){
+    stop("Specifying argument `random_seed` in `simulate_data(.) only works when `type=3`")
+  }
+
+  # Simulate conditional upon fixed and random effect estimates
+  if( type==1 ){
+    Obj$env$data$Options_list$Options['simulate_random_effects'] = FALSE
+    Return = Obj$simulate( complete=TRUE )
+  }
+
+  # Simulate new random effects and data
+  if( type==2 ){
+    Obj$env$data$Options_list$Options['simulate_random_effects'] = TRUE
+    Return = Obj$simulate( complete=TRUE )
+  }
+
+  # Simulate from predictive distribution, and then new data
+  if( type==3 ){
+    # Informative error messages
+    if( !("jointPrecision" %in% names(fit$parameter_estimates$SD)) ){
+      stop("jointPrecision not present in fit$parameter_estimates$SD; please re-run with `getJointPrecision=TRUE`")
+    }
+
+    # Sample from joint distribution
+    newpar = rmvnorm_prec( mu=Obj$env$last.par.best, prec=fit$parameter_estimates$SD$jointPrecision, n.sims=1, random_seed=random_seed )[,1]
+
+    # Simulate
+    Obj$env$data$Options_list$Options['simulate_random_effects'] = FALSE
+    Return = Obj$simulate( par=newpar, complete=TRUE )
+  }
+
+  # Return
+  return( Return )
+}
+
+
diff --git a/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.dbf b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.dbf
new file mode 100644
index 0000000..8fc04a5
Binary files /dev/null and b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.prj b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.shp b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.shp
new file mode 100644
index 0000000..19f5543
Binary files /dev/null and b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.shp differ
diff --git a/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.shx b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.shx
new file mode 100644
index 0000000..b0d7c32
Binary files /dev/null and b/inst/region_shapefiles/ATL-IBTS-Q1/Shapefile.shx differ
diff --git a/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.dbf b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.dbf
new file mode 100644
index 0000000..59b81e0
Binary files /dev/null and b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.prj b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.shp b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.shp
new file mode 100644
index 0000000..be0c590
Binary files /dev/null and b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.shp differ
diff --git a/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.shx b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.shx
new file mode 100644
index 0000000..6e81b73
Binary files /dev/null and b/inst/region_shapefiles/ATL-IBTS-Q4/Shapefile.shx differ
diff --git a/inst/region_shapefiles/BITS/Shapefile.dbf b/inst/region_shapefiles/BITS/Shapefile.dbf
new file mode 100644
index 0000000..ba67cdc
Binary files /dev/null and b/inst/region_shapefiles/BITS/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/BITS/Shapefile.prj b/inst/region_shapefiles/BITS/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/BITS/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/BITS/Shapefile.shp b/inst/region_shapefiles/BITS/Shapefile.shp
new file mode 100644
index 0000000..c8a6553
Binary files /dev/null and b/inst/region_shapefiles/BITS/Shapefile.shp differ
diff --git a/inst/region_shapefiles/BITS/Shapefile.shx b/inst/region_shapefiles/BITS/Shapefile.shx
new file mode 100644
index 0000000..8e466cb
Binary files /dev/null and b/inst/region_shapefiles/BITS/Shapefile.shx differ
diff --git a/inst/region_shapefiles/BTS-VIIa/Shapefile.dbf b/inst/region_shapefiles/BTS-VIIa/Shapefile.dbf
new file mode 100644
index 0000000..a57d3f7
Binary files /dev/null and b/inst/region_shapefiles/BTS-VIIa/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/BTS-VIIa/Shapefile.prj b/inst/region_shapefiles/BTS-VIIa/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/BTS-VIIa/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/BTS-VIIa/Shapefile.shp b/inst/region_shapefiles/BTS-VIIa/Shapefile.shp
new file mode 100644
index 0000000..84f36b6
Binary files /dev/null and b/inst/region_shapefiles/BTS-VIIa/Shapefile.shp differ
diff --git a/inst/region_shapefiles/BTS-VIIa/Shapefile.shx b/inst/region_shapefiles/BTS-VIIa/Shapefile.shx
new file mode 100644
index 0000000..37c1608
Binary files /dev/null and b/inst/region_shapefiles/BTS-VIIa/Shapefile.shx differ
diff --git a/inst/region_shapefiles/BTS/Shapefile.dbf b/inst/region_shapefiles/BTS/Shapefile.dbf
new file mode 100644
index 0000000..2cbbcc2
Binary files /dev/null and b/inst/region_shapefiles/BTS/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/BTS/Shapefile.prj b/inst/region_shapefiles/BTS/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/BTS/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/BTS/Shapefile.shp b/inst/region_shapefiles/BTS/Shapefile.shp
new file mode 100644
index 0000000..6609ca6
Binary files /dev/null and b/inst/region_shapefiles/BTS/Shapefile.shp differ
diff --git a/inst/region_shapefiles/BTS/Shapefile.shx b/inst/region_shapefiles/BTS/Shapefile.shx
new file mode 100644
index 0000000..e39e364
Binary files /dev/null and b/inst/region_shapefiles/BTS/Shapefile.shx differ
diff --git a/inst/region_shapefiles/EVHOE/Shapefile.dbf b/inst/region_shapefiles/EVHOE/Shapefile.dbf
new file mode 100644
index 0000000..0904381
Binary files /dev/null and b/inst/region_shapefiles/EVHOE/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/EVHOE/Shapefile.prj b/inst/region_shapefiles/EVHOE/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/EVHOE/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/EVHOE/Shapefile.shp b/inst/region_shapefiles/EVHOE/Shapefile.shp
new file mode 100644
index 0000000..68f2943
Binary files /dev/null and b/inst/region_shapefiles/EVHOE/Shapefile.shp differ
diff --git a/inst/region_shapefiles/EVHOE/Shapefile.shx b/inst/region_shapefiles/EVHOE/Shapefile.shx
new file mode 100644
index 0000000..8923a43
Binary files /dev/null and b/inst/region_shapefiles/EVHOE/Shapefile.shx differ
diff --git a/inst/region_shapefiles/IE-IGFS/Shapefile.dbf b/inst/region_shapefiles/IE-IGFS/Shapefile.dbf
new file mode 100644
index 0000000..5789747
Binary files /dev/null and b/inst/region_shapefiles/IE-IGFS/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/IE-IGFS/Shapefile.prj b/inst/region_shapefiles/IE-IGFS/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/IE-IGFS/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/IE-IGFS/Shapefile.shp b/inst/region_shapefiles/IE-IGFS/Shapefile.shp
new file mode 100644
index 0000000..c44b313
Binary files /dev/null and b/inst/region_shapefiles/IE-IGFS/Shapefile.shp differ
diff --git a/inst/region_shapefiles/IE-IGFS/Shapefile.shx b/inst/region_shapefiles/IE-IGFS/Shapefile.shx
new file mode 100644
index 0000000..7697cef
Binary files /dev/null and b/inst/region_shapefiles/IE-IGFS/Shapefile.shx differ
diff --git a/inst/region_shapefiles/NIGFS/Shapefile.dbf b/inst/region_shapefiles/NIGFS/Shapefile.dbf
new file mode 100644
index 0000000..888ee24
Binary files /dev/null and b/inst/region_shapefiles/NIGFS/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/NIGFS/Shapefile.prj b/inst/region_shapefiles/NIGFS/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/NIGFS/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/NIGFS/Shapefile.shp b/inst/region_shapefiles/NIGFS/Shapefile.shp
new file mode 100644
index 0000000..4b93dff
Binary files /dev/null and b/inst/region_shapefiles/NIGFS/Shapefile.shp differ
diff --git a/inst/region_shapefiles/NIGFS/Shapefile.shx b/inst/region_shapefiles/NIGFS/Shapefile.shx
new file mode 100644
index 0000000..6d5f346
Binary files /dev/null and b/inst/region_shapefiles/NIGFS/Shapefile.shx differ
diff --git a/inst/region_shapefiles/NS_IBTS/Shapefile.dbf b/inst/region_shapefiles/NS_IBTS/Shapefile.dbf
new file mode 100644
index 0000000..7205f6d
Binary files /dev/null and b/inst/region_shapefiles/NS_IBTS/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/NS_IBTS/Shapefile.prj b/inst/region_shapefiles/NS_IBTS/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/NS_IBTS/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/NS_IBTS/Shapefile.shp b/inst/region_shapefiles/NS_IBTS/Shapefile.shp
new file mode 100644
index 0000000..6e1b1cb
Binary files /dev/null and b/inst/region_shapefiles/NS_IBTS/Shapefile.shp differ
diff --git a/inst/region_shapefiles/NS_IBTS/Shapefile.shx b/inst/region_shapefiles/NS_IBTS/Shapefile.shx
new file mode 100644
index 0000000..d8b4e5e
Binary files /dev/null and b/inst/region_shapefiles/NS_IBTS/Shapefile.shx differ
diff --git a/inst/region_shapefiles/PT-IBTS/Shapefile.dbf b/inst/region_shapefiles/PT-IBTS/Shapefile.dbf
new file mode 100644
index 0000000..f8ed9bf
Binary files /dev/null and b/inst/region_shapefiles/PT-IBTS/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/PT-IBTS/Shapefile.prj b/inst/region_shapefiles/PT-IBTS/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/PT-IBTS/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/PT-IBTS/Shapefile.shp b/inst/region_shapefiles/PT-IBTS/Shapefile.shp
new file mode 100644
index 0000000..bba8535
Binary files /dev/null and b/inst/region_shapefiles/PT-IBTS/Shapefile.shp differ
diff --git a/inst/region_shapefiles/PT-IBTS/Shapefile.shx b/inst/region_shapefiles/PT-IBTS/Shapefile.shx
new file mode 100644
index 0000000..0f66369
Binary files /dev/null and b/inst/region_shapefiles/PT-IBTS/Shapefile.shx differ
diff --git a/inst/region_shapefiles/SP-ARSA/Shapefile.dbf b/inst/region_shapefiles/SP-ARSA/Shapefile.dbf
new file mode 100644
index 0000000..16a19c3
Binary files /dev/null and b/inst/region_shapefiles/SP-ARSA/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/SP-ARSA/Shapefile.prj b/inst/region_shapefiles/SP-ARSA/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/SP-ARSA/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/SP-ARSA/Shapefile.shp b/inst/region_shapefiles/SP-ARSA/Shapefile.shp
new file mode 100644
index 0000000..9352361
Binary files /dev/null and b/inst/region_shapefiles/SP-ARSA/Shapefile.shp differ
diff --git a/inst/region_shapefiles/SP-ARSA/Shapefile.shx b/inst/region_shapefiles/SP-ARSA/Shapefile.shx
new file mode 100644
index 0000000..e2be37c
Binary files /dev/null and b/inst/region_shapefiles/SP-ARSA/Shapefile.shx differ
diff --git a/inst/region_shapefiles/SP-NORTH/Shapefile.dbf b/inst/region_shapefiles/SP-NORTH/Shapefile.dbf
new file mode 100644
index 0000000..736dd43
Binary files /dev/null and b/inst/region_shapefiles/SP-NORTH/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/SP-NORTH/Shapefile.prj b/inst/region_shapefiles/SP-NORTH/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/SP-NORTH/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/SP-NORTH/Shapefile.shp b/inst/region_shapefiles/SP-NORTH/Shapefile.shp
new file mode 100644
index 0000000..b1392d9
Binary files /dev/null and b/inst/region_shapefiles/SP-NORTH/Shapefile.shp differ
diff --git a/inst/region_shapefiles/SP-NORTH/Shapefile.shx b/inst/region_shapefiles/SP-NORTH/Shapefile.shx
new file mode 100644
index 0000000..701988a
Binary files /dev/null and b/inst/region_shapefiles/SP-NORTH/Shapefile.shx differ
diff --git a/inst/region_shapefiles/SP-PORC/Shapefile.dbf b/inst/region_shapefiles/SP-PORC/Shapefile.dbf
new file mode 100644
index 0000000..65730be
Binary files /dev/null and b/inst/region_shapefiles/SP-PORC/Shapefile.dbf differ
diff --git a/inst/region_shapefiles/SP-PORC/Shapefile.prj b/inst/region_shapefiles/SP-PORC/Shapefile.prj
new file mode 100644
index 0000000..6e3eaed
--- /dev/null
+++ b/inst/region_shapefiles/SP-PORC/Shapefile.prj
@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/inst/region_shapefiles/SP-PORC/Shapefile.shp b/inst/region_shapefiles/SP-PORC/Shapefile.shp
new file mode 100644
index 0000000..d9304b4
Binary files /dev/null and b/inst/region_shapefiles/SP-PORC/Shapefile.shp differ
diff --git a/inst/region_shapefiles/SP-PORC/Shapefile.shx b/inst/region_shapefiles/SP-PORC/Shapefile.shx
new file mode 100644
index 0000000..42ea020
Binary files /dev/null and b/inst/region_shapefiles/SP-PORC/Shapefile.shx differ
diff --git a/man/Geostat_Sim.Rd b/man/Geostat_Sim.Rd
index 7183b63..e6c2fcb 100644
--- a/man/Geostat_Sim.Rd
+++ b/man/Geostat_Sim.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/Geostat_Sim.R
+% Please edit documentation in R/deprecated.R
 \name{Geostat_Sim}
 \alias{Geostat_Sim}
 \title{Simulate data}
diff --git a/man/PlotMap_Fn.Rd b/man/PlotMap_Fn.Rd
index cc543ad..7f391c2 100644
--- a/man/PlotMap_Fn.Rd
+++ b/man/PlotMap_Fn.Rd
@@ -18,6 +18,8 @@ PlotMap_Fn(MappingDetails, Mat, PlotDF, MapSizeRatio = c(`Width(in)` = 4,
 
 \item{Year_Set}{Year names for labeling panels}
 
+\item{zlim}{range for defining bounds of color scale.  If \code{zlim=NULL}, then a constant scale is inferred from the range of \code{Y_gt} and a color-legend is plotted in the last panel.  If \code{zlim=NA} then a different range is used in each panel from the range of \code{Y_gt[,t]} and a color-legend is plotted in every panel.}
+
 \item{plot_legend_fig}{Boolean, whether to plot a separate figure for the heatmap legend or not}
 
 \item{land_color}{color for filling in land (use \code{land_color=rgb(0,0,0,alpha=0)} for transparent land)}
diff --git a/man/Rerun_Fn.Rd b/man/Rerun_Fn.Rd
index 1d3a1d3..1c0537c 100644
--- a/man/Rerun_Fn.Rd
+++ b/man/Rerun_Fn.Rd
@@ -23,7 +23,7 @@ Rerun_Fn(parhat0, turnoff_pars, loc_x,
 
 \item{figname}{name for figure to plot density in counter-factual scenario}
 
-\item{Map}{OPTIONAL, a tagged list of parameters to either mirror or turn off}
+\item{Map}{OPTIONAL, a tagged list of parameters to either mirror or turn off, using standard TMB interface. This input is useful, e.g., to build a model without estimating parameters, extracting Map from the list of outputs, modifying it manually, and then passing it explicitly, \code{make_model(Map=NewMap)}}
 
 \item{MapDetails_List}{output from \code{FishStatsUtils::MapDetails_Fn}}
 
diff --git a/man/calculate_proportion.Rd b/man/calculate_proportion.Rd
index 38f6fd0..04b98f5 100644
--- a/man/calculate_proportion.Rd
+++ b/man/calculate_proportion.Rd
@@ -16,7 +16,7 @@ calculate_proportion(TmbData, Index, Expansion_cz = NULL,
 
 \item{Index}{output from \code{FishStatsUtils::plot_biomass_index}}
 
-\item{Expansion_cz}{matrix specifying how densities are expanded when calculating annual indices, with a row for each category \code{c} and two columns.  The first column specifies whether to calculate annual index for category \code{c} as the weighted-sum across density estimates, where density is weighted by area ("area-weighted expansion", \code{Expansion[c,1]=0}, the default) or where density is weighted by the expanded value for another category ("abundance weighted expansion" \code{Expansion[c1,1]=1}).  The 2nd column is only used when \code{Expansion[c1,1]=1}, and specifies the category to use for abundance-weighted expansion, where \code{Expansion[c1,2]=c2} and \code{c2} must be lower than \code{c1}.}
+\item{Expansion_cz}{matrix specifying how densities are expanded when calculating annual indices, with a row for each category \code{c} and two columns.  The first column specifies whether to calculate annual index for category \code{c} as the weighted-sum across density estimates, where density is weighted by area ("area-weighted expansion", \code{Expansion[c,1]=0}, the default), where density is weighted by the expanded value for another category ("abundance weighted expansion" \code{Expansion[c1,1]=1}), or the index is calculated as the weighted average of density weighted by the expanded value for another category ("abundance weighted-average expansion" \code{Expansion[c1,1]=2}).  The 2nd column is only used when \code{Expansion[c1,1]=1} or \code{Expansion[c1,1]=2} , and specifies the category to use for abundance-weighted expansion, where \code{Expansion[c1,2]=c2} and \code{c2} must be lower than \code{c1}.}
 
 \item{Year_Set}{Year names for labeling panels}
 
diff --git a/man/combine_lists.Rd b/man/combine_lists.Rd
index b27365a..4d70611 100644
--- a/man/combine_lists.Rd
+++ b/man/combine_lists.Rd
@@ -4,7 +4,8 @@
 \alias{combine_lists}
 \title{Combine lists}
 \usage{
-combine_lists(default, input, args_to_use = NULL)
+combine_lists(default, input, args_to_use = NULL,
+  use_partial_matching = FALSE)
 }
 \arguments{
 \item{default}{default values for combined list}
diff --git a/man/convert_shapefile.Rd b/man/convert_shapefile.Rd
new file mode 100644
index 0000000..425533e
--- /dev/null
+++ b/man/convert_shapefile.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/convert_shapefile.R
+\name{convert_shapefile}
+\alias{convert_shapefile}
+\title{Convert shapefile to extrapolation-grid}
+\usage{
+convert_shapefile(file_path, projargs = NULL, grid_dim_km = c(2, 2),
+  projargs_for_shapefile = NULL, make_plots = FALSE, quiet = TRUE,
+  area_tolerance = 0.05, ...)
+}
+\arguments{
+\item{file_path}{path for shapefile on harddrive}
+
+\item{projargs}{A character string of projection arguments used to project the shapefile prior to construction an extrapolation-grid; the arguments must be entered exactly as in the PROJ.4 documentation; Default \code{projargs=NULL} uses UTM and infers the UTM zone based on the centroid of the shapefile.}
+
+\item{projargs_for_shapefile}{projection-arguments (e.g., as parsed by \code{sp::CRS}), that are used when reading shapefile and overriding any projection arguments already saved there;  Default \code{projargs_for_shapefile=NULL} uses the projection arguments available in the shapefile}
+
+\item{make_plots}{Boolean indicating whether to visualize inputs and outputs as maps}
+}
+\value{
+extrapolation-grid
+}
+\description{
+\code{convert_shapefile} reads in a shapefile and creates an extrapolation with a chosen resolution
+}
+\examples{
+\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 )
+}
+}
diff --git a/man/fit_model.Rd b/man/fit_model.Rd
index d002148..30a98b8 100644
--- a/man/fit_model.Rd
+++ b/man/fit_model.Rd
@@ -8,10 +8,11 @@ fit_model(settings, Lat_i, Lon_i, t_iz, b_i, a_i, c_iz = rep(0,
   length(b_i)), v_i = rep(0, length(b_i)),
   working_dir = paste0(getwd(), "/"), Xconfig_zcp = NULL,
   covariate_data, formula = ~0, Q_ik = NULL, newtonsteps = 1,
-  silent = TRUE, run_model = TRUE, test_fit = TRUE, ...)
+  silent = TRUE, build_model = TRUE, run_model = TRUE,
+  test_fit = TRUE, ...)
 }
 \arguments{
-\item{settings}{Output from \code{make_settings}}
+\item{settings}{Output from \code{\link{make_settings}}}
 
 \item{Lat_i}{Latitude for each sample}
 
@@ -43,20 +44,56 @@ fit_model(settings, Lat_i, Lon_i, t_iz, b_i, a_i, c_iz = rep(0,
 
 \item{newtonsteps}{number of extra newton steps to take after optimization (alternative to \code{loopnum})}
 
+\item{build_model}{Boolean indicating whether to build the model, \code{build_model=TRUE}, or simply build the inputs, \code{build_model=FALSE}}
+
 \item{run_model}{Boolean indicating whether to run the model or simply return the inputs and built TMB object}
 
 \item{test_fit}{Boolean indicating whether to apply \code{VAST::check_fit} before calculating standard errors, to test for parameters hitting bounds etc; defaults to TRUE}
 
-\item{...}{additional arguments to pass to \code{FishStatsUtils::make_extrapolation_info}, \code{FishStatsUtils::make_spatial_info}, \code{VAST::make_data}, \code{VAST::make_model}, or \code{TMBhelper::fit_tmb}, where arguments are matched by name against each function.  If an argument doesn't match, it is still passed to \code{VAST::make_data}}
+\item{...}{additional arguments to pass to \code{\link{make_extrapolation_info}}, \code{\link{make_spatial_info}}, \code{\link[VAST]{make_data}}, \code{\link[VAST]{make_model}}, or \code{\link[TMBhelper]{fit_tmb}}, where arguments are matched by name against each function.  If an argument doesn't match, it is still passed to \code{\link[VAST]{make_data}}}
 }
 \value{
-Returns a tagged list of internal objects, the TMB object, and slot \code{parameter_estimates} containing the MLE estimates
+Object of class \code{fit_model}, containing formatted inputs and outputs from VAST
+\describe{
+\item{parameter_estimates}{Output from \code{\link[TMBhelper]{fit_tmb}}; see that documentation for definition of contents}
+\item{extrapolation_list}{Output from \code{\link{make_extrapolation_info}}; see that documentation for definition of contents}
+\item{spatial_list}{Output from \code{\link{make_spatial_info}}; see that documentation for definition of contents}
+\item{data_list}{Output from \code{\link[VAST]{make_data}}; see that documentation for definition of contents}
+\item{tmb_list}{Output from \code{\link[VAST]{make_model}}; see that documentation for definition of contents}
+\item{ParHat}{Tagged list of maximum likelihood estimatesion of fixed effects and empirical Bayes estimates of random effects, following format of initial values generated by \code{\link[VAST]{make_parameters}}; see that documentation for definition of contents}
+\item{Report}{Tagged list of VAST outputs. For example, estimated density for grid \code{g}, category \code{c}, and time \code{y} is available as \code{fit$Report$D_gcy[g,c,y]}; see Details section for description of indexing}
+}
 }
 \description{
 \code{fit_model} fits a spatio-temporal model to data
 }
 \details{
-This function is the user-interface for the functions that determine the extrapolation-grid, define spatial objects, build covariates from a formula interface, assemble data, build model, estimate parameters, and check for obvious problems with the estimates.
+This function is the user-interface for the multiple mid-level functions that
+perform separate components of a spatio-temporal analysis:
+\itemize{
+\item determine the extrapolation-grid \code{\link{make_extrapolation_info}},
+\item define spatial objects \code{\link{make_spatial_info}},
+\item build covariates from a formula interface \code{\link{make_covariates}},
+\item assemble data \code{\link{make_data}},
+\item build model \code{\link{make_model}},
+\item estimate parameters \code{\link[TMBhelper]{fit_tmb}}, and
+\item check for obvious problems with the estimates \code{\link{check_fit}}.
+}
+Please see reference documetation for each of those functions (e.g., \code{?make_extrapolation_info}) to see a list of arguments used by each mid-level function.
+
+Specifically, the mid-level functions called by \code{fit_model(.)} look for arguments in the following order of precedence (from highest to lowest precedence):
+\enumerate{
+\item \code{fit_model(.)} prioritizes using named arguments passed directly to \code{fit_model(.)}. If arguments are passed this way, they are used instead of other options below.
+\item If an argument is not passed supplied directly to \code{fit_model(.)}, then \code{fit_model(.)} looks for elements in input \code{settings}, as typically created by \code{\link{make_settings}}.
+\item If an argument is not supplied via (1) or (2) above, then each mid-level function uses default values defined in those function arguments, e.g., see \code{args(make_extrapolation_info)} for defaults for function \code{make_extrapolation_info(.)}
+}
+Collectively, this order of precedence allows users to specify inputs for a specific project via input method (1), the package author to change defaults through changes in the settings
+defined for a given purpose in \code{make_settings(.)} via input method (2), while still defaulting to package defaults via option (3).
+
+Variables are indexed internally for locations \code{g}, categories \code{c}, and times \code{y}.
+Location index \code{g} represents Longitude-Latitude \code{fit$extrapolation_list$Data_Extrap[which(fit$spatial_list$g_e==g),c('Lon','Lat')]};
+Time index \code{y} represents time \code{fit$year_labels}; and
+Category \code{g} corresponds to values in \code{fit$data_list$g_i}.
 }
 \examples{
 \dontrun{
@@ -87,6 +124,8 @@ plot_results( settings=settings, fit=fit )
 \seealso{
 \code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots
 
+\code{\link{summary.fit_model}} for methods to summarize output, including obtain a dataframe of estimated densities
+
 Other wrapper functions: \code{\link{make_settings}},
   \code{\link{plot_results}}
 }
diff --git a/man/make_extrapolation_info.Rd b/man/make_extrapolation_info.Rd
index 58e50bb..c0dbc16 100644
--- a/man/make_extrapolation_info.Rd
+++ b/man/make_extrapolation_info.Rd
@@ -6,16 +6,18 @@
 \usage{
 make_extrapolation_info(Region, projargs = NA, zone = NA,
   strata.limits = data.frame(STRATA = "All_areas"),
-  create_strata_per_region = FALSE, max_cells = Inf,
+  create_strata_per_region = FALSE, max_cells = NULL,
   input_grid = NULL, observations_LL = NULL, grid_dim_km = c(2, 2),
   maximum_distance_from_sample = NULL, grid_in_UTM = TRUE,
   grid_dim_LL = c(0.1, 0.1), region = c("south_coast", "west_coast"),
   strata_to_use = c("SOG", "WCVI", "QCS", "HS", "WCHG"),
-  survey = "Chatham_rise", surveyname = "propInWCGBTS",
-  flip_around_dateline, nstart = 100, ...)
+  epu_to_use = c("All", "Georges_Bank", "Mid_Atlantic_Bight",
+  "Scotian_Shelf", "Gulf_of_Maine", "Other")[1], survey = "Chatham_rise",
+  surveyname = "propInWCGBTS", flip_around_dateline, nstart = 100,
+  area_tolerance = 0.05, ...)
 }
 \arguments{
-\item{Region}{a character vector, where each element that is matched against potential values to determine the region for the extrapolation grid. Current options are "california_current", "west_coast_hook_and_line", "british_columbia", "eastern_bering_sea", "northern_bering_sea", "bering_sea_slope", "st_matthews_island", "aleutian_islands", "gulf_of_alaska", "northwest_atlantic", "south_africa", "gulf_of_st_lawrence", "new_zealand", "habcam", "gulf_of_mexico", "stream_network", "user", or "other"}
+\item{Region}{a character vector, where each element that is matched against potential values to determine the region for the extrapolation grid. Current options are "california_current", "west_coast_hook_and_line", "british_columbia", "eastern_bering_sea", "northern_bering_sea", "bering_sea_slope", "st_matthews_island", "aleutian_islands", "gulf_of_alaska", "northwest_atlantic", "south_africa", "gulf_of_st_lawrence", "new_zealand", "habcam", "gulf_of_mexico", "ATL-IBTS-Q1", "ATL-IBTS-Q4", "BITS", "BTS", "BTS-VIIA", "EVHOE", "IE-IGFS", "NIGFS", "NS_IBTS", "PT-IBTS", "SP-ARSA", "SP-NORTH", "SP-PORC", "stream_network", "user", "other", or the absolute path and file name for a GIS shapefile}
 
 \item{projargs}{A character string of projection arguments; the arguments must be entered exactly as in the PROJ.4 documentation; if the projection is unknown, use \code{as.character(NA)}, it may be missing or an empty string of zero length and will then set to the missing value.}
 
@@ -27,9 +29,9 @@ make_extrapolation_info(Region, projargs = NA, zone = NA,
 
 \item{max_cells}{Maximum number of extrapolation-grid cells.  If number of cells in extrapolation-grid is less than this number, then its value is ignored.  Default \code{max_cells=Inf} results in no reduction in number of grid cells from the default extrapolation-grid for a given region.  Using a lower value is particularly useful when \code{fine_scale=TRUE} and using epsilon bias-correction, such that the number of extrapolation-grid cells is often a limiting factor in estimation speed.}
 
-\item{input_grid}{a matrix with three columns (labeled 'Lat', 'Lon', and 'Area_km2') giving latitude, longitude, and area for each cell of a user-supplied grid; only used when \code{Region="other"}}
+\item{input_grid}{a matrix with three columns (labeled 'Lat', 'Lon', and 'Area_km2') giving latitude, longitude, and area for each cell of a user-supplied grid; only used when \code{Region="user"}}
 
-\item{observations_LL}{a matrix with two columns (labeled 'Lat' and 'Lon') giving latitude and longitude for each observation; only used when \code{Region="user"}}
+\item{observations_LL}{a matrix with two columns (labeled 'Lat' and 'Lon') giving latitude and longitude for each observation; only used when \code{Region="other"}}
 
 \item{grid_dim_km}{numeric-vector with length two, giving the distance in km between cells in the automatically generated extrapolation grid; only used if \code{Region="other"}}
 
@@ -43,6 +45,8 @@ make_extrapolation_info(Region, projargs = NA, zone = NA,
 
 \item{strata_to_use}{strata to include by default for the BC coast extrapolation grid; only used if \code{Region="british_columbia"}}
 
+\item{epu_to_use}{EPU to include for the Northwest Atlantic (NWA) extrapolation grid, default is "All"; only used if \code{Region="northwest_atlantic"}}
+
 \item{survey}{survey to use for New Zealand extrapolation grid; only used if \code{Region="new_zealand"}}
 
 \item{surveyname}{area of West Coast to include in area-weighted extrapolation for California Current; only used if \code{Region="california_current"}}
@@ -67,5 +71,18 @@ Tagged list used in other functions
 \code{make_extrapolation_data} builds an object used to determine areas to extrapolation densities to when calculating indices
 }
 \details{
-To do area-weighted extrapolation of estimated density for use in calculating abundance indices, it is necessary to have a precise measurement of the footprint for a given survey design. Using VAST, analysts do this by including an "extrapolation grid" where densities are predicted at the location of each grid cell and where each grid cell is associated with a known area within a given survey design. Collaborators have worked with the package author to include the extrapolation-grid for several regions automatically in FishStatsUtils, but for new regions an analyst must either detect the grid automatically using \code{Region="Other"} or input an extrapolation-grid manually using \code{Region="User"}.  The extrapolation is also used to determine where to drawn pixels when plotting predictions of density.
+To do area-weighted extrapolation of estimated density for use in calculating abundance indices,
+it is necessary to have a precise measurement of the footprint for a given survey design.
+Using VAST, analysts do this by including an "extrapolation grid" where densities are predicted
+at the location of each grid cell and where each grid cell is associated with a known area within a given survey design.
+Collaborators have worked with the package author to include the extrapolation-grid for several
+regions automatically in FishStatsUtils. For new regions an analyst can either (1) detect
+the grid automatically using \code{Region="Other"}, or (2) input an extrapolation-grid manually
+using \code{Region="User"}, or supply a GIS shapefile \code{Region="[directory_path/file_name].shp"}.
+The extrapolation is also used to determine where to drawn pixels when plotting predictions of density.
+If a user supplies a character-vector with more than one of these, then they are combined to
+assemble a combined extrapolation-grid.
+
+When supplying a shapefile, I recommend using a UTM projection for projargs, which appears to have lower
+projection errors regarding total area than rnaturalearth.
 }
diff --git a/man/make_settings.Rd b/man/make_settings.Rd
index e7c3770..9908ed3 100644
--- a/man/make_settings.Rd
+++ b/man/make_settings.Rd
@@ -8,12 +8,13 @@ make_settings(n_x, Region, purpose = "index", fine_scale = TRUE,
   strata.limits = data.frame(STRATA = "All_areas"), zone = NA,
   FieldConfig, RhoConfig, OverdispersionConfig, ObsModel, bias.correct,
   Options, use_anisotropy, vars_to_correct, Version,
-  treat_nonencounter_as_zero, n_categories, VamConfig)
+  treat_nonencounter_as_zero, n_categories, VamConfig, max_cells,
+  knot_method)
 }
 \arguments{
 \item{n_x}{the number of vertices in the SPDE mesh (determines the spatial resolution when Method="Mesh")}
 
-\item{Region}{a character vector, where each element that is matched against potential values to determine the region for the extrapolation grid. Current options are "california_current", "west_coast_hook_and_line", "british_columbia", "eastern_bering_sea", "northern_bering_sea", "bering_sea_slope", "st_matthews_island", "aleutian_islands", "gulf_of_alaska", "northwest_atlantic", "south_africa", "gulf_of_st_lawrence", "new_zealand", "habcam", "gulf_of_mexico", "stream_network", "user", or "other"}
+\item{Region}{a character vector, where each element that is matched against potential values to determine the region for the extrapolation grid. Current options are "california_current", "west_coast_hook_and_line", "british_columbia", "eastern_bering_sea", "northern_bering_sea", "bering_sea_slope", "st_matthews_island", "aleutian_islands", "gulf_of_alaska", "northwest_atlantic", "south_africa", "gulf_of_st_lawrence", "new_zealand", "habcam", "gulf_of_mexico", "ATL-IBTS-Q1", "ATL-IBTS-Q4", "BITS", "BTS", "BTS-VIIA", "EVHOE", "IE-IGFS", "NIGFS", "NS_IBTS", "PT-IBTS", "SP-ARSA", "SP-NORTH", "SP-PORC", "stream_network", "user", "other", or the absolute path and file name for a GIS shapefile}
 
 \item{purpose}{character indicating what purpose is intended for the model, and therefore what default settings are perhaps appropriate. Only currently implemented for \code{purpose="index"}, \code{purpose="condition_and_density"}, \code{purpose="MICE"}, or \code{purpose="ordination"}.}
 
@@ -47,6 +48,10 @@ make_settings(n_x, Region, purpose = "index", fine_scale = TRUE,
   \item{VamConfig[1]}{indicates the rank of the interaction matrix, indicating the number of community axes that have regulated dynamics}
   \item{VamConfig[2]}{Indicates whether interactions occur before spatio-temporal variation (\code{VamConfig[2]=0}) or after \code{VamConfig[2]=1}}
 }}
+
+\item{max_cells}{Maximum number of extrapolation-grid cells.  If number of cells in extrapolation-grid is less than this number, then its value is ignored.  Default \code{max_cells=Inf} results in no reduction in number of grid cells from the default extrapolation-grid for a given region.  Using a lower value is particularly useful when \code{fine_scale=TRUE} and using epsilon bias-correction, such that the number of extrapolation-grid cells is often a limiting factor in estimation speed.}
+
+\item{knot_method}{whether to determine location of GMRF vertices based on the location of samples \code{knot_method=`samples`} or extrapolation-grid cells within the specified strata \code{knot_method='grid'}; default \code{knot_method=NULL} is coerced to \code{knot_method=`samples`}}
 }
 \value{
 Tagged list containing default settings for a given purpose, use \code{names} on output to see list of settings.
diff --git a/man/make_spatial_info.Rd b/man/make_spatial_info.Rd
index 0a4eeef..46adb81 100644
--- a/man/make_spatial_info.Rd
+++ b/man/make_spatial_info.Rd
@@ -5,7 +5,7 @@
 \title{Build objects related to spatial information}
 \usage{
 make_spatial_info(n_x, Lon_i, Lat_i, Extrapolation_List,
-  knot_method = "samples", Method = "Mesh", grid_size_km = 50,
+  knot_method = NULL, Method = "Mesh", grid_size_km = 50,
   grid_size_LL = 1, fine_scale = FALSE, Network_sz_LL = NULL,
   iter.max = 1000, randomseed = 1, nstart = 100,
   DirPath = paste0(getwd(), "/"), Save_Results = FALSE, LON_intensity,
@@ -20,7 +20,7 @@ make_spatial_info(n_x, Lon_i, Lat_i, Extrapolation_List,
 
 \item{Extrapolation_List}{the output from \code{Prepare_Extrapolation_Data_Fn}}
 
-\item{knot_method}{whether to determine location of GMRF vertices based on the location of samples \code{knot_method=`samples`} or extrapolation-grid cells within the specified strata \code{knot_method='grid'}}
+\item{knot_method}{whether to determine location of GMRF vertices based on the location of samples \code{knot_method=`samples`} or extrapolation-grid cells within the specified strata \code{knot_method='grid'}; default \code{knot_method=NULL} is coerced to \code{knot_method=`samples`}}
 
 \item{Method}{a character of either "Grid" or "Mesh" where "Grid" is a 2D AR1 process, and "Mesh" is the SPDE method with geometric anisotropy}
 
diff --git a/man/plot_data.Rd b/man/plot_data.Rd
index 4e1bfdf..98c6fa8 100644
--- a/man/plot_data.Rd
+++ b/man/plot_data.Rd
@@ -10,7 +10,7 @@ plot_data(Extrapolation_List, Spatial_List, Data_Geostat = NULL,
   Plot1_name = "Data_and_knots.png", Plot2_name = "Data_by_year.png",
   col = "red", cex = 0.01, pch = 19, Year_Set,
   projargs = "+proj=longlat", map_resolution = "medium",
-  land_color = "grey", ...)
+  land_color = "grey", country = NULL, ...)
 }
 \arguments{
 \item{Extrapolation_List}{Output from \code{Prepare_Extrapolation_Data_Fn}}
@@ -25,6 +25,8 @@ plot_data(Extrapolation_List, Spatial_List, Data_Geostat = NULL,
 
 \item{land_color}{color for filling in land (use \code{land_color=rgb(0,0,0,alpha=0)} for transparent land)}
 
+\item{country}{a character vector of country names.}
+
 \item{...}{addition inputs to \code{plot}}
 }
 \description{
diff --git a/man/plot_factors.Rd b/man/plot_factors.Rd
index 511ebee..53dc573 100644
--- a/man/plot_factors.Rd
+++ b/man/plot_factors.Rd
@@ -7,7 +7,8 @@
 plot_factors(Report, ParHat, Data, SD = NULL, Year_Set = NULL,
   category_names = NULL, RotationMethod = "PCA",
   mapdetails_list = NULL, Dim_year = NULL, Dim_species = NULL,
-  plotdir = paste0(getwd(), "/"), land_color = "grey", ...)
+  plotdir = paste0(getwd(), "/"), land_color = "grey", zlim = NA,
+  testcutoff = 1e-04, ...)
 }
 \arguments{
 \item{Report}{output report, e.g., from \code{Report <- obj$report()}}
@@ -22,6 +23,8 @@ plot_factors(Report, ParHat, Data, SD = NULL, Year_Set = NULL,
 
 \item{category_names}{character-vector listing name for each category}
 
+\item{RotationMethod}{Method used for rotation, Options: "PCA" (recommended) or "Varimax"}
+
 \item{mapdetails_list}{output from \code{FishStatsUtils::MapDetails_Fn}}
 
 \item{Dim_year}{Plotting dimension (row,column) for plot of years (default: square with sufficient size for number of years)}
@@ -32,7 +35,11 @@ plot_factors(Report, ParHat, Data, SD = NULL, Year_Set = NULL,
 
 \item{land_color}{color for filling in land (use \code{land_color=rgb(0,0,0,alpha=0)} for transparent land)}
 
-\item{...}{additional arguments passed to \code{plot_maps(.)} when plotting spatio-temporal terms Epsilon}
+\item{zlim}{range for defining bounds of color scale.  If \code{zlim=NULL}, then a constant scale is inferred from the range of \code{Y_gt} and a color-legend is plotted in the last panel.  If \code{zlim=NA} then a different range is used in each panel from the range of \code{Y_gt[,t]} and a color-legend is plotted in every panel.}
+
+\item{testcutoff}{tolerance for numerical rounding when confirming that rotation doesn't effect results}
+
+\item{...}{additional arguments passed to \code{plot_maps(.)} and/or \code{plot_variable(.)} when plotting factor-values on a map}
 }
 \description{
 \code{plot_factors} plots factor loadings, average spatial factors, and spatio-temporal factors
diff --git a/man/plot_loadings.Rd b/man/plot_loadings.Rd
index e168150..859051c 100644
--- a/man/plot_loadings.Rd
+++ b/man/plot_loadings.Rd
@@ -4,8 +4,8 @@
 \alias{plot_loadings}
 \title{Plotting loadings matrix}
 \usage{
-plot_loadings(L_pj, whichfactor = 1, addtitle = TRUE,
-  LabelPosition = "Right", Buffer = c(0, 0.1),
+plot_loadings(L_pj, At = 1:nrow(L_pj), whichfactor = 1,
+  addtitle = TRUE, LabelPosition = "Right", Buffer = c(0, 0.1),
   Labels = rownames(L_pj), Cex = 1.2,
   legend_text = "Proportion of explained variance", ...)
 }
diff --git a/man/plot_maps.Rd b/man/plot_maps.Rd
index 797d997..3548138 100644
--- a/man/plot_maps.Rd
+++ b/man/plot_maps.Rd
@@ -66,6 +66,10 @@ plot_maps(
 
 \item{Report}{tagged list of outputs from TMB model via \code{Obj$report()}}
 
+\item{zlim}{range for defining bounds of color scale.  If \code{zlim=NULL}, then a constant scale is inferred from the range of \code{Y_gt} and a color-legend is plotted in the last panel.  If \code{zlim=NA} then a different range is used in each panel from the range of \code{Y_gt[,t]} and a color-legend is plotted in every panel.}
+
+\item{country}{a character vector of country names.}
+
 \item{...}{arguments passed to \code{FishStatsUtils::plot_variable}}
 
 \item{country}{optional list of countries to display, e.g. c("united states of america", "canada"). If maps are generating visual artefacts, please try using argument \code{country} to simplify the polygons used to represent land features.}
diff --git a/man/plot_range_edge.Rd b/man/plot_range_edge.Rd
index 7a297fe..caa9b9a 100644
--- a/man/plot_range_edge.Rd
+++ b/man/plot_range_edge.Rd
@@ -8,7 +8,7 @@ plot_range_edge(Sdreport, Obj, Year_Set = NULL, Years2Include = NULL,
   strata_names = NULL, category_names = NULL,
   working_dir = paste0(getwd(), "/"), quantiles = c(0.05, 0.95),
   n_samples = 100, interval_width = 1, width = NULL, height = NULL,
-  ...)
+  calculate_relative_to_average = FALSE, seed = 123456, ...)
 }
 \arguments{
 \item{Sdreport}{Standard deviation outputs from TMB model via \code{sdreport(Obj)}}
@@ -35,6 +35,10 @@ plot_range_edge(Sdreport, Obj, Year_Set = NULL, Years2Include = NULL,
 
 \item{height}{plot height in inches}
 
+\item{calculate_relative_to_average}{Boolean, whether to calculate edge in UTM coordinates (default), or instead calculate relative to median across all years. The latter reduces standard errors, and is appropriate when checking significance for comparison across years for a single species.  The former (default) is appropriate for checking significance for comparison across species.}
+
+\item{seed}{integer used to set random-number seed when sampling variables, as passed to \code{set.seed(.)}}
+
 \item{...}{Other inputs to `par()`}
 }
 \description{
diff --git a/man/plot_results.Rd b/man/plot_results.Rd
index 92ddfbe..9927e72 100644
--- a/man/plot_results.Rd
+++ b/man/plot_results.Rd
@@ -8,12 +8,12 @@ plot_results(fit, settings = fit$settings, plot_set = 3,
   working_dir = paste0(getwd(), "/"), year_labels = fit$year_labels,
   years_to_plot = fit$years_to_plot, use_biascorr = TRUE, map_list,
   category_names, check_residuals = TRUE, projargs = "+proj=longlat",
-  zrange, n_samples = 100, ...)
+  zrange, n_samples = 100, calculate_relative_to_average = FALSE, ...)
 }
 \arguments{
 \item{fit}{Output from \code{fit_model}}
 
-\item{settings}{Output from \code{make_settings}}
+\item{settings}{Output from \code{\link{make_settings}}}
 
 \item{plot_set}{integer-vector defining plots to create
 \describe{
@@ -45,16 +45,28 @@ plot_results(fit, settings = fit$settings, plot_set = 3,
 
 \item{n_samples}{number of samples from the joint predictive distribution for fixed and random effects.  Default is 100, which is slow.}
 
+\item{calculate_relative_to_average}{Boolean, whether to calculate edge in UTM coordinates (default), or instead calculate relative to median across all years. The latter reduces standard errors, and is appropriate when checking significance for comparison across years for a single species.  The former (default) is appropriate for checking significance for comparison across species.}
+
 \item{...}{additional settings to pass to \code{FishStatsUtils::plot_maps}}
 }
 \value{
 Invisibly returns a tagged list of outputs generated by standard plots.
+\describe{
+  \item{Q}{Output from \code{\link{plot_quantile_diagnostic}}; see that documentation for definition of optional arguments and contents}
+  \item{Enc_prob}{Output from \code{\link{plot_encounter_diagnostic}}; see that documentation for definition of optional arguments and contents}
+  \item{Index}{Output from \code{\link{plot_biomass_index}}; see that documentation for definition of optional arguments and contents}
+  \item{Proportions}{Output from \code{\link{calculate_proportion}}; see that documentation for definition of optional arguments and contents}
+  \item{Range}{Output from \code{\link{plot_range_index}}; see that documentation for definition of optional arguments and contents}
+  \item{Dens_xt}{Output from \code{\link{plot_maps}}; see that documentation for definition of optional arguments and contents}
+  \item{Edge}{Output from \code{\link{plot_range_edge}}; see that documentation for definition of optional arguments and contents}
+}
 }
 \description{
 \code{plot_results} plots diagnostics, results, and indices for a given fitted model
 }
 \details{
 This function takes a fitted VAST model and generates a standard set of diagnostic and visualization plots.
+It does this by calling a series of mid-level plotting functions; see list of functions in Value section of documentation.
 }
 \seealso{
 \code{?VAST} for general documentation, \code{?make_settings} for generic settings, \code{?fit_model} for model fitting, and \code{?plot_results} for generic plots
diff --git a/man/plot_variable.Rd b/man/plot_variable.Rd
index 823e2b9..4ca7605 100644
--- a/man/plot_variable.Rd
+++ b/man/plot_variable.Rd
@@ -7,10 +7,10 @@
 plot_variable(Y_gt, map_list, panel_labels, projargs = "+proj=longlat",
   map_resolution = "medium", file_name = "density",
   working_dir = paste0(getwd(), "/"), Format = "png", Res = 200,
-  add = FALSE, outermargintext = c("Eastings", "Northings"), zlim, col,
-  mar = c(0, 0, 2, 0), oma = c(4, 4, 0, 0), legend_x = c(0, 0.05),
-  legend_y = c(0.05, 0.45), cex.legend = 1, mfrow,
-  land_color = "grey", n_cells, xlim, ylim, country = NULL, ...)
+  add = FALSE, outermargintext = c("Eastings", "Northings"),
+  zlim = NULL, col, mar = c(0, 0, 2, 0), oma = c(4, 4, 0, 0),
+  legend_x = c(0, 0.05), legend_y = c(0.05, 0.45), cex.legend = 1,
+  mfrow, land_color = "grey", n_cells, xlim, ylim, country = NULL, ...)
 }
 \arguments{
 \item{Y_gt}{matrix where values for every column are plotted as a map}
@@ -25,7 +25,7 @@ plot_variable(Y_gt, map_list, panel_labels, projargs = "+proj=longlat",
 
 \item{outermargintext}{vector defining text to plot in outer margins of panel figure}
 
-\item{zlim}{range for defining bounds of color scale}
+\item{zlim}{range for defining bounds of color scale.  If \code{zlim=NULL}, then a constant scale is inferred from the range of \code{Y_gt} and a color-legend is plotted in the last panel.  If \code{zlim=NA} then a different range is used in each panel from the range of \code{Y_gt[,t]} and a color-legend is plotted in every panel.}
 
 \item{legend_x}{two numeric values (generally between 0 and 1, but slightly lower/higher values generate colorbars that are just outside the plotting window) giving left and right-hand location of color legend}
 
@@ -40,5 +40,5 @@ plot_variable(Y_gt, map_list, panel_labels, projargs = "+proj=longlat",
 \description{
 \code{plot_variable} plots a map and fills in regions with colors to represent intensity in an areal-interpretion of model results
 
-See \url{https://proj.org/operations/projections/index.html} for a list of projections to pass via \code{projargs}. I often prefer \code{projargs='+proj=natearth +lat_0=0 +units=km'} where argument \code{+lat_0} allows the user to center eastings on a specified latitude. If maps are generating visual artefacts, please try using argument \code{country} to simplify the polygons used to represent land features.
+See \url{https://proj.org/operations/projections/index.html} for a list of projections to pass via \code{projargs}. I often prefer \code{projargs='+proj=natearth +lon_0=0 +units=km'} where argument \code{+lon_0} allows the user to center eastings on a specified longitude. If maps are generating visual artefacts, please try using argument \code{country} to simplify the polygons used to represent land features.
 }
diff --git a/man/rotate_factors.Rd b/man/rotate_factors.Rd
index 9f43a73..34b697c 100644
--- a/man/rotate_factors.Rd
+++ b/man/rotate_factors.Rd
@@ -4,7 +4,7 @@
 \alias{rotate_factors}
 \title{Rotate results}
 \usage{
-rotate_factors(Cov_jj = NULL, L_pj = NULL, Psi = NULL,
+rotate_factors(Cov_jj = NULL, L_pj = NULL, Psi_sjt = NULL,
   RotationMethod = "PCA", testcutoff = 1e-10, quiet = FALSE)
 }
 \arguments{
@@ -12,7 +12,7 @@ rotate_factors(Cov_jj = NULL, L_pj = NULL, Psi = NULL,
 
 \item{L_pj}{Loadings matrix for `p` categories and `j` factors (calculated from \code{Cov_jj} if it is provided)}
 
-\item{Psi}{Array of factors (1st dimension: spatial knots;  2nd dimension: factors;  3rd dimension:  time)}
+\item{Psi_sjt}{Array of factors (1st dimension: spatial knots;  2nd dimension: factors;  3rd dimension:  time)}
 
 \item{RotationMethod}{Method used for rotation, Options: "PCA" (recommended) or "Varimax"}
 
diff --git a/man/sample_variable.Rd b/man/sample_variable.Rd
index 52cc160..0aa04ca 100644
--- a/man/sample_variable.Rd
+++ b/man/sample_variable.Rd
@@ -4,7 +4,8 @@
 \alias{sample_variable}
 \title{Sample from predictive distribution of a variable}
 \usage{
-sample_variable(Sdreport, Obj, variable_name, n_samples = 100)
+sample_variable(Sdreport, Obj, variable_name, n_samples = 100,
+  seed = 123456)
 }
 \arguments{
 \item{Sdreport}{TMB output from `\code{TMB::sdreport(Obj)}`}
@@ -14,6 +15,8 @@ sample_variable(Sdreport, Obj, variable_name, n_samples = 100)
 \item{variable_name}{name of variable available in report using \code{Obj$report()}}
 
 \item{n_samples}{number of samples from the joint predictive distribution for fixed and random effects.  Default is 100, which is slow.}
+
+\item{seed}{integer used to set random-number seed when sampling variables, as passed to \code{set.seed(.)}}
 }
 \description{
 \code{sample_variable} samples from the joint distribution of random effects to approximate the predictive distribution for a variable
diff --git a/man/simulate_data.Rd b/man/simulate_data.Rd
new file mode 100644
index 0000000..0aee566
--- /dev/null
+++ b/man/simulate_data.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/simulate_data.R
+\name{simulate_data}
+\alias{simulate_data}
+\title{Simulate new dynamics and sampling data}
+\usage{
+simulate_data(fit, type = 1, random_seed = NULL)
+}
+\arguments{
+\item{fit}{output form \code{fit_model(.)}}
+
+\item{type}{integer stating what type of simulation to use.  \code{type=1} simulates new data conditional upon estimated fixed and random effects.  \code{type=2} simulates new random effects conditional upon fixed effects, and new data conditional upon both.  \code{type=3} simulates new fixed and random effects from the joint precision matrix, and new data conditional upon these values.}
+
+\item{random_seed}{integer passed to \code{\link[base]{set.seed}} whenever \code{type=3}, where the default value \code{random_seed=NULL} resets the random-number seed. Argument no effect when \code{type!=3} because TMB has no interface for setting the random-number seed in C++}
+}
+\value{
+Report object containing new data and population variables including
+\describe{
+  \item{b_i}{New simulated data}
+  \item{D_gcy}{Density for each grid cell g, category c, and year y}
+  \item{Index_cyl}{Index of abundance for each category c, year y, and stratum l}
+}
+}
+\description{
+\code{simulate_data} conducts a parametric bootstrap to simulate new data and potentially simulate new population dynamics and associated variables
+}
diff --git a/man/summary.fit_model.Rd b/man/summary.fit_model.Rd
index c4eb060..22db972 100644
--- a/man/summary.fit_model.Rd
+++ b/man/summary.fit_model.Rd
@@ -7,11 +7,15 @@
 \method{summary}{fit_model}(x, what = "density", ...)
 }
 \arguments{
+\item{x}{Output from \code{\link{fit_model}}}
+
 \item{what}{Boolean indicating what to summarize; only option is `density`}
 
-\item{...}{Not used}
+\item{...}{Not used
 
-\item{fit}{Output from \code{\link{fit_model}}}
+\code{what="density"} returns a tagged list containing element \code{Density_dataframe},
+which lists the estimated density for every Latitude-Longitude-Year-Category combination
+for every modelled location in the extrapolation-grid.}
 }
 \description{
 Extract summary of spatial estimates