pbs-assess/gfranges
Groundfish Ranges

Package index
Search the pbs-assess/gfranges package
Vignettes
Functions
109
Source code
82
Man pages
47
Home
/
GitHub
/
pbs-assess/gfranges
/

In pbs-assess/gfranges: Groundfish Ranges 

knitr::opts_chunk$set(
  fig.width = 9, fig.height = 6.5, # 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
start_year <- params$start_year

bind_predictions <- params$bind_predictions
covs <- params$covs
priors <- FALSE



paste("species =", species)
paste("region =", region)
paste("ssids modelled separately =", bind_predictions)
paste("model label =", covs)
paste("priors =", priors)

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

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
}

Save predictions for full spatial grid

# fish <- readRDS(paste0("raw/bio-data-", spp, ""))

rm(adult_predictions)
rm(imm_predictions)
rm(predictions)

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

if (any(names(maturity) == "adult_density")) {
  # if (length(unique(fish$maturity_code)) > 2) {
  nd_all <- readRDS(paste0("data/predicted-DO-new.rds"))

  adult_biomass <- readRDS(paste0(
    "data/", spp,
    "/mod-mat-biomass-", spp, covs, "-all--", ssid_string, "-prior-", priors, ".rds"
  ))
  # nd$year <- as.integer(nd$year)

  sort(unique(nd_all$year))
  nd <- nd_all %>%
    filter(ssid %in% model_ssid) %>%
    filter(year %in% unique(adult_biomass$data$year)) %>%
    mutate(
      do_mlpl_scaled = (log_do - adult_biomass$data$do_mlpl_mean[1]) / (adult_biomass$data$do_mlpl_sd[1] * 2),
      do_mlpl_scaled2 = do_mlpl_scaled^2
    ) %>%
    select(-est, -est_non_rf, -est_rf, -omega_s, -epsilon_st, -zeta_s)
  nd <- na.omit(nd)
  sort(unique(nd$year))

  adult_predictions_obj <- predict(adult_biomass, newdata = nd, return_tmb_object = TRUE)

  # multiyear bioclimatic calculation
  model_pre_2015 <- readRDS(paste0(
    "data/", spp,
    "/mod-mat-biomass-", spp, covs, "-past-2015-", ssid_string, "-prior-", priors, ".rds"
  ))

  bioclim_blob <- bioclim_multi_year(2015, model_pre_2015, nd)
  bioclim_blob <- rename(bioclim_blob, bioclim_blob = bioclim)
  glimpse(bioclim_blob)
  #View(bioclim_blob)

  adult_predictions <- adult_predictions_obj$data
  adult_predictions <- left_join(adult_predictions, bioclim_blob)

  # saveRDS(adult_predictions, file = paste0(
  #   "data/", spp,
  #   "/bioclimatic-predictions-", spp, covs, "-", ssid_string, "-mat-2015-cutoff.rds"
  # ))
  # 
  # add biennial bioclimatic calculation
  model_pre_2011 <- readRDS(paste0(
    "data/", spp,
    "/mod-mat-biomass-", spp, covs, "-past-2011-", ssid_string, "-prior-", priors, ".rds"
  ))
  sort(unique(model_pre_2011$data$year))

  bioclim_2011 <- bioclim_by_year(2011, model_pre_2011, nd)
  bioclim_2011 <- rename(bioclim_2011, bioclim_2011 = bioclim)
  glimpse(bioclim_2011)

  model_pre_2013 <- readRDS(paste0(
    "data/", spp,
    "/mod-mat-biomass-", spp, covs, "-past-2013-", ssid_string, "-prior-", priors, ".rds"
  ))
  sort(unique(model_pre_2013$data$year))

  bioclim_2013 <- bioclim_by_year(2013, model_pre_2013, nd)
  bioclim_2013 <- rename(bioclim_2013, bioclim_2013 = bioclim)
  glimpse(bioclim_2013)

  model_pre_2015 <- readRDS(paste0(
    "data/", spp,
    "/mod-mat-biomass-", spp, covs, "-past-2015-", ssid_string, "-prior-", priors, ".rds"
  ))
  sort(unique(model_pre_2015$data$year))

  bioclim_2015 <- bioclim_by_year(2015, model_pre_2015, nd)
  bioclim_2015 <- rename(bioclim_2015, bioclim_2015 = bioclim)
  glimpse(bioclim_2015)

  model_pre_2017 <- readRDS(paste0(
    "data/", spp,
    "/mod-mat-biomass-", spp, covs, "-past-2017-", ssid_string, "-prior-", priors, ".rds"
  ))
  sort(unique(model_pre_2017$data$year))

  bioclim_2017 <- bioclim_by_year(2017, model_pre_2017, nd)
  bioclim_2017 <- rename(bioclim_2017, bioclim_2017 = bioclim)
  glimpse(bioclim_2017)

  adult_predictions <- left_join(adult_predictions, bioclim_2011)
  adult_predictions <- left_join(adult_predictions, bioclim_2013)
  adult_predictions <- left_join(adult_predictions, bioclim_2015)
  adult_predictions <- left_join(adult_predictions, bioclim_2017)

#  View(adult_predictions)

  saveRDS(adult_predictions, file = paste0(
    "data/", spp,
    "/bioclimatic-predictions-", spp, covs, "-", ssid_string, "-mat-2015-cutoff.rds"
  ))


  #################
  #### Immature
  #################

  imm_biomass <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, covs, "-all--", ssid_string, "-prior-", priors, ".rds"
  ))

  nd <- nd_all %>%
    filter(ssid %in% model_ssid) %>%
    filter(year %in% unique(imm_biomass$data$year)) %>%
    mutate(do_mlpl_scaled = (log_do - imm_biomass$data$do_mlpl_mean[1]) / (imm_biomass$data$do_mlpl_sd[1] * 2), do_mlpl_scaled2 = do_mlpl_scaled^2)
  nd <- na.omit(nd)

  imm_predictions_obj <- predict(imm_biomass, newdata = nd, return_tmb_object = TRUE)
  imm_predictions <- imm_predictions_obj$data

  imodel_pre_2013 <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, covs, "-past-2013-", ssid_string, "-prior-", priors, ".rds"
  ))
  sort(unique(imodel_pre_2013$data$year))

  bioclim_2013 <- bioclim_by_year(2013, imodel_pre_2013, nd)
  bioclim_2013 <- rename(bioclim_2013, bioclim_2013 = bioclim)
  glimpse(bioclim_2013)

  imodel_pre_2015 <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, covs, "-past-2015-", ssid_string, "-prior-", priors, ".rds"
  ))
  sort(unique(imodel_pre_2015$data$year))

  bioclim_2015 <- bioclim_by_year(2015, imodel_pre_2015, nd)
  bioclim_2015 <- rename(bioclim_2015, bioclim_2015 = bioclim)
  glimpse(bioclim_2015)

  imodel_pre_2017 <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, covs, "-past-2017-", ssid_string, "-prior-", priors, ".rds"
  ))
  sort(unique(imodel_pre_2017$data$year))

  bioclim_2017 <- bioclim_by_year(2017, imodel_pre_2017, nd)
  bioclim_2017 <- rename(bioclim_2017, bioclim_2017 = bioclim)
  glimpse(bioclim_2017)

  imodel_pre_2015 <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, covs, "-past-2015-", ssid_string, "-prior-", priors, ".rds"
  ))

  bioclim_blob <- bioclim_multi_year(2015, imodel_pre_2015, nd)
  bioclim_blob <- rename(bioclim_blob, bioclim_blob = bioclim)
  glimpse(bioclim_blob)
  #View(bioclim_blob)

  imm_predictions <- left_join(imm_predictions, bioclim_blob)
  imm_predictions <- left_join(imm_predictions, bioclim_2013)
  imm_predictions <- left_join(imm_predictions, bioclim_2015)
  imm_predictions <- left_join(imm_predictions, bioclim_2017)

  saveRDS(imm_predictions, file = paste0(
    "data/", spp,
    "/bioclimatic-predictions-", spp, covs, "-", ssid_string, "-imm-2015-cutoff.rds"
  ))

} else { # need to add non-split-by-age version
  # total_biomass <- readRDS(paste0("data/", spp,
  #   "/model-total-biomass-", spp, "-fixed", covs, "-", ssid_string, "-prior-", priors, ".rds"))
  #
  # nd_all <- readRDS(paste0("data/predicted-DO-new.rds"))
  # nd <- nd_all %>%
  #   filter(ssid %in% model_ssid) %>%
  #   filter(year %in% unique(total_biomass$data$year))%>%
  # mutate(do_mlpl_scaled = (log_do - total_biomass$data$do_mlpl_mean) /(total_biomass$data$do_mlpl_sd*2), do_mlpl_scaled2 = do_mlpl_scaled^2)
  # nd <- na.omit(nd)
  # nd$year <- as.integer(nd$year)
  #
  # predictions <- predict(total_biomass, newdata = nd)
  # saveRDS(predictions, file = paste0("data/", spp,
  #   "/predictions-", spp, covs, "-", ssid_string, "-total-bioclimatic-prior-", priors, ".rds"))
}

