AFSC Eastern and Northern Bering Sea

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

# 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",
  "dplyr",
  "tidyr",
  "purrr",
  "ggplot2",
  "future.callr",
  "flextable", 
  "sdmTMB"
)

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

# Set directories --------------------------------------------------------------
wd <- paste0(here::here(), "/vignettes/")
dir_out <- paste0(wd, "output/")
dir_final <- paste0(dir_out, "BS_0results/")

Run scenarios

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

spp_list <- data.frame(
  srvy = "BS",
  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 + depth_m",
  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"))
catch <- surveyresamplr::noaa_afsc_catch |>
  dplyr::filter(srvy %in% c("NBS", "EBS"))

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

# grid_yrs <- replicate_df(surveyresamplr::noaa_afsc_GOA_pred_grid_depth, "year", unique(catch$year))
pred_grid_depth <- surveyresamplr::noaa_afsc_bs_pred_grid_depth
#### 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.
crs_latlon <- "+proj=longlat +datum=WGS84" # decimal degrees
ebs_only <- setdiff(names(terra::unwrap(coldpool::ebs_bottom_temperature)), names(terra::unwrap(coldpool::nbs_ebs_bottom_temperature)))
grid_yrs_temperature <- dplyr::full_join(
  dplyr::bind_cols(
    pred_grid_depth[, c("longitude_dd", "latitude_dd", "depth_m", "area_km2")],
    terra::unwrap(coldpool::ebs_bottom_temperature) |> # TOLEDO
      terra::subset(ebs_only) |>
      terra::project(crs_latlon) |>
      terra::extract(pred_grid_depth[, c("longitude_dd", "latitude_dd")])
  ),
  dplyr::bind_cols(
    pred_grid_depth[, c("longitude_dd", "latitude_dd", "depth_m")],
    terra::unwrap(coldpool::nbs_ebs_bottom_temperature) |>
      terra::project(crs_latlon) |>
      terra::extract(pred_grid_depth[, c("longitude_dd", "latitude_dd")])
  )
)

grid_yrs_depth_temperature <- grid_yrs <- grid_yrs_temperature |>
  tidyr::pivot_longer(
    names_to = "year",
    values_to = "bottom_temperature_c",
    cols = names(grid_yrs_temperature)[6:ncol(grid_yrs_temperature)]
  ) |>
  dplyr::mutate(year = as.numeric(year))
# save(grid_yrs_depth_temperature, file = paste0("grids/grid_yr_temperature/noaa_afsc_bs_pred_grid_depth_temperature.rdata"))

# # test you extracted correct
# 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 <- "BS"
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 --------------------------------------------------------------

out <- plot_results(srvy = srvy, dir_out = dir_out, dir_final = dir_final)

out$plots

Parameter output:

i <- 1
print(names(out$tables)[i])
head(out$tables[i][[1]]) |>
  flextable::flextable()

i <- 1 + i
print(names(out$tables)[i])
head(out$tables[i][[1]]) |>
  flextable::flextable()

i <- 1 + i
print(names(out$tables)[i])
head(out$tables[i][[1]]) |>
  flextable::flextable()

i <- 1 + i
print(names(out$tables)[i])
head(out$tables[i][[1]]) |>
  flextable::flextable()