AFSC Eastern Bering Sea using wrapper_surveyindex()

# get rid of memory limits
options(future.globals.maxSize = 1 * 1024^4) # Allow up to 1 TB for globals

# Set directories --------------------------------------------------------------
library(here)

wd <- paste0(here::here(), "/vignettes/")
dir_out <- paste0(wd, "output/")
crs_latlon <- "+proj=longlat +datum=WGS84" # decimal degrees

# Install Libraries ------------------------------------------------------------

# Here we list all the packages we will need for this whole process
# We'll also use this in our works cited page.
PKG <- c(
  "surveyresamplr",

  # tidyverse
  "dplyr",
  "tidyr",
  "viridis",
  "ggplot2",
  "tibble",
  "janitor",
  "data.table",

  # parallelizing
  "forcats",
  "purrr",
  "furrr",
  "doParallel",

  # sampling
  "sampling",

  # modeling
  "arrow",
  "future.apply",
  "future.callr",
  "sdmTMB", # install.packages("remotes"",; remotes::install_github("pbs-assess/sdmTMBextra", dependencies = TRUE",
  "DTUAqua", # remotes::install_github("DTUAqua/DATRAS/DATRAS")
  "surveyIndex", # remotes::install_github("casperwberg/surveyIndex/surveyIndex")
  "Matrix",
  "MASS",
  "cluster",
  "TMB",
  "INLA"
)

pkg_install <- function(p) {
  if (!require(p, character.only = TRUE)) {
    install.packages(p)
  }
  require(p, character.only = TRUE)
}
base::lapply(unique(PKG), pkg_install)

### Define study species -------------------------------------------------------

spp_list <- data.frame(
  srvy = "EBS",
  common_name = c(
    "walleye pollock", "snow crab", "Pacific cod",
    "red king crab", "blue king crab",
    "yellowfin sole", "Pacific halibut",
    "Alaska plaice", "flathead sole", "northern rock sole", "arrowtooth flounder"
  ),
  species_code = as.character(c(
    21740, 68580, 21720,
    69322, 69323,
    10210, 10120,
    10285, 10130, 10261, 10110
  )),
  filter_lat_lt = NA,
  filter_lat_gt = NA,
  filter_depth = NA,
  model_fn = "total_catch_wt_kg ~ 0 + factor(year) + bottom_temperature_c",
  model_family = "delta_gamma",
  model_anisotropy = TRUE,
  model_spatiotemporal = "iid, iid"
) |>
  dplyr::mutate(
    file_name = gsub(pattern = " ", replacement = "_", x = (tolower(common_name)))
  )

### Load survey data -----------------------------------------------------------

# source(paste0(wd, "code/data_dl_ak.r"))

load(file = paste0(wd, "/data/noaa_afsc_catch.rda"))
catch <- noaa_afsc_catch |> dplyr::filter(srvy == "EBS")

### Load grid data -------------------------------------------------------------

load(paste0(wd, "grids/noaa_afsc_ebs_pred_grid_depth.rdata"), verbose = TRUE)

#### Add temperature: Coldpool temperature data
# Data that varies over space and time (bottom temperature)
# Here, bottom temperature, and thereby the cold pool extent, have been show to drive the distribution of many species. This is especially true for walleye pollock.
# For this we are going to lean on our in-house prepared validated and pre-prepared [{coldpool} R package](https://github.com/afsc-gap-products/coldpool) (S. Rohan, L. Barnett, and N. Charriere). This data interpolates over the whole area of the survey so there are no missing data.
grid_yrs <-
  dplyr::bind_cols(
    pred_grid_depth[, c("longitude_dd", "latitude_dd", "depth_m")],
    terra::unwrap(coldpool::ebs_bottom_temperature) |>
      terra::project(crs_latlon) |>
      terra::extract(pred_grid_depth[, c("longitude_dd", "latitude_dd")])
  )
grid_yrs <- grid_yrs |>
  tidyr::pivot_longer(
    names_to = "year",
    values_to = "bottom_temperature_c",
    cols = names(grid_yrs_temperature)[4:ncol(grid_yrs_temperature)]
  )
save(grid_yrs_depth_temperature, file = paste0("grids/grid_yr_temperature/noaa_afsc_ebs_pred_grid_depth_temperature.rdata"))

# # test you extracted correctkt
# ggplot(data = grid_yrs |>
#          dplyr::filter(year %in% c(2022:2024)),
#        mapping = ggplot2::aes(x = longitude_dd, y = latitude_dd, color = bottom_temperature_c)) +
#   geom_point() +
#   facet_wrap(facets = "year")

### Variables ------------------------------------------------------------------

srvy <- "EBS"
seq_from <- 0.2
seq_to <- 1.0
seq_by <- 0.2
tot_dataframes <- 13
replicate_num <- 3

### Run ------------------------------------------------------------------------

purrr::map(
  1:nrow(spp_list),
  ~ clean_and_resample(spp_list[.x, ],
    catch, seq_from, seq_to, seq_by,
    tot_dataframes, replicate_num, grid_yrs, dir_out,
    test = TRUE
  )
)

### Plot indices ---------------------------------------------------------------

plot_results(srvy = srvy, dir_out = dir_out)