Load model predictions

rm(adult_predictions)
rm(imm_predictions)
rm(predictions)

rm(pred1n3)
rm(pred4)
rm(pred16)

rm(pred1n3_imm)
rm(pred4_imm)
rm(pred16_imm)

rm(pred1n3_all)
rm(pred4_all)
rm(pred16_all)

if (region != "All synoptic surveys") {
  try(adult_predictions <- readRDS(paste0(
    "data/", spp,
    "/bioclimatic-predictions-", spp, covs, "-", ssid_string, "-mat-2015-cutoff.rds"
  )) %>% filter(ssid %in% model_ssid))
  try(imm_predictions <- readRDS(paste0(
    "data/", spp,
    "/bioclimatic-predictions-", spp, covs, "-", ssid_string, "-imm-2015-cutoff.rds"
  )) %>% filter(ssid %in% model_ssid))

  if (!exists("adult_predictions")) {
    paste("Maturity split biomass predictions are not available for this survey.")

    try(predictions <- readRDS(paste0(
      "data/", spp,
      "/predictions-", spp, covs, "-1n3n4n16-total-bioclimatic-prior-", priors, ".rds"
    )) %>% filter(ssid %in% model_ssid))

    if (!exists("predictions")) {
      paste("Total biomass predictions for this survey are also not available.")
    }
  }
} else {
  try({
    adult_predictions <- readRDS(paste0(
      "data/", spp,
      "/bioclimatic-predictions-", spp, covs, "-", ssid_string, "-mat-2015-cutoff.rds"
    )) %>% filter(ssid %in% model_ssid)
  })
  try({
    imm_predictions <- readRDS(paste0(
      "data/", spp,
      "/bioclimatic-predictions-", spp, covs, "-", ssid_string, "-imm-2015-cutoff.rds"
    )) %>% filter(ssid %in% model_ssid)
  })

  if (!exists("adult_predictions")) {
    paste("Maturity split biomass predictions are not available.")

    try({
      predictions <- readRDS(paste0(
        "data/", spp,
        "/predictions-", spp, covs, "-", ssid_string, "-total-bioclimatic-prior-", priors, ".rds"
      ))
    })

    if (!exists("predictions")) paste("Total biomass predictions are not available.")
  }
}

