In pbs-assess/gfranges: Groundfish Ranges

# install.packages("yaImpute") # install package for k-nearest neighbour (kNN) search
# install.packages("raster") # required for calcslope + all raster functions
# install.package(gfplot) # required for some plotting and data functions

library(dplyr)
library(ggplot2)
library(sdmTMB)
library(gfranges)

Retrieve raw data

all_depth <- readRDS("../tmb-sensor-explore/data/dat-sensor-trawl-all-depth.rds")

# add and scale predictors after filtering for survey(s) of interest
data <- all_depth$data %>%
  filter(depth_max > 10) %>% # removes data if sensor stayed too near surface
  # filter(ssid != 1) %>%
  # filter(ssid != 3) %>%
  # filter(ssid != 4) %>%
  # filter(ssid != 16) %>%
  gfplot:::scale_survey_predictors()

data <- data[data$fishing_event_id != 481861, ]
# add any new predictors
# data <- data %>%
#   dplyr::mutate(
#   DOY = lubridate::yday(date),
#   shallow = ifelse(depth > 35, 0, 1)
# )

Or use old data

data <- readRDS("../tmb-sensor-explore/data/sensor-data-processed.rds")
data <- data %>%
  mutate(depth = depth_m) %>%
  filter(depth_m > 10) %>% # removes data if sensor stayed too near surface
  gfplot:::scale_survey_predictors()

data <- data[data$fishing_event_id != 481861, ]

Create prediction grids for each survey area so that only years with data are included

# choose base year(s) to create grid from
dummy_year <- c(2004, 2005)

# create grids for each ssid separately
ssid <- 1
survey_abbrev <- "SYN QCS"
nd_1 <- spatiotemporal_grid(data, ssid, survey_abbrev, dummy_year)
unique(nd_1$year)

ssid <- 3
survey_abbrev <- "SYN HS"
nd_3 <- spatiotemporal_grid(data, ssid, survey_abbrev, dummy_year)
unique(nd_3$year)

ssid <- 4
survey_abbrev <- "SYN WCVI"
nd_4 <- spatiotemporal_grid(data, ssid, survey_abbrev, dummy_year)
unique(nd_4$year)

ssid <- 16
survey_abbrev <- "SYN WCHG"
nd_16 <- spatiotemporal_grid(data, ssid, survey_abbrev, dummy_year = 2006)
unique(nd_16$year)

nd_all <- rbind(nd_1, nd_3, nd_4, nd_16)
nd <- nd_all %>% dplyr::mutate(shallow = ifelse(depth > 35, 0, 1))

Run sdmTMB model (if not done previously).

# # choose the spatial range to model
dat <- data #%>% dplyr::filter(year > 2003) # %>% select(year, ssid, X, Y, temperature_c, depth_scaled, depth_scaled2)
#nd <- nd_1 # %>% filter(ssid == 4)

spde <- sdmTMB::make_spde(dat$X, dat$Y, n_knots = 500)
sdmTMB::plot_spde(spde)

m_temp <- sdmTMB::sdmTMB(dat,
  temperature_c ~ 0 + as.factor(year),
  time_varying = ~ 0 + depth_scaled + depth_scaled2,
  time = "year", spde = spde,
  family = gaussian(link = "identity"),
  ar1_fields = TRUE, # maybe TRUE is better for all areas combined?
  include_spatial = TRUE,
  silent = FALSE
)

# stopifnot(m_temp$model$convergence == 0L)
# m_temp
# # Warning messages:
# #   1: In doTryCatch(return(expr), name, parentenv, handler) :
# #   restarting interrupted promise evaluation
#
# saveRDS(m_temp, file = "data/m_temp_post2003.rds")
saveRDS(m_temp, file = "data/m_temp_allyears.rds")

Calculate predicted values

m_temp <- readRDS("../VOCC/data/m_temp_allyears.rds")

# # must be filtered to match years in model
# m_temp <- readRDS("../VOCC/data/m_temp_allpost2003.rds")
# nd <- nd %>% dplyr::filter(year > 2003)

temp <- predict(m_temp, newdata = nd)

Define subset of the spatial range

temperature <- temp %>% 
  # filter(ssid!=3) %>% 
  filter(ssid!=4) %>% 
  filter(ssid!=16) %>% 
  filter(year>2004)

#glimpse(temperature)

temp <-plot_facet_map(temperature, "est", transform_col = no_trans) + 
  labs(fill="kg/ha") +
  ggtitle("Bottom temperature")
print(temp)

Biannual changes starting 2005 and ending 2017

out1 <- make_vector_data(temperature,
  ssid = c(1,3),
  start_time = 2005,
  end_time = 2007,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  thresholds = c(0.5) # vector of plus/minus threshold(s) to define climate match.
)
gvocc1 <- plot_vocc(out1,
  low_col = "white",
  mid_col = "white",
  high_col = "grey87",
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = ".",
  fill_col = "var_1_e",
  fill_label = "2007",
  raster_alpha = 1,
  vec_alpha = 0.25,
  axis_lables = FALSE,
  viridis_option = "D",
  transform_col = no_trans,
  raster_limits = c(5, 13)
)
#gvocc1 <- gvocc1 + ggtitle("VOCC vectors for <1°C change in bottom temperature")
gvocc1

out2 <- make_vector_data(temperature,
  ssid = c(1,3),
  start_time = 2007,
  end_time = 2009,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  thresholds = c(0.5) # vector of plus/minus threshold(s) to define climate match.
)
gvocc2 <- plot_vocc(out2,
  low_col = "white",
  mid_col = "white",
  high_col = "grey87",
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = ".",
  fill_col = "var_1_e",
  fill_label = "2009",
  raster_alpha = 1,
  vec_alpha = 0.25,
  axis_lables = FALSE,
  viridis_option = "D",
  transform_col = no_trans,
  raster_limits = c(5, 13)
)

