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("region =", region)
paste("ssids modelled separately =", bind_predictions)
paste("model label =", covs)
paste("priors =", priors)
spp <- gsub(" ", "-", gsub("\\/", "-", tolower(species)))
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")
}
if (region == "both odd year surveys") {
survey <- c("SYN QCS", "SYN HS")
model_ssid <- c(1, 3)
ssid_string <- paste0(model_ssid, collapse = "n")
}
if (region == "West Coast Vancouver Island") {
survey <- c("SYN WCVI")
model_ssid <- c(4)
ssid_string <- paste0(model_ssid, collapse = "n")
}
if (region == "West Coast Haida Gwaii") {
survey <- c("SYN WCHG")
model_ssid <- c(16)
ssid_string <- paste0(model_ssid, collapse = "n")
}
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,
"/predictions-", spp, covs, "-1n3n4n16-mature-biomass-prior-", priors, ".rds"
)) %>% filter(ssid %in% model_ssid))
try(imm_predictions <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-1n3n4n16-imm-biomass-prior-", priors, ".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-biomass-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,
"/predictions-", spp, covs, "-", ssid_string, "-mature-biomass-prior-", priors, ".rds"
))
})
try({imm_predictions <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-", ssid_string, "-imm-biomass-prior-", priors, ".rds"
))
})
# OR bind together predictions from different models
if (bind_predictions) {
try({pred1n3 <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-1n3-mature-biomass-prior-", priors, ".rds"
))
})
try({pred4 <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-4-mature-biomass-prior-", priors, ".rds"
))
})
try({pred16 <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-16-mature-biomass-prior-", priors, ".rds"
))
})
if (!exists("pred1n3")) pred1n3 <- NULL
if (!exists("pred16")) pred16 <- NULL
if (!exists("pred4")) pred4 <- NULL
adult_predictions <- rbind(pred1n3, pred4, pred16)
try({pred1n3_imm <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-1n3-imm-biomass-prior-", priors, ".rds"
))
})
try({pred4_imm <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-4-imm-biomass-prior-", priors, ".rds"))
})
try({pred16_imm <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-16-imm-biomass-prior-", priors, ".rds"))
})
if (!exists("pred1n3_imm")) pred1n3_imm <- NULL
if (!exists("pred16_imm")) pred16_imm <- NULL
if (!exists("pred4_imm")) pred4_imm <- NULL
imm_predictions <- rbind(pred1n3_imm, pred4_imm, pred16_imm)
}
}
if (!exists("adult_predictions")) {
paste("Maturity split biomass predictions are not available.")
try({predictions <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-", ssid_string, "-total-biomass-prior-", priors, ".rds"))
})
if (!exists("predictions")) {
try({pred1n3_all <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-1n3-total-biomass-prior-", priors, ".rds"))
})
try({pred4_all <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-4-total-biomass-prior-", priors, ".rds"))
})
try({pred16_all <- readRDS(paste0("data/", spp,
"/predictions-", spp, covs, "-16-total-biomass-prior-", priors, ".rds"))
})
try({
predictions <- rbind(pred1n3_all, pred4_all, pred16_all)
})
if (!exists("predictions")) paste("Total biomass predictions are not available.")
}
}
Caluclate variables of interest from predictions
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
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-biomass-prior-", priors, ".rds"
))
saveRDS(imm_predictions, file = paste0(
"data/", spp,
"/all-predictions-", spp, covs, "-imm-biomass-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$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-biomass-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-biomass-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-biomass-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-biomass-priors-", priors, ".rds"))
} else {
out_total <- out_mature
rm(out_mature)
saveRDS(out_total, file = paste0("data/", spp, "/trends-",
spp, covs, "-", ssid_string, "-total-biomass-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,
legend_position = c(0.2, 0.25)
) + ggtitle(paste0(
"", species, " ", min(unique(adult_predictions[["year"]])), "-", max(unique(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(predictions[["year"]])), "-", max(unique(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-biomass-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-biomass-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, "-priors-", priors, ".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, "/biomass-", spp, covs, "-", ssid_string, "-priors-", priors, ".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, "/exploratory-", spp, covs, "-", ssid_string, "-priors-", priors, ".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 <- plot_facet_map(adult_predictions, "epsilon_st",
# white_zero = TRUE
# #raster_limits = c(0, max_raster),
# #legend_position = "none",
# #transform_col = fourth_root_power
# ) +
# #labs(fill = "kg/ha") +
# ggtitle(paste0("", species, "mature spatiotemporal variation"))
# 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)
}
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, "/biotic-residuals-", spp, covs, "-", ssid_string, "-priors-", priors, ".png"),
res = 600,
units = "in",
width = 16,
height = 8.5
)
gridExtra::grid.arrange(
grobs = p_st,
nrow = 2
)
dev.off()
pbs-assess/gfranges documentation built on Dec. 13, 2021, 4:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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("region =", region) paste("ssids modelled separately =", bind_predictions) paste("model label =", covs) paste("priors =", priors) spp <- gsub(" ", "-", gsub("\\/", "-", tolower(species))) 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") } if (region == "both odd year surveys") { survey <- c("SYN QCS", "SYN HS") model_ssid <- c(1, 3) ssid_string <- paste0(model_ssid, collapse = "n") } if (region == "West Coast Vancouver Island") { survey <- c("SYN WCVI") model_ssid <- c(4) ssid_string <- paste0(model_ssid, collapse = "n") } if (region == "West Coast Haida Gwaii") { survey <- c("SYN WCHG") model_ssid <- c(16) ssid_string <- paste0(model_ssid, collapse = "n") }
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, "/predictions-", spp, covs, "-1n3n4n16-mature-biomass-prior-", priors, ".rds" )) %>% filter(ssid %in% model_ssid)) try(imm_predictions <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-1n3n4n16-imm-biomass-prior-", priors, ".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-biomass-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, "/predictions-", spp, covs, "-", ssid_string, "-mature-biomass-prior-", priors, ".rds" )) }) try({imm_predictions <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-", ssid_string, "-imm-biomass-prior-", priors, ".rds" )) }) # OR bind together predictions from different models if (bind_predictions) { try({pred1n3 <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-1n3-mature-biomass-prior-", priors, ".rds" )) }) try({pred4 <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-4-mature-biomass-prior-", priors, ".rds" )) }) try({pred16 <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-16-mature-biomass-prior-", priors, ".rds" )) }) if (!exists("pred1n3")) pred1n3 <- NULL if (!exists("pred16")) pred16 <- NULL if (!exists("pred4")) pred4 <- NULL adult_predictions <- rbind(pred1n3, pred4, pred16) try({pred1n3_imm <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-1n3-imm-biomass-prior-", priors, ".rds" )) }) try({pred4_imm <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-4-imm-biomass-prior-", priors, ".rds")) }) try({pred16_imm <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-16-imm-biomass-prior-", priors, ".rds")) }) if (!exists("pred1n3_imm")) pred1n3_imm <- NULL if (!exists("pred16_imm")) pred16_imm <- NULL if (!exists("pred4_imm")) pred4_imm <- NULL imm_predictions <- rbind(pred1n3_imm, pred4_imm, pred16_imm) } } if (!exists("adult_predictions")) { paste("Maturity split biomass predictions are not available.") try({predictions <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-", ssid_string, "-total-biomass-prior-", priors, ".rds")) }) if (!exists("predictions")) { try({pred1n3_all <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-1n3-total-biomass-prior-", priors, ".rds")) }) try({pred4_all <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-4-total-biomass-prior-", priors, ".rds")) }) try({pred16_all <- readRDS(paste0("data/", spp, "/predictions-", spp, covs, "-16-total-biomass-prior-", priors, ".rds")) }) try({ predictions <- rbind(pred1n3_all, pred4_all, pred16_all) }) if (!exists("predictions")) paste("Total biomass predictions are not available.") } }
Caluclate variables of interest from predictions
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 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-biomass-prior-", priors, ".rds" )) saveRDS(imm_predictions, file = paste0( "data/", spp, "/all-predictions-", spp, covs, "-imm-biomass-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$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-biomass-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-biomass-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-biomass-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-biomass-priors-", priors, ".rds")) } else { out_total <- out_mature rm(out_mature) saveRDS(out_total, file = paste0("data/", spp, "/trends-", spp, covs, "-", ssid_string, "-total-biomass-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, legend_position = c(0.2, 0.25) ) + ggtitle(paste0( "", species, " ", min(unique(adult_predictions[["year"]])), "-", max(unique(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(predictions[["year"]])), "-", max(unique(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-biomass-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-biomass-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, "-priors-", priors, ".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, "/biomass-", spp, covs, "-", ssid_string, "-priors-", priors, ".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, "/exploratory-", spp, covs, "-", ssid_string, "-priors-", priors, ".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 <- plot_facet_map(adult_predictions, "epsilon_st", # white_zero = TRUE # #raster_limits = c(0, max_raster), # #legend_position = "none", # #transform_col = fourth_root_power # ) + # #labs(fill = "kg/ha") + # ggtitle(paste0("", species, "mature spatiotemporal variation")) # 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) }
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, "/biotic-residuals-", spp, covs, "-", ssid_string, "-priors-", priors, ".png"), res = 600, units = "in", width = 16, height = 8.5 ) gridExtra::grid.arrange( grobs = p_st, nrow = 2 ) dev.off()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.