Caluclate variables of interest from predictions

rm(max_raster)
rm(max_adult)
rm(max_bio)

if (exists("adult_predictions")) {
  # to convert from kg/m2 to kg/hectare multiply by 10000
  adult_predictions$est_exp <- exp(adult_predictions$est) * 10000
  imm_predictions$est_exp <- exp(imm_predictions$est) * 10000
  adult_predictions$total_bio <- imm_predictions$est_exp + adult_predictions$est_exp

  # check that both dataframes are arranged similarly
  if (adult_predictions[30000, 6] == imm_predictions[30000, 6]) {
    imm_predictions$prop_imm <-
      imm_predictions$est_exp / (adult_predictions$est_exp + imm_predictions$est_exp)
  }

  saveRDS(adult_predictions, file = paste0(
    "data/", spp,
    "/all-predictions-", spp, covs, "-mat-bioclimatic-prior-", priors, ".rds"
  ))
  saveRDS(imm_predictions, file = paste0(
    "data/", spp,
    "/all-predictions-", spp, covs, "-imm-bioclimatic-prior-", priors, ".rds"
  ))

  adult_predictions <- filter(adult_predictions, year >= start_year)
  imm_predictions <- filter(imm_predictions, year >= start_year)

  max_raster <- quantile(adult_predictions$est_exp, 0.999)
  max_adult <- signif(max(adult_predictions$est_exp), digits = 2)

  model_ssid <- unique(adult_predictions$ssid)
} else {
  predictions <- filter(predictions, year >= start_year)

  predictions$est_exp <- exp(predictions$est) * 10000

  # predictions <- predictions %>% group_by(X,Y) %>% arrange(year) %>% mutate(lag_epsilon = lag(epsilon_st), bioclimatic = est - epsilon_st + lag_epsilon)
  #
  predictions$total_bio <- exp(predictions$est) * 10000
  max_raster <- quantile(predictions$total_bio, .999)
  # max_bio <- round(max(predictions$total_bio))
  max_bio <- signif(max(predictions$est_exp), digits = 2)

  saveRDS(predictions, file = paste0(
    "data/", spp,
    "/all-predictions-", spp, covs, "-total-bioclimatic-prior-", priors, ".rds"
  ))

  model_ssid <- unique(predictions$ssid)
}

Mature (or total) biomass predictions

if (exists("adult_predictions")) {
  p_adult_all <- plot_facet_map(adult_predictions, "est_exp",
    raster_limits = c(0, max_raster),
    legend_position = "none",
    transform_col = fourth_root_power
  ) +
    labs(fill = "kg/ha") +
    ggtitle(paste0("", species, " mature biomass \n(max = ", max_adult, " kg/ha)"))
  print(p_adult_all)
} else {
  p_adult_all <- plot_facet_map(predictions, "total_bio",
    transform_col = fourth_root_power
  ) +
    labs(fill = "kg/ha") +
    ggtitle(paste0("", species, " total biomass \n(max = ", max_bio, " kg/ha)"))

  print(p_adult_all)
}

Split by survey (SSID)

if (exists("adult_predictions")) predictions <- adult_predictions

if (1 %in% model_ssid) {
  adult_predictions1 <- predictions %>% filter(ssid == 1)
} else {
  adult_predictions1 <- NULL
}

if (3 %in% model_ssid) {
  adult_predictions3 <- predictions %>%
    filter(ssid != 16) %>%
    filter(ssid != 4)
} else {
  adult_predictions3 <- NULL
}

if (4 %in% model_ssid) {
  adult_predictions4 <- predictions %>% filter(ssid == 4)
} else {
  adult_predictions4 <- NULL
}

if (16 %in% model_ssid) {
  if (any(predictions$ssid == 16)) {
    adult_predictions16 <- predictions %>% filter(ssid == 16)
  } else {
    adult_predictions16 <- NULL
  }
} else {
  adult_predictions16 <- NULL
}

ad_prediction_list <- list(adult_predictions1, adult_predictions3, adult_predictions4, adult_predictions16)
ad_prediction_list <- ad_prediction_list[!sapply(ad_prediction_list, is.null)]

p_adult <- list()
if (exists("adult_predictions")) {
  for (i in seq_len(length(model_ssid))) {
    p_adult[[i]] <- plot_facet_map(ad_prediction_list[[i]], "est_exp",
      raster_limits = c(0, max_raster),
      # legend_position = "none"
      transform_col = fourth_root_power
    ) +
      labs(fill = "kg/ha") +
      ggtitle(paste0("", species, " mature biomass \n(max = ", max_adult, " kg/ha)"))
    p_adult[[i]]
  }
} else {
  for (i in seq_len(length(model_ssid))) {
    p_adult[[i]] <- plot_facet_map(ad_prediction_list[[i]], "est_exp",
      raster_limits = c(0, max_raster),
      # legend_position = "none"
      transform_col = fourth_root_power
    ) +
      labs(fill = "kg/ha") +
      ggtitle(paste0("", species, " total biomass \n(max = ", max_bio, " kg/ha)"))
    p_adult[[i]]
  }
}
print(p_adult)

rm(out_mature)
rm(out_total)
# try(out_mature <- readRDS(paste0("data/", spp, "/trends-", spp, covs, "-", ssid_string, "-mature-bioclimatic-priors-", priors, ".rds")))