#gvocc2 <- gvocc2 + ggtitle(" ")
gvocc2

out3 <- make_vector_data(temperature,
  ssid = c(1,3),
  start_time = 2009,
  end_time = 2011,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  thresholds = c(0.5) # vector of plus/minus threshold(s) to define climate match.
)
gvocc3 <- plot_vocc(out3,
  low_col = "white",
  mid_col = "white",
  high_col = "grey87",
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = ".",
  fill_col = "var_1_e",
  fill_label = "2011",
  raster_alpha = 1,
  vec_alpha = 0.25,
  axis_lables = FALSE,
  viridis_option = "D",
  transform_col = no_trans,
  raster_limits = c(5, 13)
)
#gvocc3 <- gvocc3 + ggtitle(" ")
gvocc3

out4 <- make_vector_data(temperature,
  ssid = c(1,3),
  start_time = 2011,
  end_time = 2013,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  thresholds = c(0.5) # vector of plus/minus threshold(s) to define climate match.
)
gvocc4 <- plot_vocc(out4,
  low_col = "white",
  mid_col = "white",
  high_col = "grey87",
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = ".",
  fill_col = "var_1_e",
  fill_label = "2013",
  raster_alpha = 1,
  vec_alpha = 0.25,
  axis_lables = FALSE,
  viridis_option = "D",
  transform_col = no_trans,
  raster_limits = c(5, 13)
)
#gvocc4 <- gvocc4 + ggtitle(" ")
gvocc4

out5 <- make_vector_data(temperature,
  ssid = c(1,3),
  start_time = 2013,
  end_time = 2015,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  thresholds = c(0.5) # vector of plus/minus threshold(s) to define climate match.
)
gvocc5 <- plot_vocc(out5,
  low_col = "white",
  mid_col = "white",
  high_col = "grey87",
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = ".",
  fill_col = "var_1_e",
  fill_label = "2015",
  raster_alpha = 1,
  vec_alpha = 0.25,
  axis_lables = FALSE,
  viridis_option = "D",
  transform_col = no_trans,
  raster_limits = c(5, 13)
)
#gvocc5 <- gvocc5 + ggtitle(" ")
gvocc5

out6 <- make_vector_data(temperature,
  ssid = c(1,3),
  start_time = 2015,
  end_time = 2017,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  thresholds = c(0.5) # vector of plus/minus threshold(s) to define climate match.
)
gvocc6 <- plot_vocc(out6,
  low_col = "white",
  mid_col = "white",
  high_col = "grey87",
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = ".",
  fill_col = "var_1_e",
  fill_label = "2017",
  raster_alpha = 1,
  vec_alpha = 0.25,
  axis_lables = FALSE,
  viridis_option = "D",
  transform_col = no_trans,
  raster_limits = c(5, 13)
)
#gvocc6 <- gvocc6 + ggtitle(" ")
gvocc6

png(file = "biannual_temp_change.png",   # The directory you want to save the file in
res = 600,
units = 'in',
width = 10, # The width of the plot in inches
height = 8) # The height of the plot in inches

gridExtra::grid.arrange(gvocc1, gvocc2, gvocc3, gvocc4, gvocc5, gvocc6, nrow = 2, top = grid::textGrob("VOCC vectors for <0.5°C biannual changes in bottom temperature"))#,gp=gpar(fontsize=20,font=3)) 
dev.off()

Cumulative change between last two decades VOCC for 0.5 degree change and 1 degree change

out0.5d <- make_vector_data(temperature,
  ssid = c(1,3),
  start_time = 2005,
  #skip_time = 2011,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 10,
  delta_t_step = 2,
  indices = c(1, 1, 1, 2, 2, 2, 2),
  thresholds = c(0.5) # vector of plus/minus threshold(s) to define climate match.
)

gvocc <- plot_vocc(out0.5,
  low_col = "white",
  mid_col = "white",
  high_col = "grey87",
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 100,
  fill_col = "var_1_e",
  fill_label = "2013-2017",
  raster_alpha = 1,
  vec_alpha = 0.25,
  axis_lables = FALSE,
  viridis_option = "D",
  transform_col = no_trans
)
gvocc <- gvocc + ggtitle("VOCC vectors for <0.5°C change")
gvocc

out1d <- make_vector_data(temperature,
  ssid = c(1,3),
  start_time = 2005,
  #skip_time = 2011,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 10,
  delta_t_step = 2,
  indices = c(1, 1, 1, 2, 2, 2, 2),
  thresholds = c(1) # vector of plus/minus threshold(s) to define climate match.
)

gvocc1d <- plot_vocc(out1d,
  low_col = "white",
  mid_col = "white",
  high_col = "grey87",
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 100,
  fill_col = "var_1_e",
  fill_label = "2013-2017",
  raster_alpha = 1,
  vec_alpha = 0.25,
  axis_lables = FALSE,
  viridis_option = "D",
  transform_col = no_trans
)
gvocc1d <- gvocc1d + ggtitle("VOCC vectors for <1°C change")
gvocc1d

gtrend <- plot_vocc(out1d,
  fill_col = "units_per_decade",
  fill_label = "°C per decade",
  raster_alpha = 1,
  vec_aes = NULL,
  transform_col = no_trans
)
gtrend <- gtrend + ggtitle("Bottom temperature trend 2005-2017")
gtrend

png(file = "temp-2005-2017-1n3.png",   # The directory you want to save the file in
res = 600,
units = 'in',
width = 10, # The width of the plot in inches
height = 7) # The height of the plot in inches

gridExtra::grid.arrange(gtrend, gvocc, gvocc1d, nrow = 1)
dev.off()