diff --git a/.github/.gitignore b/.github/.gitignore
index 2d19fc7..0c817ed 100644
--- a/.github/.gitignore
+++ b/.github/.gitignore
@@ -1 +1,2 @@
 *.html
+*.rds
diff --git a/.github/workflows/call-doc-and-style-r.yml b/.github/workflows/call-doc-and-style-r.yml
new file mode 100644
index 0000000..2e93be3
--- /dev/null
+++ b/.github/workflows/call-doc-and-style-r.yml
@@ -0,0 +1,10 @@
+# document and style R code using a reusable workflow
+name: call-doc-and-style-r
+# on specifies the build triggers. See more info at https://docs.github.com/en/actions/learn-github-actions/events-that-trigger-workflows
+on:
+  push:
+    branches: [main]
+jobs:
+  call-workflow:
+    uses: nmfs-fish-tools/ghactions4r/.github/workflows/doc-and-style-r.yml@main
+  
diff --git a/R/calc_index_areas.R b/R/calc_index_areas.R
index bab9877..7b662df 100644
--- a/R/calc_index_areas.R
+++ b/R/calc_index_areas.R
@@ -44,13 +44,14 @@
 #' list("test" = c(39, 38))
 #' # run the function
 #' index <- calc_index_areas(
-#'   data, fit, grid, dir = getwd(), boundaries = list("test" = c(39, 38))
+#'   data, fit, grid,
+#'   dir = getwd(), boundaries = list("test" = c(39, 38))
 #' )
 #' # look at the index, which will be in "test" because that is what we
 #' # named the boundary
 #' index[["test"]][["index"]]
 #' }