if (!exists("out_mature")) {
  if (exists("adult_predictions")) {
    predictions <- adult_predictions
  } else {
    try(out_total <- readRDS(paste0("data/", spp, "/trends-", spp, covs, "-", ssid_string, "-total-bioclimatic-priors-", priors, ".rds")))
  }

  if (!exists("out_total")) {
    if (1 %in% model_ssid) {
      .predictions <- filter(predictions, ssid == 1)
      if (length(unique(.predictions$year)) == 5) indices <- c(1, 1, 1, 2, 2)
      if (length(unique(.predictions$year)) == 6) indices <- c(1, 1, 1, 2, 2, 2)
      if (length(unique(.predictions$year)) == 7) indices <- c(1, 1, 1, 2, 2, 2, 2)
      if (length(unique(.predictions$year)) == 8) indices <- c(1, 1, 1, 1, 2, 2, 2, 2)

      out_mature1 <- make_trend_data(predictions,
        ssid = c(1),
        start_time = start_year,
        skip_time = 2004,
        input_cell_size = 2,
        scale_fac = 1,
        delta_t_total = 10,
        delta_t_step = 2,
        indices = indices,
        variable_names = "est_exp"
      )
    } else {
      out_mature1 <- NULL
    }

    if (3 %in% model_ssid) {
      .predictions <- filter(predictions, ssid == 3)
      if (length(unique(.predictions$year)) == 5) indices <- c(1, 1, 1, 2, 2)
      if (length(unique(.predictions$year)) == 6) indices <- c(1, 1, 1, 2, 2, 2)
      if (length(unique(.predictions$year)) == 7) indices <- c(1, 1, 1, 2, 2, 2, 2)
      if (length(unique(.predictions$year)) == 8) indices <- c(1, 1, 1, 1, 2, 2, 2, 2)

      out_mature3 <- make_trend_data(predictions,
        ssid = c(3),
        start_time = start_year,
        input_cell_size = 2,
        scale_fac = 1,
        delta_t_total = 10,
        delta_t_step = 2,
        indices = indices,
        variable_names = "est_exp"
      )
    } else {
      out_mature3 <- NULL
    }


    if (4 %in% model_ssid) {
      .predictions <- filter(predictions, ssid == 4)
      if (length(unique(.predictions$year)) == 5) indices <- c(1, 1, 1, 2, 2)
      if (length(unique(.predictions$year)) == 6) indices <- c(1, 1, 1, 2, 2, 2)
      if (length(unique(.predictions$year)) == 7) indices <- c(1, 1, 1, 2, 2, 2, 2)
      if (length(unique(.predictions$year)) == 8) indices <- c(1, 1, 1, 1, 2, 2, 2, 2)

      out_mature4 <- make_trend_data(predictions,
        ssid = c(4),
        start_time = start_year,
        input_cell_size = 2,
        scale_fac = 1,
        delta_t_total = 10,
        delta_t_step = 2,
        indices = indices,
        variable_names = "est_exp"
      )
    } else {
      out_mature4 <- NULL
    }

    if (16 %in% model_ssid) {
      .predictions <- filter(predictions, ssid == 16)
      # unique(.predictions$year)
      if (length(unique(.predictions$year)) == 5) indices <- c(1, 1, 1, 2, 2)
      if (length(unique(.predictions$year)) == 6) indices <- c(1, 1, 1, 2, 2, 2)
      if (length(unique(.predictions$year)) == 7) indices <- c(1, 1, 1, 2, 2, 2, 2)
      if (length(unique(.predictions$year)) == 8) indices <- c(1, 1, 1, 1, 2, 2, 2, 2)

      out_mature16 <- make_trend_data(predictions,
        ssid = c(16),
        start_time = start_year,
        skip_time = 2007,
        input_cell_size = 2,
        scale_fac = 1,
        delta_t_total = 10,
        delta_t_step = 2,
        indices = indices,
        variable_names = "est_exp"
      )
    } else {
      out_mature16 <- NULL
    }

    out_mature <- rbind(out_mature1, out_mature3, out_mature4, out_mature16)

    if (exists("adult_predictions")) {
      saveRDS(out_mature, file = paste0(
        "data/", spp, "/trends-",
        spp, covs, "-", ssid_string, "-mature-bioclimatic-priors-", priors, ".rds"
      ))
    } else {
      out_total <- out_mature
      rm(out_mature)
      saveRDS(out_total, file = paste0(
        "data/", spp, "/trends-",
        spp, covs, "-", ssid_string, "-total-bioclimatic-priors-", priors, ".rds"
      ))
    }
  }
}

Recent biomass with trawl footprint

if (!exists("out_mature")) {
  out_recent <- out_total
} else {
  out_recent <- out_mature
}

trawl_footprint <- sf::st_read(dsn = "data/trawl-footprint", layer = "Trawl_footprint")
proj.to <- "+proj=utm +zone=9 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
trawl_footprint <- sf::st_transform(trawl_footprint, crs = proj.to)

current_biomass <- plot_vocc(out_recent,
  fill_col = "est_exp_e",
  fill_label = "Mean\nbiomass (kg/ha)\npost-2010",
  raster_alpha = 1,
  viridis_option = "C",
  vec_aes = NULL,
  transform_col = fourth_root_power
)
current_biomass <- current_biomass +
  geom_sf(data = trawl_footprint, aes(geometry = geometry / 1000), colour = "gray", fill = "gray", alpha = 0.2, lwd = 0.25, inherit.aes = FALSE) +
  coord_sf(xlim = c(180, 820), ylim = c(5350, 6050), datum = NULL) +
  xlab("") +
  ylab("") +
  theme(panel.grid.major = element_line(colour = "transparent")) +
  ggtitle(paste("Predicted", species, "density with trawl footprint overlain"))
current_biomass

Trends in biomass density

