In pbs-assess/gfranges: Groundfish Ranges

knitr::opts_chunk$set(
  fig.width = 11, fig.height = 8.5,
  # fig.path=paste0("figs/", spp, "/"),
  echo = FALSE, warning = FALSE, message = FALSE
)
library(dplyr)
library(ggplot2)
library(gfplot)
library(gfdata)
library(sdmTMB)
library(gfranges)

species <- params$species
region <- params$region
covariates <- params$covariates
covs <- params$covs
priors <- FALSE

paste("region =", region)
paste("model covariates =", covariates)
paste("model label =", covs)
paste("priors =", priors)

spp <- gsub(" ", "-", gsub("\\/", "-", tolower(species)))

# folder to hold figs for this species
dir.create(file.path("figs", spp))
dir.create(file.path("data", spp))

if (region == "Both odd year surveys") {
  survey <- c("SYN QCS", "SYN HS")
  model_ssid <- c(1, 3)
  ssid_string <- paste0(model_ssid, collapse = "n")
  years <- NULL
}

if (region == "West Coast Vancouver Island") {
  survey <- c("SYN WCVI")
  model_ssid <- c(4)
  ssid_string <- paste0(model_ssid, collapse = "n")
  years <- NULL
}

if (region == "West Coast Haida Gwaii") {
  survey <- c("SYN WCHG")
  model_ssid <- c(16)
  ssid_string <- paste0(model_ssid, collapse = "n")
  years <- NULL
}

if (region == "All synoptic surveys") {
  survey <- c("SYN QCS", "SYN HS", "SYN WCVI", "SYN WCHG")
  model_ssid <- c(1, 3, 4, 16)
  ssid_string <- paste0(model_ssid, collapse = "n")
  years <- NULL
}

Build spatiotemporal density model

Load and filter data

biomass <- readRDS(paste0(
  "data/", spp, "/data-by-maturity-", spp, "-", ssid_string, ".rds"
))

if(nrow(biomass)<4000) {stop("Need to recalculate split by maturity!")}

covars <- readRDS("data/event-covariates.rds")
data <- dplyr::left_join(biomass, covars)

# scale predictors before filtering to ensure mean and SD are global
data <- data %>% mutate(raw_depth = depth, depth = log(raw_depth))
data <- gfranges::scale_predictors(data, # predictors = c(quo(depth)))
  predictors = c(quo(depth), quo(mixed), quo(muddy), quo(sandy), quo(rocky), quo(any_rock))
)
data <- data %>%
  mutate(depth = raw_depth) %>%
  filter(ssid %in% model_ssid) %>%
  filter(year < 2019)

Model without environmental covariates

max_depth_found <- max(data[data$present == 1, ]$raw_depth, na.rm = TRUE)
rm(depth_model_list)
rm(adult_biomass)
rm(imm_biomass)
rm(total_biomass)
rm(depth_plots)
#max_depth_found <- 800

try({
  adult_biomass <- readRDS(paste0(
    "data/", spp,
    "/mod-mat-biomass-", spp, covs, "-", ssid_string, "-prior-", priors, ".rds"
  ))
})
# adult_biomass<-sdmTMB:::update_model(adult_biomass)

try({
  imm_biomass <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, covs, "-1n3n4n16-prior-FALSE.rds"
  ))
# imm_biomass<-sdmTMB:::update_model(imm_biomass)
  # depth_model_list <- list(adult = adult_biomass, imm = imm_biomass)
})

if (!exists("adult_biomass")) {
  try({
    total_biomass <- readRDS(paste0(
      "data/", spp,
      "/model-total-biomass-", spp, covs, "-", ssid_string, "-prior-", priors, ".rds"
    ))
  })
}


if (!exists("adult_biomass")) adult_biomass <- NULL
if (!exists("imm_biomass")) imm_biomass <- NULL
if (!exists("total_biomass")) total_biomass <- NULL

depth_model_list <- list(adult = adult_biomass, imm = imm_biomass, total = total_biomass)
depth_model_list <- depth_model_list[!sapply(depth_model_list, is.null)]

depth_model_list

write_grads <- function(x, macro, ...) {
  paste0("\\newcommand{\\", macro, "}{", x, "}") %>%
    readr::write_lines("model-grads.tex", append = TRUE)
}

write_grads(signif(max(depth_model_list$adult$gradients), digits = 8), paste0(spp, "-mat"))
write_grads(signif(max(depth_model_list$imm$gradients), digits = 8), paste0(spp, "-imm"))
write_grads(signif(max(depth_model_list$total$gradients), digits = 8), paste0(spp, "-total"))

try({
  if(max(depth_model_list$adult$gradients)>0.001) {
    adult_biomass <- sdmTMB:::run_extra_optimization(depth_model_list$adult, 
      nlminb_loops = 2, newton_steps = 2)
    write_grads(signif(max(adult_biomass$gradients), digits = 8), paste0(spp, "-mat2"))
  }
})

try({
  if(max(depth_model_list$imm$gradients)>0.001) {
    imm_biomass <- sdmTMB::run_extra_optimization(depth_model_list$imm, 
      nlminb_loops = 2, newton_steps = 2)
    write_grads(signif(max(imm_biomass$gradients), digits = 8), paste0(spp, "-imm2"))
  }
})

try({
  if(max(depth_model_list$total$gradients)>0.001) {
    total_biomass <- sdmTMB::run_extra_optimization(depth_model_list$total, 
      nlminb_loops = 2, newton_steps = 2)
    write_grads(signif(max(total_biomass$gradients), digits = 8), paste0(spp, "-total2"))
  }
})

 if (any(names(data) == "adult_density")) {

    saveRDS(adult_biomass,
      file = paste0(
        "data/", spp,
        "/mod-mat-biomass-", spp, covs, "-", ssid_string, "-prior-", priors, ".rds"
      )
    )

    saveRDS(imm_biomass,
      file = paste0("data/", spp, "/mod-imm-biomass-", 
        spp, covs, "-", ssid_string, "-prior-", priors, ".rds"))


  } else {
    saveRDS(total_biomass, file = paste0(
      "data/", spp,
      "/model-total-biomass-", spp, covs, "-", ssid_string, "-prior-", priors, ".rds"
    ))
  }