In pbs-assess/gfranges: Groundfish Ranges

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

Load predicted values

pred_temp_do <- readRDS(here::here("analysis/tmb-sensor-explore/models/predicted-DO.rds"))
predicted1 <- pred_temp_do %>% filter(ssid == 4)
predicted <- list(predicted1, predicted1)

Define subset of the spatial range

do <- plot_facet_map(predicted[[1]], "do_est", transform_col = no_trans) +
  labs(fill = "ml/L") +
  ggtitle("Bottom DO")
print(do)

Biannual changes starting 2008 and ending 2018

out1 <- make_vector_data(predicted,
  variable_names = c("do_est", "temp"),
  ssid = c(4),
  start_time = 2008,
  end_time = 2010,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),  
  # match_logic = c(">=", "=="),
  # plus_minus = c(1, 1)) # vector of plus/minus threshold(s) to define climate match.
  min_thresholds = c(0.5, 3),
  max_thresholds = c(Inf, 1),
  round_fact = 10
  )
out1$do_est.units_per_decade <- out1$do_est.units_per_decade / 5
out1$temp.units_per_decade <- out1$temp.units_per_decade / 5

gvocc1a <- plot_vocc(out1,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "do_est.units_per_decade",
  fill_label = "ml/L DO",
  high_fill = "Steel Blue 4",
  low_fill = "Red 3",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  raster_limits = c(-2, 2))

gvocc1b <- plot_vocc(out1,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "temp.units_per_decade",
  fill_label = "°C",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  # transform_col = fourth_root_power,
  raster_limits = c(-2, 2))


gvocc1 <- gridExtra::grid.arrange(gvocc1b, gvocc1a, nrow = 1, top = "Change between 2008 and 2010")

out2 <- make_vector_data(predicted,
  variable_names = c("do_est", "temp"),
  ssid = c(4),
  start_time = 2010,
  end_time = 2012,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  min_thresholds = c(0.5, 3),
  max_thresholds = c(Inf, 1),
  round_fact = 10)
out2$do_est.units_per_decade <- out2$do_est.units_per_decade / 5
out2$temp.units_per_decade <- out2$temp.units_per_decade / 5

gvocc2a <- plot_vocc(out2,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "do_est.units_per_decade",
  fill_label = "ml/L DO",
  high_fill = "Steel Blue 4",
  low_fill = "Red 3",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  raster_limits = c(-2, 2))

gvocc2b <- plot_vocc(out2,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "temp.units_per_decade",
  fill_label = "°C",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  # transform_col = fourth_root_power,
  raster_limits = c(-2, 2)
)


gvocc2 <- gridExtra::grid.arrange(gvocc2b, gvocc2a, nrow = 1, top = "Change between 2010 and 2012")

out3 <- make_vector_data(predicted,
  variable_names = c("do_est", "temp"),
  ssid = c(4),
  start_time = 2012,
  end_time = 2014,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  min_thresholds = c(0.5, 3),
  max_thresholds = c(Inf, 1),
  round_fact = 10)
out3$do_est.units_per_decade <- out3$do_est.units_per_decade / 5
out3$temp.units_per_decade <- out3$temp.units_per_decade / 5
View(out3)

gvocc3a <- plot_vocc(out3,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "do_est.units_per_decade",
  fill_label = "ml/L DO",
  high_fill = "Steel Blue 4",
  low_fill = "Red 3",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  raster_limits = c(-2, 2))

gvocc3b <- plot_vocc(out3,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "temp.units_per_decade",
  fill_label = "°C",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  # transform_col = fourth_root_power,
  raster_limits = c(-2, 2))

gvocc3 <- gridExtra::grid.arrange(gvocc3b, gvocc3a, nrow = 1, top = "Change between 2012 and 2014")

out4 <- make_vector_data(predicted,
  variable_names = c("do_est", "temp"),
  ssid = c(4),
  start_time = 2014,
  end_time = 2016,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  min_thresholds = c(0.5, 3),
  max_thresholds = c(Inf, 1),
  round_fact = 10)
out4$do_est.units_per_decade <- out4$do_est.units_per_decade / 5
out4$temp.units_per_decade <- out4$temp.units_per_decade / 5


gvocc4a <- plot_vocc(out4,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "do_est.units_per_decade",
  fill_label = "ml/L DO",
  high_fill = "Steel Blue 4",
  low_fill = "Red 3",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE 
)

gvocc4b <- plot_vocc(out4,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "temp.units_per_decade",
  fill_label = "°C",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  # transform_col = fourth_root_power,
  raster_limits = c(-2, 2))

gvocc4 <- gridExtra::grid.arrange(gvocc4b, gvocc4a, nrow = 1, top = "Change between 2014 and 2016")

out5 <- make_vector_data(predicted,
  variable_names = c("do_est", "temp"),
  ssid = c(4),
  start_time = 2016,
  #skip_time = 2016,
  end_time = 2018,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 2,
  delta_t_step = 2,
  indices = c(1, 2),
  min_thresholds = c(0.5, 3),
  max_thresholds = c(Inf, 1),
  round_fact = 10)
out5$do_est.units_per_decade <- out5$do_est.units_per_decade / 5
out5$temp.units_per_decade <- out5$temp.units_per_decade / 5


gvocc5a <- plot_vocc(out5,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "do_est.units_per_decade",
  fill_label = "ml/L DO",
  high_fill = "Steel Blue 4",
  low_fill = "Red 3",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE 
)