if (exists("out_mature")) {
  data <- out_mature
  biotrend_ad <- plot_vocc(data,
    fill_col = "units_per_decade",
    fill_label = "Mature \nbiomass \ntrend \n(kg/ha/decade)",
    raster_alpha = 1,
    vec_aes = NULL,
    white_zero = TRUE,
    high_fill = "deepskyblue3", # "lightseagreen", # limegreen",
    low_fill = "red", # "deeppink3",# "steel blue 4",
    transform_col = fourth_root_power,
    # raster_limits = c(0, quantile(data$units_per_decade, 0.99999)),
    legend_position = c(0.2, 0.25)
  ) + ggtitle(paste0(
    "", species, " ", min(unique(as.numeric(adult_predictions[["year"]]))), "-", max(unique(as.numeric(adult_predictions[["year"]])))
  ))
  biotrend_ad
} else {
  biotrend_ad <- plot_vocc(out_total,
    fill_col = "units_per_decade",
    fill_label = "Biomass \ntrend \n(kg/ha/decade)",
    raster_alpha = 1,
    vec_aes = NULL,
    white_zero = TRUE,
    high_fill = "deepskyblue3", # "greenyellow", # limegreen",
    low_fill = "red", # "steel blue 4",
    transform_col = fourth_root_power,
    legend_position = c(0.2, 0.25)
  ) + ggtitle(paste0(
    "", species, " ", min(unique(as.numeric(predictions[["year"]]))), "-", max(unique(as.numeric(predictions[["year"]])))
  ))
  biotrend_ad
}

Immature biomass density

if available

if (exists("imm_predictions")) {
  p_imm_all <- plot_facet_map(imm_predictions, "est_exp",
    transform_col = fourth_root_power,
    raster_limits = c(0, max_raster)
  ) +
    labs(fill = "kg/ha") +
    ggtitle(paste0("", species, " immature biomass \n "))

  p_imm_all
} else {
  p_imm_all <- grid::grid.rect(gp = grid::gpar(col = "white"))
}

if (exists("adult_predictions")) {
  if (1 %in% model_ssid) {
    imm_predictions1 <- imm_predictions %>% filter(ssid == 1)
  } else {
    imm_predictions1 <- NULL
  }

  if (3 %in% model_ssid) {
    imm_predictions3 <- imm_predictions %>%
      filter(ssid != 16) %>%
      filter(ssid != 4)
  } else {
    imm_predictions3 <- NULL
  }

  if (4 %in% model_ssid) {
    imm_predictions4 <- imm_predictions %>% filter(ssid == 4)
  } else {
    imm_predictions4 <- NULL
  }

  if (16 %in% model_ssid) {
    imm_predictions16 <- imm_predictions %>% filter(ssid == 16)
  } else {
    imm_predictions16 <- NULL
  }

  im_prediction_list <- list(imm_predictions1, imm_predictions3, imm_predictions4, imm_predictions16)
  im_prediction_list <- im_prediction_list[!sapply(im_prediction_list, is.null)]

  p_imm <- list()
  for (i in seq_len(length(model_ssid))) {
    p_imm[[i]] <- plot_facet_map(im_prediction_list[[i]], "est_exp",
      raster_limits = c(0, max_raster),
      # legend_position = "none",
      transform_col = fourth_root_power
    ) +
      labs(fill = "kg/ha") +
      ggtitle(paste0("", species, " immature biomass \n(max = ", max_adult, " kg/ha)"))

    p_imm[[i]]
  }
  print(p_imm)
}

# try(out_imm <- readRDS(paste0("data/", spp, "/trends-", spp, covs, "-", ssid_string, "-immature-bioclimatic-priors-", priors, ".rds")))
rm(out_imm)
if (!exists("out_imm")) {
  if (exists("imm_predictions")) {
    predictions <- imm_predictions

    if (1 %in% model_ssid) {
      .predictions <- filter(predictions, ssid == 1)
      if (length(unique(.predictions$year)) == 5) indices <- c(1, 1, 1, 2, 2)
      if (length(unique(.predictions$year)) == 6) indices <- c(1, 1, 1, 2, 2, 2)
      if (length(unique(.predictions$year)) == 7) indices <- c(1, 1, 1, 2, 2, 2, 2)
      if (length(unique(.predictions$year)) == 8) indices <- c(1, 1, 1, 1, 2, 2, 2, 2)
      out_im1 <- make_trend_data(predictions,
        ssid = c(1),
        start_time = start_year,
        skip_time = 2004,
        input_cell_size = 2,
        scale_fac = 1,
        delta_t_total = 10,
        delta_t_step = 2,
        indices = indices,
        variable_names = "est_exp"
      )
    } else {
      out_im1 <- NULL
    }

    if (3 %in% model_ssid) {
      .predictions <- filter(predictions, ssid == 3)
      if (length(unique(.predictions$year)) == 5) indices <- c(1, 1, 1, 2, 2)
      if (length(unique(.predictions$year)) == 6) indices <- c(1, 1, 1, 2, 2, 2)
      if (length(unique(.predictions$year)) == 7) indices <- c(1, 1, 1, 2, 2, 2, 2)
      if (length(unique(.predictions$year)) == 8) indices <- c(1, 1, 1, 1, 2, 2, 2, 2)

      out_im3 <- make_trend_data(predictions,
        ssid = c(3),
        start_time = start_year,
        input_cell_size = 2,
        scale_fac = 1,
        delta_t_total = 10,
        delta_t_step = 2,
        indices = indices,
        variable_names = "est_exp"
      )
    } else {
      out_im3 <- NULL
    }

    if (4 %in% model_ssid) {
      .predictions <- filter(predictions, ssid == 4)
      unique(.predictions$year)
      if (length(unique(.predictions$year)) == 5) indices <- c(1, 1, 1, 2, 2)
      if (length(unique(.predictions$year)) == 6) indices <- c(1, 1, 1, 2, 2, 2)
      if (length(unique(.predictions$year)) == 7) indices <- c(1, 1, 1, 2, 2, 2, 2)
      if (length(unique(.predictions$year)) == 8) indices <- c(1, 1, 1, 1, 2, 2, 2, 2)

      out_im4 <- make_trend_data(predictions,
        ssid = c(4),
        start_time = start_year,
        input_cell_size = 2,
        scale_fac = 1,
        delta_t_total = 10,
        delta_t_step = 2,
        indices = indices,
        variable_names = "est_exp"
      )
    } else {
      out_im4 <- NULL
    }

    if (16 %in% model_ssid) {
      .predictions <- filter(predictions, ssid == 16)
      if (length(unique(.predictions$year)) == 5) indices <- c(1, 1, 1, 2, 2)
      if (length(unique(.predictions$year)) == 6) indices <- c(1, 1, 1, 2, 2, 2)
      if (length(unique(.predictions$year)) == 7) indices <- c(1, 1, 1, 2, 2, 2, 2)
      if (length(unique(.predictions$year)) == 8) indices <- c(1, 1, 1, 1, 2, 2, 2, 2)

      out_im16 <- make_trend_data(predictions,
        ssid = c(16),
        start_time = start_year,
        input_cell_size = 2,
        scale_fac = 1,
        delta_t_total = 10,
        delta_t_step = 2,
        indices = indices,
        variable_names = "est_exp"
      )
    } else {
      out_im16 <- NULL
    }

    out_imm <- rbind(out_im1, out_im3, out_im4, out_im16)
    saveRDS(out_imm, file = paste0(
      "data/", spp, "/trends-", spp, covs, "-", ssid_string, "-immature-bioclimatic-priors-", priors, ".rds"
    ))
  } else {
    rm(out_imm)
  }
}