-#' 
+#'
 calc_index_areas <- function(data,
                              fit,
                              prediction_grid,
diff --git a/R/data.R b/R/data.R
index db49a36..433ed3c 100644
--- a/R/data.R
+++ b/R/data.R
@@ -33,4 +33,4 @@
 #' * South of Cape Mendocino
 #' * North of Monterey Bay
 #' * South of Monterey Bay
-"boundaries_data"
\ No newline at end of file
+"boundaries_data"
diff --git a/R/diagnose.R b/R/diagnose.R
index 30baf34..8dc9171 100644
--- a/R/diagnose.R
+++ b/R/diagnose.R
@@ -81,7 +81,7 @@ diagnose <- function(dir,
     fit[["data"]][["residuals2"]] <- stats::residuals(
       fit,
       model = 2,
-      type = "mle-mvn"  
+      type = "mle-mvn"
     )
   }
 
diff --git a/R/filter_grid.R b/R/filter_grid.R
index 13aa93c..b774ede 100644
--- a/R/filter_grid.R
+++ b/R/filter_grid.R
@@ -6,7 +6,7 @@
 #' @return A data frame with the same structure as `grid`.
 #' @author Kelli F. Johnson
 filter_grid <- function(grid, boundary_north, boundary_south) {
-  sub_grid <- grid %>% 
+  sub_grid <- grid %>%
     dplyr::filter(latitude < boundary_north) %>%
-    dplyr::filter(latitude >= boundary_south)  
+    dplyr::filter(latitude >= boundary_south)
 }
diff --git a/R/filter_odfw.R b/R/filter_odfw.R
index 0c0de7a..58741ff 100644
--- a/R/filter_odfw.R
+++ b/R/filter_odfw.R
@@ -1,13 +1,11 @@
-
-
 # ##### AT SEA OBSERVER CPUE DATA PREP/FILTERING
 
-# # created by A. Whitman 
-# # updated 06/21/2021 
+# # created by A. Whitman
+# # updated 06/21/2021
 
 # # last update to exclude 2020 due to sample size.
 
-# # purpose: explore drift level data for the At Sea charter observer CPUE index 
+# # purpose: explore drift level data for the At Sea charter observer CPUE index
 
 # # load some libraries  - not sure if I'll need all these.
 
@@ -51,29 +49,29 @@
 # dat<-dat[!is.na(dat$Total_lingcod),]
 # summary(dat)
 
-# # adding a p/a variable for lingcod 
+# # adding a p/a variable for lingcod
 # dat$ling.pa<-ifelse(dat$Total_lingcod==0,0,1)
 
-# # adding a megaregion variable 
+# # adding a megaregion variable
 # dat$Megaregion<-as.factor(ifelse(dat$Port<30,"North","South"))
 
-# # positive rate of overall dataset 
-# dim(dat)[1] # total records 
-# length(which(dat$ling.pa==1)) # total positive drifts 
-# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate 
+# # positive rate of overall dataset
+# dim(dat)[1] # total records
+# length(which(dat$ling.pa==1)) # total positive drifts
+# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate
 
 # # positive rate by year without any filtering
-# pos_rate_year<-dat %>% 
+# pos_rate_year<-dat %>%
 #   dplyr::group_by(Year)%>%
-#   dplyr::summarise(No.Drifts=n(),Pos.Drifts=round((sum(ling.pa)/n()),3)) %>% 
+#   dplyr::summarise(No.Drifts=n(),Pos.Drifts=round((sum(ling.pa)/n()),3)) %>%
 #   as.data.frame()
-# pos_rate_year  
+# pos_rate_year
 
-# # ouch on 2020, might need to throw that year out, otherwise, annual rate 
+# # ouch on 2020, might need to throw that year out, otherwise, annual rate
 # # ranges from approx. 0.203 to 0.611
 
-# #looking at potential filters 
-# # drift minutes, drifts without depth, depth bins, closed season or depths 
+# #looking at potential filters
+# # drift minutes, drifts without depth, depth bins, closed season or depths
 # # look at drifts by megaregion and by depth bin, filter by percent midwater
 
 # ### FILTER: DRIFT MINUTES
@@ -91,7 +89,7 @@
 # dat<-dat[dat$DriftMinutes<=28,]
 
 # #'
-# #' 
+# #'
 # filter_odfw_drift_quantile <- function(x, percent = 0.95, verbose = FALSE) {
 #   limit <- stats::quantile(x, percent)
 #   remove <- which(x > limit)
@@ -108,48 +106,48 @@
 # }
 # test %>%
 #   dplyr::mutate(filter_drift = filter_odfw_drift_quantile(DriftMinutes))
-# # positive rate of overall dataset 
-# dim(dat)[1] # total records 
-# length(which(dat$ling.pa==1)) # total positive drifts 
-# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate 
+# # positive rate of overall dataset
+# dim(dat)[1] # total records
+# length(which(dat$ling.pa==1)) # total positive drifts
+# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate
 
 # ### FILTER: DEPTHS
 
 # summary(dat$SDEPTH_fm)
-# # remove the NAs 
+# # remove the NAs
 # dat<-dat[!is.na(dat$SDEPTH_fm),]
 # histogram(dat$SDEPTH_fm)
-# quantile(dat$SDEPTH_fm,0.99) #99 percentile is 60, so I will cut it off there 
+# quantile(dat$SDEPTH_fm,0.99) #99 percentile is 60, so I will cut it off there
 # dat<-dat[dat$SDEPTH_fm<=60,]
 
-# # positive rate of overall dataset 
-# dim(dat)[1] # total records 
-# length(which(dat$ling.pa==1)) # total positive drifts 
-# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate 
+# # positive rate of overall dataset
+# dim(dat)[1] # total records
+# length(which(dat$ling.pa==1)) # total positive drifts
+# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate
 
 # #creating a depth bin from the start depth in FMs
-# dat<- dat %>% 
-#   mutate(DepthBin=cut(dat$SDEPTH_fm,breaks = c(seq(0,60,by = 5)))) %>% 
+# dat<- dat %>%
+#   mutate(DepthBin=cut(dat$SDEPTH_fm,breaks = c(seq(0,60,by = 5)))) %>%
 #   as.data.frame()
 
-# # look at drifts by depth bin 
+# # look at drifts by depth bin
 # drifts_by_depth<-dat %>%
-#   group_by(DepthBin) %>% 
+#   group_by(DepthBin) %>%
 #   summarise(drifts=n())
 # drifts_by_depth
 # histogram(dat$DepthBin)
 
-# ### FILTER: PERCENT MIDWATER FISH 
+# ### FILTER: PERCENT MIDWATER FISH
 
 # summary(dat$percent_midwater)
 # quantile(dat$percent_midwater)  # so 50% is 0.75 midwater, and 75% is 1.0! this would be a pretty heavy filter
 # histogram(dat$percent_midwater) # lots of zeros and lots of 1.0, no clear break in between
-# length(which(dat$percent_midwater<=0.75))  
+# length(which(dat$percent_midwater<=0.75))
 
-# #just to check, are there lingcod caught on any 100% midwater trips? 
+# #just to check, are there lingcod caught on any 100% midwater trips?
 # length(which(dat$percent_midwater==1.0 & dat$ling.pa==1)) # yay
 
-# # if I filter out all the 1.0's, what do the quantiles look like then? 
+# # if I filter out all the 1.0's, what do the quantiles look like then?
 # dat_test<-dat[dat$percent_midwater<1.0,]
 # dim(dat_test)[1]
 # quantile(dat_test$percent_midwater,0.99) # that's 0.92 percent midwater, so I think that's reasonable
@@ -157,61 +155,61 @@
 
 # dat<-dat[dat$percent_midwater<=0.92,]
 
-# # positive rate of overall dataset 
-# dim(dat)[1] # total records 
-# length(which(dat$ling.pa==1)) # total positive drifts 
-# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate 
+# # positive rate of overall dataset
+# dim(dat)[1] # total records
+# length(which(dat$ling.pa==1)) # total positive drifts
+# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate
 
-# #### CHECKING OPEN DEPTHS 
+# #### CHECKING OPEN DEPTHS
 # summary(dat$GF_OpenDepth)
 
-# # so lingcod is not a LL species, so those get recategorized 
+# # so lingcod is not a LL species, so those get recategorized
 # levels(dat$GF_OpenDepth)[levels(dat$GF_OpenDepth)=="<30 fm & >40 fm (LL)"] <- "<30 fm"
 # levels(dat$GF_OpenDepth)[levels(dat$GF_OpenDepth)=="<40 fm & >40 fm (LL)"] <- "<40 fm"
 # histogram(dat$GF_OpenDepth)
 # table(dat$GF_OpenDepth,dat$ling.pa > 0) # a few trips outside of 40fm
 
-# #checking if fishing in legal depths 
+# #checking if fishing in legal depths
 # #add a legal depth from the GF_OpenDepth
 # dat$legal.dep<-as.numeric(ifelse(dat$GF_OpenDepth=="All Depth",1000,substring(dat$GF_OpenDepth,2,3)))
 # dat$dep.diff<-dat$legal.dep-dat$SDEPTH_fm # so if it's negative, they were fishing in depths that were deeper than legal
 # length(which(dat$dep.diff<0))
 
-# #### FILTER: FISHING OUTSIDE LEGAL DEPTHS 
+# #### FILTER: FISHING OUTSIDE LEGAL DEPTHS
 
-# # remove drifts outside legal depths, with a buffer 
+# # remove drifts outside legal depths, with a buffer
 # dat<-dat[dat$dep.diff>-5,]
-  
-# # positive rate of overall dataset 
-# dim(dat)[1] # total records 
-# length(which(dat$ling.pa==1)) # total positive drifts 
-# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate 
+
+# # positive rate of overall dataset
+# dim(dat)[1] # total records
+# length(which(dat$ling.pa==1)) # total positive drifts
+# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate
 
 # # positive rate by year without any filtering
-# pos_rate_year<-dat %>% 
+# pos_rate_year<-dat %>%
 #   dplyr::group_by(Year)%>%
-#   dplyr::summarise(No.Drifts=n(),Pos.Drifts=sum(ling.pa),Pos.Perc=round((sum(ling.pa)/n()),3)) %>% 
+#   dplyr::summarise(No.Drifts=n(),Pos.Drifts=sum(ling.pa),Pos.Perc=round((sum(ling.pa)/n()),3)) %>%
 #   as.data.frame()
-# pos_rate_year  
+# pos_rate_year
 
-# ### FILTER - REMOVE 2020 BECAUSE OF SAMPLE SIZE 
+# ### FILTER - REMOVE 2020 BECAUSE OF SAMPLE SIZE
 
 # dat<-dat[dat$Year %in% c(2001,2003:2019),]
 
-# # positive rate of overall dataset 
-# dim(dat)[1] # total records 
-# length(which(dat$ling.pa==1)) # total positive drifts 
-# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate 
+# # positive rate of overall dataset
+# dim(dat)[1] # total records
+# length(which(dat$ling.pa==1)) # total positive drifts
+# length(which(dat$ling.pa==1))/dim(dat)[1] # positive rate
 
-# pos_rate_year<-dat %>% 
+# pos_rate_year<-dat %>%
 #   dplyr::group_by(Year)%>%
-#   dplyr::summarise(No.Drifts=n(),Pos.Drifts=sum(ling.pa),Pos.Perc=round((sum(ling.pa)/n()),3)) %>% 
+#   dplyr::summarise(No.Drifts=n(),Pos.Drifts=sum(ling.pa),Pos.Perc=round((sum(ling.pa)/n()),3)) %>%
 #   as.data.frame()
 # pos_rate_year
 
-# # okay, I think that's it for filters. 
+# # okay, I think that's it for filters.
 
-# # exporting final dataset and positive rate table 
+# # exporting final dataset and positive rate table
 
 # #write.csv(dat,"Final_filtered_ORAtSea_updated.csv",row.names = F)
 # #write.csv(pos_rate_year,"PositivesByYear_updated.csv",row.names = F)
diff --git a/R/format_data.R b/R/format_data.R
index d9e0a0c..df5863a 100644
--- a/R/format_data.R
+++ b/R/format_data.R
@@ -46,7 +46,7 @@
 format_data <- function(data, ...) {
   if (all(
     c("Latitude_dd", "cpue_kg_per_ha_der", "Datetime_utc_iso") %in%
-    colnames(data)
+      colnames(data)
   )) {
     class(data) <- c("nwfscSurvey", class(data))
   }
diff --git a/R/format_strata.R b/R/format_strata.R
index 5783d57..dea58dc 100644
--- a/R/format_strata.R
+++ b/R/format_strata.R
@@ -9,7 +9,7 @@
 #' format_strata(nwfscSurvey::GetStrata.fn("sablefish"))
 format_strata <- function(strata = nwfscSurvey::GetStrata.fn("coast"),
                           min_depth = 55,
-                          max_depth = 1280){
+                          max_depth = 1280) {
   strata_temp <- strata %>%
     tidyr::separate(name, into = c("dname", "STRATA"), sep = "_") %>%
     dplyr::group_by(STRATA) %>%
diff --git a/R/lookup_grid.R b/R/lookup_grid.R
index fa783e1..701d704 100644
--- a/R/lookup_grid.R
+++ b/R/lookup_grid.R
@@ -50,7 +50,7 @@ lookup_grid <- function(x,
 
   out <- dplyr::mutate(
     .data = data,
-    area_km2 = {{column}},
+    area_km2 = {{ column }},
     vessel_year = "0",
     depth_scaled = scale(depth, center = mean_depth, scale = sd_depth),
     depth_scaled_squared = depth_scaled^2
diff --git a/R/map_base.R b/R/map_base.R
index ccd8c33..732f8b9 100644
--- a/R/map_base.R
+++ b/R/map_base.R
@@ -27,11 +27,11 @@ map_base <- function(x_min = -155,
   # Crop the polygon for plotting and efficiency:
   # sf::st_bbox(map_data) to find the rough coordinates
   coast <- suppressWarnings(suppressMessages(
-      sf::st_crop(
-        map_data,
-        c(xmin = x_min, ymin = y_min, xmax = x_max, ymax = y_max)
-      )
-    ))
+    sf::st_crop(
+      map_data,
+      c(xmin = x_min, ymin = y_min, xmax = x_max, ymax = y_max)
+    )
+  ))
   coast_proj <- sf::st_transform(coast, crs = utm_zone_10)
   gg <- ggplot2::ggplot(coast_proj) +
     ggplot2::geom_sf() +
diff --git a/R/map_density.R b/R/map_density.R
index d0216f1..3a98284 100644
--- a/R/map_density.R
+++ b/R/map_density.R
@@ -34,7 +34,6 @@ map_density <- function(predictions,
                         n_col = 2,
                         column_grep = "^est[2]*$",
                         tile_size = c(2000, 2000)) {
-
   column <- grep(column_grep, colnames(predictions), value = TRUE)
   predictions <- dplyr::rename(
     .data = predictions,
@@ -66,7 +65,7 @@ map_density <- function(predictions,
   )
   names(x) <- split_names
 
- gg <- map_base() +
+  gg <- map_base() +
     ggplot2::geom_tile(
       data = purrr::list_rbind(x, names_to = "year"),
       mapping = aes(x * 1000, y * 1000, fill = layer)
diff --git a/R/map_effects.R b/R/map_effects.R
index 44607f5..254336b 100644
--- a/R/map_effects.R
+++ b/R/map_effects.R
@@ -10,51 +10,53 @@
 #'
 #' @author Chantel R. Wetzel
 #' @export
-map_effects <- function(predictions, dir, verbose = FALSE){
-	column <- ifelse(
-		"est_non_rf2" %in% colnames(predictions),
-		"est_non_rf2",
-		ifelse("est_non_rf" %in% colnames(predictions),
-		"est_non_rf", 
-		"skip")
-	)
+map_effects <- function(predictions, dir, verbose = FALSE) {
+  column <- ifelse(
+    "est_non_rf2" %in% colnames(predictions),
+    "est_non_rf2",
+    ifelse("est_non_rf" %in% colnames(predictions),
+      "est_non_rf",
+      "skip"
+    )
+  )
 
-	if(column != 'skip'){
-	p1 <- map_nwfsc(predictions, exp(predictions[,column])) +
-  		scale_fill_viridis_c(trans = "sqrt", name = "Fixed \nEffects") +
-  		labs(x = "Longitude", y = "Latitude") + 
-  		ggtitle("Fixed pass effects")		
-	} else {
-		if(verbose){
-			message('The est_non_rf column not found in the predictions. Fixed effects map not created.')
-		}
-	}
+  if (column != "skip") {
+    p1 <- map_nwfsc(predictions, exp(predictions[, column])) +
+      scale_fill_viridis_c(trans = "sqrt", name = "Fixed \nEffects") +
+      labs(x = "Longitude", y = "Latitude") +
+      ggtitle("Fixed pass effects")
+  } else {
+    if (verbose) {
+      message("The est_non_rf column not found in the predictions. Fixed effects map not created.")
+    }
+  }
 
 
-  	column <- ifelse(
-		"omega_s2" %in% colnames(predictions),
-		"omega_s2",
-		ifelse("omega_s" %in% colnames(predictions),
-		"omega_s", 
-		"skip")
-	)
+  column <- ifelse(
+    "omega_s2" %in% colnames(predictions),
+    "omega_s2",
+    ifelse("omega_s" %in% colnames(predictions),
+      "omega_s",
+      "skip"
+    )
+  )
 
-  	if(column != 'skip'){
-  	p2 <- map_nwfsc(predictions, predictions[,column]) +
-  		scale_fill_viridis_c(name = "Random \nEffects") +
-  		labs(x = "Longitude", y = "Latitude") +
-  		ggtitle("Spatial random effects")  		
-  	} else {
-  		if(verbose){
-			message('The omega_s column not found in the predictions. Spatial random effects map not created.')
-		}
-  	}
+  if (column != "skip") {
+    p2 <- map_nwfsc(predictions, predictions[, column]) +
+      scale_fill_viridis_c(name = "Random \nEffects") +
+      labs(x = "Longitude", y = "Latitude") +
+      ggtitle("Spatial random effects")
+  } else {
+    if (verbose) {
+      message("The omega_s column not found in the predictions. Spatial random effects map not created.")
+    }
+  }
 
-  	g <- gridExtra::grid.arrange(p1, p2, nrow = 1)
+  g <- gridExtra::grid.arrange(p1, p2, nrow = 1)
 
-  	suppressMessages(ggsave(
-  		plot = g,
-    	filename = file.path(dir, paste0("fixed_and_spatial_effects.png")), 
-    	width = 14, height = 10, units = 'in'
-    ))
-}
\ No newline at end of file
+  suppressMessages(ggsave(
+    plot = g,
+    filename = file.path(dir, paste0("fixed_and_spatial_effects.png")),
+    width = 14, height = 10, units = "in"
+  ))
+}
diff --git a/R/map_residuals.R b/R/map_residuals.R
index 3ada6c6..7f09e98 100644
--- a/R/map_residuals.R
+++ b/R/map_residuals.R
@@ -14,7 +14,6 @@ map_residuals <- function(data,
                           n_row = 1,
                           n_col = 2,
                           save_prefix = file.path(getwd(), "data_")) {
-
   # data_extent <- raster::extent(data[, c("x", "y")])
   # data_raster <- raster::raster(data_extent,
   #   ncol = floor((slot(data_extent, "xmax") - slot(data_extent, "xmin")) / 2),
diff --git a/R/map_year_re.R b/R/map_year_re.R
index ec61b31..8d66e11 100644
--- a/R/map_year_re.R
+++ b/R/map_year_re.R
@@ -1,9 +1,9 @@
-#' Plot the annual random effects predictions by year 
+#' Plot the annual random effects predictions by year
 #' for the fit model
 #'
 #'
 #' @template predictions
-#' @template dir 
+#' @template dir
 #' @template nrow
 #' @template ncol
 #' @template verbose
@@ -14,64 +14,67 @@
 #'
 #' @examples
 #' \dontrun{
-#'		predictions <- predict(
-#'		    fit , 
-#'		    newdata = prediction_grid
-#'		) 
-#'		predictions[, c("lon", "lat")] <- round(predictions[, c("Lon", "Lat")], 1)
-#'		map_year_re(
-#'			predictions = predictions, 
-#'			dir = file.path(sppdir, survey, "index", obs)
-#'		)
-#'  }
+#' predictions <- predict(
+#'   fit,
+#'   newdata = prediction_grid
+#' )
+#' predictions[, c("lon", "lat")] <- round(predictions[, c("Lon", "Lat")], 1)
+#' map_year_re(
+#'   predictions = predictions,
+#'   dir = file.path(sppdir, survey, "index", obs)
+#' )
+#' }
 #'
-#'
-map_year_re <- function(predictions, dir, nrow = 3, ncol = 4, verbose = FALSE){
+map_year_re <- function(predictions, dir, nrow = 3, ncol = 4, verbose = FALSE) {
+  column <- ifelse(
+    "epsilon_st2" %in% colnames(predictions),
+    "epsilon_st2",
+    ifelse("epsilon_st" %in% colnames(predictions),
+      "epsilon_st",
+      "skip"
+    )
+  )
 
-	column <- ifelse(
-		"epsilon_st2" %in% colnames(predictions),
-		"epsilon_st2",
-		ifelse("epsilon_st" %in% colnames(predictions),
-		"epsilon_st", 
-		"skip")
-	)
+  if (column != "skip") {
+    num_years <- sort(unique(predictions$Year))
+    g <- split(
+      num_years,
+      ceiling(seq_along(num_years) / (ncol * nrow))
+    )
 
-	if(column != 'skip'){
-		num_years <- sort(unique(predictions$Year))
-		g <- split(
-			num_years, 
-			ceiling(seq_along(num_years) / (ncol * nrow))
-		)
-	
-		for(page in 1:length(g)) {
-  			map_nwfsc(predictions[predictions$Year %in% g[[page]], ], 
-  				predictions[predictions$Year %in% g[[page]], column]) +
-  				#scale_fill_gradient2() +
-  				scale_fill_viridis_c(
-  					  name = "Spatial \neffects",
-  					  # trim extreme high values to make spatial variation more visible
-  					  na.value = "yellow", 
-  					  limits = c(quantile(predictions[, column], 0.005), 
-  					  	quantile(predictions[, column], 0.995))
-  					) +
-  				facet_wrap(~Year, ncol = ncol, nrow = nrow) +
-  				labs(x = "Longitude", y = "Latitude")   + 
-  				ggtitle("Spatio-temporal random effects of the catch rate model")
-  			
-  			height <- ifelse(
-    			length(g[[page]]) == nrow * ncol, 10, 7)
+    for (page in 1:length(g)) {
+      map_nwfsc(
+        predictions[predictions$Year %in% g[[page]], ],
+        predictions[predictions$Year %in% g[[page]], column]
+      ) +
+        # scale_fill_gradient2() +
+        scale_fill_viridis_c(
+          name = "Spatial \neffects",
+          # trim extreme high values to make spatial variation more visible
+          na.value = "yellow",
+          limits = c(
+            quantile(predictions[, column], 0.005),
+            quantile(predictions[, column], 0.995)
+          )
+        ) +
+        facet_wrap(~Year, ncol = ncol, nrow = nrow) +
+        labs(x = "Longitude", y = "Latitude") +
+        ggtitle("Spatio-temporal random effects of the catch rate model")
 
-  			suppressMessages(ggsave(
-    			filename = file.path(dir, paste0("year_random_effects_", page, ".png")), 
-    			width = 10, 
-    			height = height, 
-    			units = 'in'
-    		))
-  		}
-  	} else {
-  		if(verbose){
-			message('The espsilon_st column not found in the predictions. Random effects map not created.')
-		}
-  	}
+      height <- ifelse(
+        length(g[[page]]) == nrow * ncol, 10, 7
+      )
 
+      suppressMessages(ggsave(
+        filename = file.path(dir, paste0("year_random_effects_", page, ".png")),
+        width = 10,
+        height = height,
+        units = "in"
+      ))
+    }
+  } else {
+    if (verbose) {
+      message("The espsilon_st column not found in the predictions. Random effects map not created.")
+    }
+  }
 }
diff --git a/R/plot_indices.R b/R/plot_indices.R
index a8512bb..fd65cca 100644
--- a/R/plot_indices.R
+++ b/R/plot_indices.R
@@ -49,9 +49,9 @@ plot_indices <- function(data,
     ) +
     ggplot2::theme_bw() +
     ggplot2::theme(
-    # If we want legend on the figure
-    #   legend.justification = c(0, 1),
-    #   legend.direction = "horizontal"
+      # If we want legend on the figure
+      #   legend.justification = c(0, 1),
+      #   legend.direction = "horizontal"
       legend.position = legend_loc,
     ) +
     ggplot2::scale_colour_viridis_d() +
diff --git a/R/plot_qq.R b/R/plot_qq.R
index ab3f088..a085ff8 100644
--- a/R/plot_qq.R
+++ b/R/plot_qq.R
@@ -75,12 +75,13 @@ plot_qq <- function(fit, file_name = "qq.png") {
     ggplot2::ylab("Sample quantiles") +
     ggplot2::labs(title = "Normal quantile-quantile plot") +
     ggplot2::theme_bw() +
-    ggplot2::theme(legend.position =
-      if (length(unique(data[["model"]])) > 1) {
-        c(0.15, 0.9)
-      } else {
-        "none"
-      }
+    ggplot2::theme(
+      legend.position =
+        if (length(unique(data[["model"]])) > 1) {
+          c(0.15, 0.9)
+        } else {
+          "none"
+        }
     ) +
     ggplot2::geom_abline(
       mapping = ggplot2::aes(
diff --git a/R/run.R b/R/run.R
index 78ddfcb..a0e2e4b 100644
--- a/R/run.R
+++ b/R/run.R
@@ -36,7 +36,7 @@ run <- function(data,
       "year", "fyear", "survey_name", "common_name",
       "catch_weight", "effort", "x", "y"
     ) %in%
-    colnames(data)
+      colnames(data)
   ))
 
   # Objects
diff --git a/R/run_sdmtmb.R b/R/run_sdmtmb.R
index 2e066b6..b4a7b51 100644
--- a/R/run_sdmtmb.R
+++ b/R/run_sdmtmb.R
@@ -32,7 +32,7 @@ run_sdmtmb <- function(dir_main = getwd(),
       "year", "fyear", "survey_name", "common_name",
       "catch_weight", "effort", "x", "y"
     ) %in%
-    colnames(data)
+      colnames(data)
   ))
   # TODO:
   # * check if random effects in formula and if converged, if yes and no, then
@@ -67,15 +67,15 @@ run_sdmtmb <- function(dir_main = getwd(),
 
   # Create prediction grid
   ranges <- data %>%
-      dplyr::filter(catch_weight > 0) %>%
-      dplyr::summarize(
-        dplyr::across(
-          dplyr::matches("tude"),
-          .fns = list("max" = ~ max(.) + 0.1, "min" = ~ min(.) - 0.1)
-        ),
-        # depth_min = min(abs(depth), na.rm = TRUE),
-        depth_max = min(depth, na.rm = TRUE)
-      )
+    dplyr::filter(catch_weight > 0) %>%
+    dplyr::summarize(
+      dplyr::across(
+        dplyr::matches("tude"),
+        .fns = list("max" = ~ max(.) + 0.1, "min" = ~ min(.) - 0.1)
+      ),
+      # depth_min = min(abs(depth), na.rm = TRUE),
+      depth_max = min(depth, na.rm = TRUE)
+    )
   data_truncated <- data %>%
     dplyr::filter(
       latitude > ranges[["latitude_min"]] & latitude < ranges[["latitude_max"]],
diff --git a/data-raw/configuration.R b/data-raw/configuration.R
index e5d8848..a144e15 100644
--- a/data-raw/configuration.R
+++ b/data-raw/configuration.R
@@ -54,21 +54,21 @@ indices <- dir(
       load(x)
       if (exists("index_areas")) {
         return(data.frame(i = x, index_areas))
-        } else {
-          return(NULL)
-        }
+      } else {
+        return(NULL)
+      }
     }
   )
 gg <- ggplot2::ggplot(
- data = indices %>%
-  dplyr::mutate(
-    dist = basename(dirname(dirname(i)))
-  ) %>%
-  dplyr::filter(area == "coastwide"),
+  data = indices %>%
+    dplyr::mutate(
+      dist = basename(dirname(dirname(i)))
+    ) %>%
+    dplyr::filter(area == "coastwide"),
   ggplot2::aes(x = year, y = est, lty = dist, col = dist, group = i)
 ) +
-ggplot2::geom_line() +
-ggplot2::theme_bw()
+  ggplot2::geom_line() +
+  ggplot2::theme_bw()
 ggsave(gg, filename = "indexwc_copper_rockfish.png")
 
 
diff --git a/data-raw/grid_nwfscSurvey.R b/data-raw/grid_nwfscSurvey.R
index 65d7e4e..23ef8f9 100644
--- a/data-raw/grid_nwfscSurvey.R
+++ b/data-raw/grid_nwfscSurvey.R
@@ -26,12 +26,12 @@ california_current_grid <- california_current_grid %>%
     )
   )
 
-  california_current_grid <- suppressWarnings(sdmTMB::add_utm_columns(
-    california_current_grid,
-    utm_crs = utm_zone_10,
-    ll_names = c("Lon", "Lat"),
-    utm_names = c("x", "y")
-  )) %>%
+california_current_grid <- suppressWarnings(sdmTMB::add_utm_columns(
+  california_current_grid,
+  utm_crs = utm_zone_10,
+  ll_names = c("Lon", "Lat"),
+  utm_names = c("x", "y")
+)) %>%
   dplyr::rename(
     longitude = "Lon",
     latitude = "Lat"
diff --git a/man/calc_index_areas.Rd b/man/calc_index_areas.Rd
index 5990c5e..082a1f4 100644
--- a/man/calc_index_areas.Rd
+++ b/man/calc_index_areas.Rd
@@ -62,7 +62,8 @@ boundaries_data[4]
 list("test" = c(39, 38))
 # run the function
 index <- calc_index_areas(
-  data, fit, grid, dir = getwd(), boundaries = list("test" = c(39, 38))
+  data, fit, grid,
+  dir = getwd(), boundaries = list("test" = c(39, 38))
 )
 # look at the index, which will be in "test" because that is what we
 # named the boundary
diff --git a/man/map_year_re.Rd b/man/map_year_re.Rd
index 858d344..050efad 100644
--- a/man/map_year_re.Rd
+++ b/man/map_year_re.Rd
@@ -27,17 +27,16 @@ for the fit model
 }
 \examples{
 \dontrun{
-	predictions <- predict(
-	    fit , 
-	    newdata = prediction_grid
-	) 
-	predictions[, c("lon", "lat")] <- round(predictions[, c("Lon", "Lat")], 1)
-	map_year_re(
-		predictions = predictions, 
-		dir = file.path(sppdir, survey, "index", obs)
-	)
- }
-
+predictions <- predict(
+  fit,
+  newdata = prediction_grid
+)
+predictions[, c("lon", "lat")] <- round(predictions[, c("Lon", "Lat")], 1)
+map_year_re(
+  predictions = predictions,
+  dir = file.path(sppdir, survey, "index", obs)
+)
+}
 
 }
 \author{