gvocc5b <- plot_vocc(out5,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "temp.units_per_decade",
  fill_label = "°C",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  # transform_col = fourth_root_power,
  raster_limits = c(-2, 2))

gvocc5 <- gridExtra::grid.arrange(gvocc5b, gvocc5a, nrow = 1, top = "Change between 2016 and 2018")

out6 <- make_vector_data(predicted,
  variable_names = c("do_est", "temp"),
  ssid = c(4),
  start_time = 2014,
  skip_time = 2016,
  end_time = 2018,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 4,
  delta_t_step = 2,
  indices = c(1, 2),
  min_thresholds = c(0.5, 3),
  max_thresholds = c(Inf, 1),
  round_fact = 10)
out6$do_est.units_per_decade <- out6$do_est.units_per_decade / 5
out6$temp.units_per_decade <- out6$temp.units_per_decade / 5


gvocc6a <- plot_vocc(out6,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "do_est.units_per_decade",
  fill_label = "ml/L DO",
  high_fill = "Steel Blue 4",
  low_fill = "Red 3",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  raster_limits = c(-2, 2))

gvocc6b <- plot_vocc(out6,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 50,
  NA_label = "*",
  fill_col = "temp.units_per_decade",
  fill_label = "°C",
  raster_alpha = 1,
  vec_alpha = 0.5,
  axis_lables = FALSE,
  raster_limits = c(-2, 2))

gvocc6 <- gridExtra::grid.arrange(gvocc6b, gvocc6a, nrow = 1, top = "Change between 2014 and 2018")

png(
  file = "figs/biannual_temp-do-_change-4-asym3-inf.png", 
  res = 600,
  units = "in",
  width = 8, 
  height = 20
) 
gridExtra::grid.arrange(gvocc1, gvocc2, gvocc3, gvocc4, gvocc5, gvocc6, nrow = 6, top = grid::textGrob("VOCC vectors for <0.5 ml/L biannual decrease in bottom DO and -3 to +1 °C")) 
dev.off()

Cumulative change between last two decades

VOCC for declining levels of DO

unique(predicted1$year)
starttime1 <- Sys.time()

declining_var <- "do_est"
min_threshold <- 0.2

out0.5d <- make_vector_data(predicted1,
  variable_names = c(declining_var),
  #ssid = c(4),
  start_time = 2008,
  skip_time = 2016,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 6,
  delta_t_step = 2,
  indices = c(1, 1, 1, 2, 2),
  min_thresholds = c(min_threshold),
  max_thresholds = c(Inf),
  round_fact = 10
)
endtime1 <- Sys.time()
time1 <- round(starttime1 - endtime1)

saveRDS(out0.5d, file = paste0("data/", 
      "vocc-", ssid_string, declining_var, min_threshold, "decline.rds"))

do4 <- plot_vocc(out0.5d,
  fill_col = "var_1_s",
  fill_label = "ml/L",
  raster_alpha = 1,
  vec_aes = NULL,
  #viridis_option = "A",
  raster_limits = c(0, 5 )
  )

do4 <- do4 + ggtitle("Mean DO 2008-2012")
do4

voccd4 <- plot_vocc(out0.5d,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 100,
  fill_col = "units_per_decade",
  fill_label = "ml/L DO\nper decade",
  raster_alpha = 1,
  vec_alpha = 0.75,
  axis_lables = FALSE,
  viridis_option = "C",
  viridis_dir = -1,
  NA_label = "*",
  high_fill = "Steel Blue 4",
  low_fill = "Red 3",
  min_vec_plotted = 2,
  raster_limits = c(-1.75, 1.75))
voccd4 <- voccd4 + ggtitle("VOCC <1 ml/L decline (excludes 2016)")
voccd4

increasing_var <- "temp" 
max_threshold <- 0.5

out0.5t <- make_vector_data(predicted1,
  variable_names = c(increasing_var),
  #ssid = c(4),
  start_time = 2008,
  skip_time = 2016,
  input_cell_size = 2,
  scale_fac = 1,
  delta_t_total = 6,
  delta_t_step = 2,
  indices = c(1, 1, 1, 2, 2),
  min_thresholds = c(Inf),
  max_thresholds = c(max_threshold),
  round_fact = 10
)


saveRDS(out0.5t, file = paste0("data/", 
      "vocc-", ssid_string, increasing_var, max_threshold, "increase.rds"))

temp4 <- plot_vocc(out0.5t,
  fill_col = "var_1_s",
  fill_label = "°C ",
  raster_alpha = 1,
  vec_aes = NULL,
  viridis_option = "C",
  raster_limits = c(4.5, 11)
  )

temp4 <- temp4 + ggtitle("Mean temperature 2008-2012")
temp4

vocct4 <- plot_vocc(out0.5t,
  vec_aes = "distance",
  vec_lwd_range = c(0.2, 0.5),
  max_vec_plotted = 100,
  fill_col = "units_per_decade",
  fill_label = "°C per \ndecade",
  vec_alpha = 0.75,
  raster_alpha = 1,
  min_vec_plotted = 2,
  white_zero = TRUE,
  raster_limits = c(-0.5, 2))
vocct4 <- vocct4 + ggtitle("VOCC <0.5°C increase (excludes 2016)")
vocct4

png(
  file = "figs/vocc-temp-do-2008-2018-4-asym-1.png", 
  res = 600,
  units = "in",
  width = 7, 
  height = 7
) 

gridExtra::grid.arrange(temp4, do4, vocct4, voccd4, nrow = 2)
dev.off()