if (exists("out_imm")) {
  data <- out_imm
  biotrend_im <- plot_vocc(data,
    fill_col = "units_per_decade",
    fill_label = "Immature \nbiomass \ntrend \n(kg/ha/decade)",
    raster_alpha = 1,
    vec_aes = NULL,
    white_zero = TRUE,
    high_fill = "deepskyblue3", # "greenyellow", # limegreen",
    low_fill = "red", # "steel blue 4",
    transform_col = fourth_root_power,
    legend_position = c(0.2, 0.25)
  ) + ggtitle(paste0(
    "", species, " ", min(unique(imm_predictions[["year"]])), "-", max(unique(imm_predictions[["year"]]))
  ))
  biotrend_im
} else {
  biotrend_im <- grid::grid.rect(gp = grid::gpar(col = "white"))
}

if (exists("imm_predictions")) {
  p_prop_imm <- plot_facet_map(imm_predictions, "prop_imm") +
    labs(fill = "% \nby mass") +
    ggtitle(paste0("", species, " proportion immature"))
  print(p_prop_imm)
} else {
  p_prop_imm <- grid::grid.rect(gp = grid::gpar(col = "white"))
}

Save trend plots

png(
  file = paste0("figs/", spp, "/trends-", spp, covs, "-", ssid_string, "-fixed-everything-noAR1.png"),
  res = 600,
  units = "in",
  width = 11,
  height = 8.5
)
gridExtra::grid.arrange(
  # p_adult, p_imm,
  biotrend_ad, biotrend_im,
  nrow = 1
)
dev.off()

Save spatiotemporal biomass plots

if (exists("adult_predictions")) {
  biomass_plots <- c(p_adult, p_imm)
} else {
  biomass_plots <- c(p_adult)
}

png(
  file = paste0("figs/", spp, "/bioclimatic-biomass-", spp, covs, "-", ssid_string, "-fixed-everything-noAR1.png"),
  res = 600,
  units = "in",
  width = 16,
  height = 8.5
)
gridExtra::grid.arrange(
  grobs = biomass_plots,
  nrow = 2
)
dev.off()

Save plot medley

png(
  file = paste0("figs/", spp, "/bioclimatic-exploratory-", spp, covs, "-", ssid_string, "-fixed-everything-noAR1.png"),
  res = 600,
  units = "in",
  width = 9.5,
  height = 12.5
)

gridExtra::grid.arrange(
  current_biomass,
  p_prop_imm,
  p_adult_all, p_imm_all,
  biotrend_ad, biotrend_im,
  nrow = 3
)
dev.off()

# plot_map_raster <- function(dat, column = "est") {
#     ggplot(dat, aes_string("X", "Y", fill = column)) +
#     geom_raster() +
#     coord_fixed()+
# gfplot::theme_pbs() + theme(axis.title.x = element_blank(), axis.title.y = element_blank())
# }

# plot_map_raster(adult_predictions, "epsilon_st") +
#  ggtitle("Spatiotemporal random effects only") +
#  facet_wrap(~year) +
#  scale_fill_gradient2( high = "Steel Blue 4",
#  low = "Red 3") + coord_fixed()

if (exists("adult_predictions")) {
  # p_st <- list()
  # for (i in seq_len(length(model_ssid))) {
  #   p_st[[i]] <- plot_map_raster(ad_prediction_list[[i]], "epsilon_st") +
  #     ggtitle("Spatiotemporal random effects only") +
  #     facet_wrap(~year) +
  #     scale_fill_gradient2(high = "#66C2A5", #mid = "black",
  #       low = "#5E4FA2"#, trans = fourth_root_power
  #     )
  #   p_st[[i]]
  #  }
  # print(p_st)

  p_adult_st <- list()
  for (i in seq_len(length(model_ssid))) {
    p_adult_st[[i]] <- plot_facet_map(ad_prediction_list[[i]], "do_est",
      # white_zero = TRUE,
      # high_fill = "#66C2A5",
      # low_fill = "#5E4FA2"
      # # raster_limits = c(0, max_raster),
      # legend_position = "none",
      # transform_col = fourth_root_power
    ) +
      # labs(fill = "kg/ha") +
      ggtitle(paste0("", species, " mature spatiotemporal variation"))
    p_adult_st[[i]]
  }
  print(p_adult_st)
  # } else {
  # p_adult_st <- plot_facet_map(predictions, "epsilon_st",
  #   transform_col = fourth_root_power
  # ) +
  #   #labs(fill = "kg/ha") +
  #   ggtitle(paste0("", species, "total spatiotemporal variation"))
  #
  # print(p_adult_st)
}

adult_predictions$year.f <- as.factor(adult_predictions$year)

glimpse(adult_predictions)


data <- mutate(adult_predictions, cell_id = paste0(X, Y)) %>%
  group_by(cell_id) %>%
  mutate(cell_max = max(bioclimatic)) %>%
  group_by(year.f) %>%
  mutate(ann_biomass = sum(est_exp))
unique(data$ann_biomass)
glimpse(data)

# start_biomass <- filter(adult_predictions, year == min(as.numeric(adult_predictions$year)))
#     bio5perc <- sum(adult_predictions$start_density)*0.05
#     s <- sort(annual_biomass$start_density)
#     bio_sum <- cumsum(s)
#     lower_density_threshold <- s[which(bio_sum>bio5perc)[1]]

threshold <- quantile(data$cell_max, 0.2)
library(forcats)
violin_data <- data %>%
  filter(cell_max > threshold) %>%
  mutate(year = fct_reorder(year.f, desc(ann_biomass)))
glimpse(violin_data)


ggplot(violin_data) +
  geom_hline(aes(yintercept = 0)) +
  geom_violin(aes(fct_reorder(year.f, (ann_biomass)), epsilon_st, fill = year.f)) +
  geom_text(aes(fct_reorder(year.f, (ann_biomass)), y = 2.8, label = round(ann_biomass)), check_overlap = TRUE) +
  scale_fill_viridis_d(option = "C") +
  ylab("Spatiotemporal variation") +
  xlab("Year") +
  # coord_flip() +
  theme(legend.position = "none") +
  ggtitle(paste("Local", species, "abundance deviations not explained by climate and overall abundance"))

# data <- adult_predictions %>% group_by (year) %>%
#   mutate(mean_temp = mean(temp)) %>%
#   filter (est_exp > quantile(est_exp, 0.05)) %>%
#   filter (do_est > 1.5)
#
# ggplot(data, aes(temp, epsilon_st)) +
#   geom_point(alpha = 0.025) +
#   geom_smooth(method="lm", formula = y~poly(x,2)) +
#   facet_wrap(~year) +
#   scale_color_viridis_d(option = "C") +
#   guides(color = FALSE) +
#   ylab("Spatiotemporal variation") +
#   ggtitle(paste("Local", species, "abundance deviations not explained by climate and overall abundance"))

Save spatiotemporal residual plots

# if(exists("adult_predictions")) { resid_plots <- c(p_st, p_imm) } else {
#   biomass_plots <- c(p_adult)
# }

png(
  file = paste0("figs/", spp, "/bioclimatic-residuals-", spp, covs, "-", ssid_string, "-fixed-everything-noAR1.png"),
  res = 600,
  units = "in",
  width = 16,
  height = 8.5
)
gridExtra::grid.arrange(
  grobs = p_adult_st,
  nrow = 2
)
dev.off()

if (exists("adult_predictions")) {
  p_bioclim <- list()
  for (i in seq_len(length(model_ssid))) {
    p_bioclim[[i]] <- plot_facet_map(ad_prediction_list[[i]], "bioclimatic",
      raster_limits = c(0, max_raster),
      # legend_position = "none",
      transform_col = fourth_root_power
    ) +
      labs(fill = "kg/ha") +
      ggtitle(paste0("", species, " predicted biomass based on climate \nand fixed spatial variation"))
    p_bioclim[[i]]
  }
}
print(p_bioclim)

Save bioclimatic predictions plots

# if(exists("adult_predictions")) { resid_plots <- c(p_st, p_imm) } else {
#   biomass_plots <- c(p_adult)
# }

png(
  file = paste0("figs/", spp, "/bioclimatic-predictions-", spp, covs, "-", ssid_string, "-priors-", priors, ".png"),
  res = 600,
  units = "in",
  width = 16,
  height = 8.5
)
gridExtra::grid.arrange(
  grobs = p_bioclim,
  nrow = 2
)
dev.off()

if (exists("imm_predictions")) {
  p_prop_imm <- plot_facet_map(imm_predictions, "prop_imm") +
    labs(fill = "% \nby mass") +
    ggtitle(paste0("", species, " proportion immature"))
  print(p_prop_imm)
} else {
  p_prop_imm <- grid::grid.rect(gp = grid::gpar(col = "white"))
}

# Velocity of climate change ---------------------------------------------

d <- ad_prediction_list[[2]]

d$layer <- d$bioclimatic

# Just to confirm the trends we're seeing are real,
# mod = mgcv::gam(est ~ s(year, k=5) + te(X,Y) + log(depth),
#   data=d)
# plot(mod)

# This parameter controls how the original 2-km projection is aggregated.
# for example, a value of 5 means that the raster would be reprojected to 10km grid
scale_fac <- 1

# create a RasterBrick
# raster for each year
rlist <- list()
for (i in 1:length(unique(d$year))) {
  rlist[[i]] <- raster::rasterFromXYZ(dplyr::filter(d, year == unique(d$year)[i]) %>%
    dplyr::select(X, Y, layer))
  rlist[[i]] <- raster::aggregate(rlist[[i]], fact = scale_fac)
}
# stack rasters into layers -> rasterbrick
rstack <- raster::stack(rlist[[1]], rlist[[2]])
for (i in 3:length(rlist)) {
  rstack <- raster::stack(rstack, rlist[[i]])
}
rbrick <- raster::brick(rstack)

# Then calculate the trend per pixel:
slopedat <- vocc::calcslope(rbrick)
# Then get the mean temperature for a time period and calculate the spatial gradient:
allyears <- rep(1, raster::nlayers(rbrick))
mnsst <- raster::stackApply(rbrick, indices = allyears, fun = mean)
spatx <- vocc::spatialgrad(mnsst)
# Now we can calculate the VoCC:
velodf <- vocc::calcvelocity(spatx, slopedat, y_dist = 1)

# Mapping it again is straightforward:
rtrend <- rgrad_lon <- rgrad_lat <- rvocc <- angle <- magn <- raster::raster(rbrick)
rgrad_lat[spatx$icell] <- spatx$NS # latitude shift, NS
rgrad_lon[spatx$icell] <- spatx$WE # longitude shift, WE
rtrend[slopedat$icell] <- -1 * slopedat$slope
rvocc[velodf$icell] <- velodf$velocity

# convert to data frames for ggplot
rtrend_df <- as.data.frame(raster::rasterToPoints(rtrend)) %>%
  dplyr::rename(trend = layer)
rgradlat_df <- as.data.frame(raster::rasterToPoints(rgrad_lat)) %>%
  dplyr::rename(gradNS = layer)
rgradlon_df <- as.data.frame(raster::rasterToPoints(rgrad_lon)) %>%
  dplyr::rename(gradWE = layer)
rvocc_df <- as.data.frame(raster::rasterToPoints(rvocc)) %>%
  dplyr::rename(velocity = layer)

# create ggquiver plots. need dataframe of lon, lat, delta_lon, delta_lat, trend, velocity
df <- dplyr::left_join(rtrend_df, rgradlat_df, by = c("x", "y")) %>%
  dplyr::left_join(rgradlon_df, by = c("x", "y")) %>%
  dplyr::left_join(rvocc_df, by = c("x", "y"))

gtrend <- ggplot(df, aes(x, y, fill = -trend)) +
  geom_raster() +
  scale_fill_gradient2(low = scales::muted("blue"), high = scales::muted("red")) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Trend") + coord_fixed()

# spatial gradient plot
quantile_cutoff <- 0.05 # for plotting
df <- dplyr::mutate(df,
  u_velo = trend / gradWE,
  v_velo = trend / gradNS,
  ulow = quantile(u_velo, quantile_cutoff),
  uhi = quantile(u_velo, 1 - quantile_cutoff),
  vlow = quantile(v_velo, quantile_cutoff),
  vhi = quantile(v_velo, 1 - quantile_cutoff),
  u_velo = ifelse(u_velo < ulow, ulow, u_velo),
  u_velo = ifelse(u_velo > uhi, uhi, u_velo),
  v_velo = ifelse(v_velo < vlow, vlow, v_velo),
  v_velo = ifelse(v_velo > vhi, vhi, v_velo)
) %>%
  dplyr::select(-ulow, -uhi, -vlow, -vhi)

# gradient plot with ggquiver
ggrad <- ggplot(df) +
  ggquiver::geom_quiver(aes(x, y, u = gradNS, v = gradWE)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Gradient") + coord_fixed()

# velocity plot
gvocc <- ggplot(mutate(df, `Trend\n(kg/decade)` = -trend)) +
  ggquiver::geom_quiver(aes(x, y,
    u = u_velo, v = v_velo,
    colour = `Trend\n(kg/decade)`
  ),
  vecsize = 2
  ) +
  scale_colour_gradient2(low = scales::muted("blue"), high = scales::muted("red")) +
  xlab("UTM") + ylab("UTM") + # ggtitle("Velocity") +
  gfplot::theme_pbs()

# coast <- gfplot:::load_coastline(range(surv$longitude) + c(-1, 1),
#   range(surv$latitude) + c(-1, 1),
#   utm_zone = 9
# )

# isobath <- gfplot:::load_isobath(range(surv$longitude) + c(-5, 5),
#   range(surv$latitude) + c(-5, 5),
#   bath = c(100, 200, 300, 500), utm_zone = 9
# )

# library(ggnewscale)

gvocc <- gvocc +
  labs(colour = "Local\nclimate trend\n(°C/decade)") +
  ggnewscale::new_scale_color() +
  # geom_path(
  #   data = isobath, aes_string(
  #     x = "X", y = "Y",
  #     group = "paste(PID, SID)", colour = "PID"
  #   ),
  #   inherit.aes = FALSE, lwd = 0.4, alpha = 0.4
  # ) +
  scale_colour_continuous(low = "grey70", high = "grey10") +
  # labs(colour = "Depth") +
  guides(colour = FALSE) +
  coord_fixed(xlim = range(df$x) + c(-3, 3), ylim = range(df$y) + c(-3, 3))

# gvocc <- gvocc + geom_polygon(
#   data = coast, aes_string(x = "X", y = "Y", group = "PID"),
#   fill = "grey87", col = "grey70", lwd = 0.2
# )


pbs-assess/gfranges documentation built on Dec. 13, 2021, 4:50 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close