In pbs-assess/gfvelocities: Spatiotemporal modelling of groundfish biomass densities in response to climate
Fit spatial temporal models to bottom temperature
library(tidyverse)
library(sdmTMB)
# library(future)
library(dplyr)
library(gfvelocities)
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 = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()
)
}
Model diagnostics function
get_diag <- function(m, response = "temperature_c", variable = "depth_scaled", colour_var = "depth", start_year = 2007, ssids = c(1,3,4,16)) {
predictions <- predict(m)
predictions$residuals <- residuals(m)
predictions <- predictions %>% filter (year >= start_year)
plot_map <- function(dat, column = "est") {
ggplot(dat, aes_string("X", "Y", colour = column)) +
geom_point() +
coord_fixed()+
gfplot::theme_pbs()
}
g <- plot_map(predictions, "est") +
scale_colour_viridis_c() +
ggtitle("Prediction (fixed effects + all random effects)")
print(g)
g <- plot_map(predictions, "est_non_rf") +
ggtitle("Prediction (fixed effects only)") +
scale_colour_viridis_c()
print(g)
g <- plot_map(predictions, "est_rf") +
ggtitle("All random effects only") +
scale_colour_gradient2()
print(g)
g <- plot_map(predictions, "omega_s") +
ggtitle("Spatial random effects only") +
scale_colour_gradient2()
print(g)
g <- plot_map(predictions, "epsilon_st") +
ggtitle("Spatiotemporal random effects only") +
facet_wrap(~year) +
scale_colour_gradient2()
print(g)
g <- ggplot(predictions, aes_string("est", response)) +
geom_point(alpha = 0.2) +
facet_wrap(~year) +
coord_fixed() +
geom_abline()+
gfplot::theme_pbs()
print(g)
g <- ggplot(predictions, aes(est, residuals)) +
geom_point(alpha = 0.2) +
geom_smooth()+
gfplot::theme_pbs()
print(g)
g <- ggplot(predictions, aes_string(variable, "residuals", colour = colour_var)) +
geom_point(alpha = 0.4, size = 1.2) +
geom_smooth(colour="grey", size = 1.2) +
scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7,
limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt
# facet_wrap(~year, scales = "free_x")+
gfplot::theme_pbs() + theme(
axis.text.y = element_text(size = 14),
axis.text.x = element_blank(), axis.ticks.x = element_blank())
print(g)
g <- ggplot(predictions, aes_string(variable, "residuals", colour = colour_var)) +
geom_point(alpha = 0.4, size = 0.7) +
geom_smooth(colour="grey", size = 0.7) +
scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7,
limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt
facet_wrap(~year, scales = "free_x")+
gfplot::theme_pbs() + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
print(g)
g <- ggplot(filter(predictions, ssid %in% ssids), aes_string(variable, "residuals", colour = colour_var)) +
geom_point(alpha = 0.4, size = 1) +
geom_smooth(colour="grey", size = 1) +
scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7,
limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt
facet_wrap(~ssid, scales = "free_x")+
gfplot::theme_pbs()
print(g)
g <- ggplot(predictions, aes(X, Y, colour = residuals)) +
geom_point(alpha = 0.2) +
coord_fixed() +
scale_color_gradient2() +
facet_wrap(~year)+
gfplot::theme_pbs()
print(g)
aic <- function(m) {
k <- length(m$model$par)
nll <- m$model$objective
2 * k - 2 * (-nll)
}
print("R^2:")
r2 <- cor(predictions$est, predictions[[response]])^2
print(r2)
print("")
print("AIC:")
print(aic(m))
}
Load data if not using SST
all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-processed.rds")) # no SST
all_depth <- all_sensor %>%
# filter(depth_max>10) %>% # maybe we should filter incase missing depth values also indicate a sensor fail?
# filter missing location data and trial year
dplyr::filter(!is.na(latitude), !is.na(longitude)) %>%
# convert lat and lon to UTMs
dplyr::mutate(X = longitude, Y = latitude) %>%
gfplot:::ll2utm(., utm_zone = 9) %>%
# interpolate missing depth values
dplyr::rename(depth = depth_m) %>%
gfplot:::interp_survey_bathymetry()
d_trawl <- all_depth$data %>%
# remove obvious sensor fails, d_trawl %>% filter(depth_max<10), fishing event ids = 481861, 2179089
filter(depth_max>10) %>%
gfplot:::scale_survey_predictors()
d_trawl <- d_trawl[!is.na(d_trawl$depth), ]
d_trawl <- d_trawl[!is.na(d_trawl$temperature_c), ]
d_trawl <- mutate(d_trawl,
DOY = lubridate::yday(date),
DOY_scaled = arm::rescale(lubridate::yday(date)),
DOY_scaled2 = arm::rescale(lubridate::yday(date))^2)
# 1 redbanded fish found at too shallow a depth... probably leftover in net?
# slender sole also has an outlier here, so maybe should throw out this event for everyone?
# could mean there's a location or depth error? currently still in climate models
# d_trawl <- d_trawl[d_trawl$fishing_event_id!=2536031,]
EXLORE TEMPERATURE at DEPTH
ggplot(d_trawl, aes_string("X", "Y", colour = "temperature_c")) +
geom_point() +
facet_wrap(~year) +
coord_fixed() +
scale_color_viridis_c()
ggplot(d_trawl, aes(DOY, temperature_c, colour = -depth, alpha=0.5)) + geom_point() + facet_wrap(~year)
Current temperature model
spde <- sdmTMB::make_mesh(d_trawl, c("X", "Y"), n_knots = 800, type = "kmeans")
plot(spde)
# spde <- make_spde(d_trawl$X, d_trawl$Y, n_knots = 800)
m_temp_rw <- sdmTMB(
temperature_c ~ 0 + as.factor(year),
time_varying = ~ 0 + depth_scaled + depth_scaled2,
data = d_trawl,
ar1_fields = TRUE, # changed AR1 to true
include_spatial = TRUE,
time = "year", spde = spde,
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = FALSE
)
saveRDS(m_temp_rw, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-2.rds")))
New CV model code
# m_temp_cv <- sdmTMB_cv(
# temperature_c ~ 0 + as.factor(year),
# time_varying = ~ 0 + depth_scaled + depth_scaled2,
# data = d_trawl,
# ar1_fields = TRUE, # changed AR1 to true
# include_spatial = TRUE,
# time = "year", spde = spde,
# family = gaussian(link = "identity"),
# k_folds = 5,
# # nlminb_loops = 2,
# # newton_steps = 1,
# silent = FALSE
# )
# saveRDS(m_temp_cv, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-cv.rds")))
m_temp_cv <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-cv.rds")))
m_temp_cv$fold_loglik
m_temp_cv$sum_loglik
m_temp_cv$max_gradients
Check predictions against raw data
m_temp_rw <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-2.rds"))
# m_temp_rw <- sdmTMB:::update_model(m_temp_rw)
predictions <- predict(m_temp_rw)
predictions$residuals <- residuals(m_temp_rw)
cor(predictions$est, predictions$temperature_c)^2
sum(abs(predictions$residuals))
saveRDS(predictions, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-predictions.rds")))
# 2908.879-2789.322
# 119.557/2908.879 # 4% reduction in abs errors
# test against ROMS
events_roms <- readRDS("data/roms_temp_by_events2.rds") %>% select(fishing_event_id, roms, year)
predictions <- left_join(predictions, events_roms)
cor(predictions$roms, predictions$est, use = "pairwise.complete.obs")
cor(predictions$roms, predictions$temperature_c, use = "pairwise.complete.obs")
ggplot(predictions, aes(temperature_c, est, colour = year)) + geom_point(alpha=0.3) +
# facet_wrap(~year) +
coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 1) +
scale_colour_viridis_c(option = "A", end = 0.75) +
gfplot::theme_pbs()
ggplot(predictions, aes(roms, est, colour = -depth)) + geom_point(alpha=0.3) +
facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 1) +
gfplot::theme_pbs()
ggplot(predictions, aes(roms, temperature_c, colour = DOY)) + geom_point(alpha=0.3) +
facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 1) +
scale_colour_viridis_c() +
gfplot::theme_pbs()
ggplot(predictions, aes(DOY, (roms-temperature_c), colour = -depth)) + geom_point(alpha=0.3) +
facet_wrap(~year) + #coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 0) +
scale_colour_viridis_c() +
gfplot::theme_pbs()
ggplot(predictions, aes(DOY, (est-roms), colour = -depth)) + geom_point(alpha=0.3) +
facet_wrap(~year) + #coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 0) +
scale_colour_viridis_c() +
gfplot::theme_pbs()
p_nd <- readRDS((here::here("analysis/VOCC/data/predicted-DO-2020-06-20-more2016-roms.rds")))
cor(p_nd$roms_temp, p_nd$temp, use = "pairwise.complete.obs")
ggplot(p_nd, aes(roms_temp, temp, colour = -depth)) + geom_point(alpha=0.04) +
facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 1) +
gfplot::theme_pbs()
get_diag(m_temp_rw, response = "temperature_c", variable = "depth_scaled")
WOULD SST HELP?
Load and prep SST data
all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-meanSST-all.rds"))
d_trawl2 <-all_sensor
# all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-processed.rds")) # no SST
# glimpse(all_sensor)
#
# all_depth <- all_sensor %>%
# # filter(depth_max>10) %>% # maybe we should filter incase missing depth values also indicate a sensor fail?
# # filter missing location data and trial year
# dplyr::filter(!is.na(latitude), !is.na(longitude)) %>%
#
# # convert lat and lon to UTMs
# dplyr::mutate(X = longitude, Y = latitude) %>%
# gfplot:::ll2utm(., utm_zone = 9) %>%
#
# # interpolate missing depth values
# dplyr::rename(depth = depth_m) %>%
# gfplot:::interp_survey_bathymetry()
# d_trawl <- all_depth$data %>%
# # remove obvious sensor fails, d_trawl %>% filter(depth_max<10), fishing event ids = 481861, 2179089
# filter(depth_max>10) %>%
# gfplot:::scale_survey_predictors()
d_trawl2 <- d_trawl2[!is.na(d_trawl2$SST), ]
d_trawl2 <- d_trawl2[!is.na(d_trawl2$meanSST), ]
d_trawl2 <- mutate(d_trawl2,
depth = depth_m,
DOY = lubridate::yday(date),
DOY_scaled = arm::rescale(DOY),
DOY_scaled2 = DOY_scaled^2,
depth_scaled = arm::rescale(log(depth)),
SST_scaled = arm::rescale(meanSST),
depth_scaled2 = depth_scaled^2,
shallow = ifelse(depth>35,0,1),
deep = ifelse(depth<35,0,1)
)
d_trawl2 <- d_trawl2[!is.na(d_trawl2$depth_scaled), ]
# d_trawl2[is.na(d_trawl2$shallow), ]
# d_trawl2[is.na(d_trawl2$SST_scaled), ]
Explore SST on survey day
ggplot(d_trawl2, aes_string("X", "Y", colour = "SST")) +
geom_point() +
facet_wrap(~year) +
coord_fixed() +
scale_color_viridis_c()
Is bottom temp related to SST in shallower samples?
ggplot(dplyr::filter(d_trawl2, depth <35), aes(SST,temperature_c, colour = -depth)) + geom_point()
ggplot(dplyr::filter(d_trawl2, depth <35), aes(SST,temperature_c, colour = DOY)) + geom_point()
Explore mean SST
ggplot(d_trawl2, aes_string("X", "Y", colour = "meanSST")) +
geom_point() +
facet_wrap(~year) +
coord_fixed() +
scale_color_viridis_c()
ggplot(d_trawl2, aes_string("depth_m", "meanSST", colour = "meanSST")) +
geom_point() +
facet_wrap(~ssid, scales = "free") +
# coord_fixed() +
scale_color_viridis_c()
ggplot(dplyr::filter(d_trawl2, shallow==1), aes(SST_scaled,temperature_c, colour = -depth)) + geom_point()
Try BOOSTED REGRESSION MODEL to check for non-linearity
library(gbm)
d_trawl2$log_depth <- log(d_trawl2$depth_m)
m <- gbm( temperature_c ~ log_depth + DOY + X + Y, data = d_trawl2, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02)
m2 <- gbm( temperature_c ~ SST_scaled + log_depth + DOY + X + Y, data = d_trawl2, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02)
m3 <- gbm( temperature_c ~ SST_scaled + depth_m + DOY + X + Y, data = d_trawl2, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02)
plot(m,i.var=1) # exp(3.4)= ~30 m is when depth effect begins
plot(m,i.var=2) # strange pattern
plot(m,i.var=3)
plot(m,i.var=1:2) # 50 m depth is when effect of DOY disappears entirely
plot(m2,i.var=1)
plot(m2,i.var=2)
plot(m2,i.var=3)
plot(m2,i.var=4)
plot(m2,i.var=c(2,1))
plot(m2,i.var=c(2,3)) # only < 160 DOY show cooler temps at <30m depths
plot(m2,i.var=c(1,3))
plot(m3,i.var=1:2)
d_trawl2$r <- predict(m, n.trees = 2000) - d_trawl2$temperature_c
ggplot(d_trawl2, aes(depth_scaled, r)) +
geom_point(alpha=0.4) +
ylim(-8,8) +
geom_smooth()
# complex interaction with depth could prove useful at depths shallower than 35m, but for vast majority of depths sampled, DOY and SST should have little effect so test simplier linear models below
complex interaction with depth could prove useful at depths shallower than 35m, but for vast majority of depths sampled, DOY and SST should have little effect so test simplier linear models below
Rerun main model with SST dataset for comparisons
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800)
m_temp_new <- sdmTMB(
temperature_c ~ 0 + as.factor(year),
time_varying = ~ 0 + depth_scaled + depth_scaled2,
data = d_trawl2, spde = spde,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = FALSE
)
saveRDS(m_temp_new, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-sstdata.rds")))
m_temp_new
AIC(m_temp_new)
try just DOY
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800)
m_temp_doy_only <- sdmTMB(
temperature_c ~ 0 + as.factor(year) + DOY_scaled,
data = d_trawl2, spde = spde,
time_varying = ~ 0 + depth_scaled + depth_scaled2,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = TRUE
)
saveRDS(m_temp_doy_only, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-DOY-only.rds")))
tidy(m_temp_doy_only, conf.int = T)
try DOY and SST
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800)
m_temp_sst <- sdmTMB(
temperature_c ~ 0 + as.factor(year) + DOY_scaled + SST_scaled,
data = d_trawl2, spde = spde,
time_varying = ~ 0 + depth_scaled + depth_scaled2,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = TRUE
)
saveRDS(m_temp_sst, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-DOY-only.rds")))
tidy(m_temp_sst, conf.int = T)
compare models
m_temp_new <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-sstdata.rds")))
m_temp_doy_only <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-DOY-only.rds")))
m_temp_sst <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-SST.rds")))
AIC(m_temp_new, m_temp_doy_only, m_temp_sst)
Plotting predictions from depth only temperature model
nd_all <- readRDS(here::here("analysis/VOCC/data/nd_all_synoptic.rds")) %>% filter(year<2019)
unique(nd_all$year)
glimpse(nd_all)
p_nd <- predict(m_temp_rw, newdata = nd_all)
# saveRDS(p_nd, file = "../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds")
plot_map_raster(p_nd, "est") +
scale_fill_viridis_c(trans = sqrt) +
facet_wrap(~year) +
ggtitle("Prediction (fixed effects + all random effects)")
plot_facet_map(p_nd, "est", viridis_option = "D",
raster_limits = c(min(p_nd$est), quantile(p_nd$est, 0.99)),
transform_col = sqrt
)
p_nd1 <- p_nd %>% filter(ssid!=16) %>% filter(ssid!=4)
breaks <- scales::trans_breaks(no_trans[["transform"]], no_trans[["inverse"]],n = 6)
labels <- function(x) format(x, digits = 2, scientific = FALSE)
pred_1 <- plot_map_raster(p_nd1, "est") +
scale_fill_viridis_c(labels = labels,
option = "C",
limits = c(4, 11),
na.value = "red") +
facet_wrap(~year) +
labs(fill = "°C") +
ggtitle("Prediction (fixed effects + all random effects)")
pred_1
Time-varying depth
# pd <- expand.grid(
# depth_scaled = seq(min(m_temp_rw$data$depth_scaled),
# max(m_temp_rw$data$depth_scaled), length.out = 50),
# year = unique(m_temp_rw$data$year)
# )
# pd$depth_scaled2 <- pd$depth_scaled^2
#
# # p <- predict(m_temp_rw, newdata = pd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))
# p <- predict(m_temp_new, newdata = pd, se_fit = TRUE, re_form = NA)
#
# p$Depth <- exp(p$depth_scaled*0.661 + 5.07)
# p$Temperature <- p$est
# #
# p$lowCI <- (p$est - 1.96 * p$est_se)
# p$highCI <- (p$est + 1.96 * p$est_se)
# #
# saveRDS(p, here::here("analysis/tmb-sensor-explore/data/depth-tv-m-temp-rw-800-kn.rds"))
p <- readRDS(here::here("analysis/tmb-sensor-explore/data/depth-tv-m-temp-rw-800-kn.rds"))
ggplot(filter(p, Depth >15 & year >2007), aes(Depth, Temperature,
ymin = lowCI,
ymax = highCI,
group = as.factor(year), fill=year, colour = year)) +
geom_line(size = 0.5, alpha =0.85) +
# geom_smooth(method= "loess", size = 0.4, alpha =0.85, se =F ) +
# geom_smooth(span = 0.5, size = 0.4, alpha =0.85, se =F ) +
scale_fill_viridis_c(option = "C") +
scale_colour_viridis_c(option = "C") +
# ylab("DO") +
geom_ribbon(alpha = 0.2, colour = NA) +
# scale_y_continuous(limits = c(0,4.5)) +
scale_x_reverse(limits = c(600,20)) +
ylim( 2, 15) +
coord_flip() +
# scale_x_continuous(limits = c(18,210)) +
gfplot::theme_pbs() +
theme(legend.position = c(.8,.4))
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/depth-tv-m-temp-rw-800-kn.png"), width = 3, height = 2.5)
# ggplot(filter(p, Depth >15 & year == 2018), aes(Depth, Temperature,
# ymin = (Temperature - 1.96 * est_se),
# ymax = (Temperature + 1.96 * est_se))) +
# geom_line(size = 1.5, alpha =0.85) +
# # scale_color_viridis_c(option = "C") +
# geom_ribbon(alpha=0.2) +
# scale_x_reverse(limits = c(210,20)) +
# ylim(5, 15) +
# coord_flip() +
# gfplot::theme_pbs() +
# theme(legend.position = c(.8,.4))
p_nd4 <- p_nd %>% filter(ssid==4)
pred_4 <- plot_map_raster(p_nd4, "est") +
scale_fill_viridis_c(labels = labels,
option = "C",
limits = c(5, 11),
na.value = "black")+
facet_wrap(~year) +
theme(legend.position = c(.85,.12)) +
labs(fill = "°C")
pred_4
plot_map_raster(p_nd, "est_non_rf") +
ggtitle("Prediction (fixed effects only)") +
facet_wrap(~year) +
scale_fill_viridis_c()
plot_map_raster(p_nd, "est_rf") +
ggtitle("Prediction (all random effects only)") +
facet_wrap(~year) +
scale_fill_viridis_c()
(omega_s_all <- plot_map_raster(p_nd, "omega_s") +
ggtitle("Spatial random effects only") +
scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5),
na.value = "red") + labs(fill = "") + theme(legend.position = c(0.2, 0.3))
)
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/omegas-temp-rw-800-kn.png"), width = 3, height = 3)
omega_s_all
plot_map_raster(p_nd, "epsilon_st") +
ggtitle("Spatiotemporal random effects only") +
facet_wrap(~year) + scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5),
na.value = "red")
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd1 <- p_nd %>% filter(ssid!=16) %>% filter(ssid!=4) %>% filter(year>2007)
(st_1n3 <- plot_map_raster(p_nd1, "epsilon_st") +
ggtitle("Spatiotemporal random effects") +
facet_grid(~year) +
scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5), #trans = sqrt,
na.value = "red") +
labs(fill = "epsilon") +
theme(
legend.position = "none"))
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-1n3.png"), width =6, height = 2)
# png(
# file = "figs/temp-predictions-1n3.png",
# res = 400,
# units = "in",
# width = 8,
# height = 8
# )
# st_1n3
# dev.off()
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd4 <- p_nd %>% filter(ssid==4) %>% filter(year>2007)#%>% filter(ssid!=4)
(st_4 <- plot_map_raster(p_nd4, "epsilon_st") +
#ggtitle(" ") +
facet_grid(~year) +
scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5), #trans = sqrt,
na.value = "red") +
labs(fill = " ") + theme( legend.position = "none") )
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-4.png"), width = 6, height = 2)
st_4
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd16 <- p_nd %>% filter(ssid==16) %>% filter(year>2007)#%>% filter(ssid!=4)
(st_16 <- plot_map_raster(p_nd16, "epsilon_st") +
#ggtitle(" ") +
facet_grid(~year) +
scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5),
na.value = "red") +
xlim(20,360) + ylim(5850,6050)+
# xlim(65,345) + ylim(5800,6100)+
theme(panel.border = element_blank(), legend.position = "none") )
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-16.png"), width = 6, height = 1)
st_16
png(
file = "figs/temp-predictions.png",
res = 600,
units = "in",
width = 11,
height = 10
)
gridExtra::grid.arrange(
grobs = list(pred_1, st_1n3, omega_s_all, pred_4, st_4, st_16 ),
widths = c(2, 2, 1.75),
heights = c(2, 1.35),
layout_matrix = rbind(c(1, 2, 3),
c(4, 5, 6)),
top = grid::textGrob("Bottom temperature estimates from sdmTMB model")
)
# gridExtra::grid.arrange( nrow = 2, heights = c(2,1.5),
# top = grid::textGrob("Dissolved oxygen estimates from sdmTMB model"))
dev.off()
pbs-assess/gfvelocities documentation built on Dec. 22, 2021, 6:45 a.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.
Fit spatial temporal models to bottom temperature
library(tidyverse) library(sdmTMB) # library(future) library(dplyr) library(gfvelocities) 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 = element_blank(), axis.ticks = element_blank(), axis.text = element_blank() ) }
Model diagnostics function
get_diag <- function(m, response = "temperature_c", variable = "depth_scaled", colour_var = "depth", start_year = 2007, ssids = c(1,3,4,16)) { predictions <- predict(m) predictions$residuals <- residuals(m) predictions <- predictions %>% filter (year >= start_year) plot_map <- function(dat, column = "est") { ggplot(dat, aes_string("X", "Y", colour = column)) + geom_point() + coord_fixed()+ gfplot::theme_pbs() } g <- plot_map(predictions, "est") + scale_colour_viridis_c() + ggtitle("Prediction (fixed effects + all random effects)") print(g) g <- plot_map(predictions, "est_non_rf") + ggtitle("Prediction (fixed effects only)") + scale_colour_viridis_c() print(g) g <- plot_map(predictions, "est_rf") + ggtitle("All random effects only") + scale_colour_gradient2() print(g) g <- plot_map(predictions, "omega_s") + ggtitle("Spatial random effects only") + scale_colour_gradient2() print(g) g <- plot_map(predictions, "epsilon_st") + ggtitle("Spatiotemporal random effects only") + facet_wrap(~year) + scale_colour_gradient2() print(g) g <- ggplot(predictions, aes_string("est", response)) + geom_point(alpha = 0.2) + facet_wrap(~year) + coord_fixed() + geom_abline()+ gfplot::theme_pbs() print(g) g <- ggplot(predictions, aes(est, residuals)) + geom_point(alpha = 0.2) + geom_smooth()+ gfplot::theme_pbs() print(g) g <- ggplot(predictions, aes_string(variable, "residuals", colour = colour_var)) + geom_point(alpha = 0.4, size = 1.2) + geom_smooth(colour="grey", size = 1.2) + scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7, limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt # facet_wrap(~year, scales = "free_x")+ gfplot::theme_pbs() + theme( axis.text.y = element_text(size = 14), axis.text.x = element_blank(), axis.ticks.x = element_blank()) print(g) g <- ggplot(predictions, aes_string(variable, "residuals", colour = colour_var)) + geom_point(alpha = 0.4, size = 0.7) + geom_smooth(colour="grey", size = 0.7) + scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7, limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt facet_wrap(~year, scales = "free_x")+ gfplot::theme_pbs() + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) print(g) g <- ggplot(filter(predictions, ssid %in% ssids), aes_string(variable, "residuals", colour = colour_var)) + geom_point(alpha = 0.4, size = 1) + geom_smooth(colour="grey", size = 1) + scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7, limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt facet_wrap(~ssid, scales = "free_x")+ gfplot::theme_pbs() print(g) g <- ggplot(predictions, aes(X, Y, colour = residuals)) + geom_point(alpha = 0.2) + coord_fixed() + scale_color_gradient2() + facet_wrap(~year)+ gfplot::theme_pbs() print(g) aic <- function(m) { k <- length(m$model$par) nll <- m$model$objective 2 * k - 2 * (-nll) } print("R^2:") r2 <- cor(predictions$est, predictions[[response]])^2 print(r2) print("") print("AIC:") print(aic(m)) }
Load data if not using SST
all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-processed.rds")) # no SST all_depth <- all_sensor %>% # filter(depth_max>10) %>% # maybe we should filter incase missing depth values also indicate a sensor fail? # filter missing location data and trial year dplyr::filter(!is.na(latitude), !is.na(longitude)) %>% # convert lat and lon to UTMs dplyr::mutate(X = longitude, Y = latitude) %>% gfplot:::ll2utm(., utm_zone = 9) %>% # interpolate missing depth values dplyr::rename(depth = depth_m) %>% gfplot:::interp_survey_bathymetry() d_trawl <- all_depth$data %>% # remove obvious sensor fails, d_trawl %>% filter(depth_max<10), fishing event ids = 481861, 2179089 filter(depth_max>10) %>% gfplot:::scale_survey_predictors() d_trawl <- d_trawl[!is.na(d_trawl$depth), ] d_trawl <- d_trawl[!is.na(d_trawl$temperature_c), ] d_trawl <- mutate(d_trawl, DOY = lubridate::yday(date), DOY_scaled = arm::rescale(lubridate::yday(date)), DOY_scaled2 = arm::rescale(lubridate::yday(date))^2) # 1 redbanded fish found at too shallow a depth... probably leftover in net? # slender sole also has an outlier here, so maybe should throw out this event for everyone? # could mean there's a location or depth error? currently still in climate models # d_trawl <- d_trawl[d_trawl$fishing_event_id!=2536031,]
EXLORE TEMPERATURE at DEPTH
ggplot(d_trawl, aes_string("X", "Y", colour = "temperature_c")) + geom_point() + facet_wrap(~year) + coord_fixed() + scale_color_viridis_c()
ggplot(d_trawl, aes(DOY, temperature_c, colour = -depth, alpha=0.5)) + geom_point() + facet_wrap(~year)
Current temperature model
spde <- sdmTMB::make_mesh(d_trawl, c("X", "Y"), n_knots = 800, type = "kmeans") plot(spde)
# spde <- make_spde(d_trawl$X, d_trawl$Y, n_knots = 800) m_temp_rw <- sdmTMB( temperature_c ~ 0 + as.factor(year), time_varying = ~ 0 + depth_scaled + depth_scaled2, data = d_trawl, ar1_fields = TRUE, # changed AR1 to true include_spatial = TRUE, time = "year", spde = spde, family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = FALSE ) saveRDS(m_temp_rw, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-2.rds")))
New CV model code
# m_temp_cv <- sdmTMB_cv( # temperature_c ~ 0 + as.factor(year), # time_varying = ~ 0 + depth_scaled + depth_scaled2, # data = d_trawl, # ar1_fields = TRUE, # changed AR1 to true # include_spatial = TRUE, # time = "year", spde = spde, # family = gaussian(link = "identity"), # k_folds = 5, # # nlminb_loops = 2, # # newton_steps = 1, # silent = FALSE # ) # saveRDS(m_temp_cv, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-cv.rds"))) m_temp_cv <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-cv.rds"))) m_temp_cv$fold_loglik m_temp_cv$sum_loglik m_temp_cv$max_gradients
Check predictions against raw data
m_temp_rw <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-2.rds")) # m_temp_rw <- sdmTMB:::update_model(m_temp_rw) predictions <- predict(m_temp_rw) predictions$residuals <- residuals(m_temp_rw) cor(predictions$est, predictions$temperature_c)^2 sum(abs(predictions$residuals)) saveRDS(predictions, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-predictions.rds"))) # 2908.879-2789.322 # 119.557/2908.879 # 4% reduction in abs errors
# test against ROMS events_roms <- readRDS("data/roms_temp_by_events2.rds") %>% select(fishing_event_id, roms, year) predictions <- left_join(predictions, events_roms) cor(predictions$roms, predictions$est, use = "pairwise.complete.obs") cor(predictions$roms, predictions$temperature_c, use = "pairwise.complete.obs") ggplot(predictions, aes(temperature_c, est, colour = year)) + geom_point(alpha=0.3) + # facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 1) + scale_colour_viridis_c(option = "A", end = 0.75) + gfplot::theme_pbs() ggplot(predictions, aes(roms, est, colour = -depth)) + geom_point(alpha=0.3) + facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 1) + gfplot::theme_pbs() ggplot(predictions, aes(roms, temperature_c, colour = DOY)) + geom_point(alpha=0.3) + facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 1) + scale_colour_viridis_c() + gfplot::theme_pbs() ggplot(predictions, aes(DOY, (roms-temperature_c), colour = -depth)) + geom_point(alpha=0.3) + facet_wrap(~year) + #coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 0) + scale_colour_viridis_c() + gfplot::theme_pbs() ggplot(predictions, aes(DOY, (est-roms), colour = -depth)) + geom_point(alpha=0.3) + facet_wrap(~year) + #coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 0) + scale_colour_viridis_c() + gfplot::theme_pbs() p_nd <- readRDS((here::here("analysis/VOCC/data/predicted-DO-2020-06-20-more2016-roms.rds"))) cor(p_nd$roms_temp, p_nd$temp, use = "pairwise.complete.obs") ggplot(p_nd, aes(roms_temp, temp, colour = -depth)) + geom_point(alpha=0.04) + facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 1) + gfplot::theme_pbs()
get_diag(m_temp_rw, response = "temperature_c", variable = "depth_scaled")
WOULD SST HELP?
Load and prep SST data
all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-meanSST-all.rds")) d_trawl2 <-all_sensor # all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-processed.rds")) # no SST # glimpse(all_sensor) # # all_depth <- all_sensor %>% # # filter(depth_max>10) %>% # maybe we should filter incase missing depth values also indicate a sensor fail? # # filter missing location data and trial year # dplyr::filter(!is.na(latitude), !is.na(longitude)) %>% # # # convert lat and lon to UTMs # dplyr::mutate(X = longitude, Y = latitude) %>% # gfplot:::ll2utm(., utm_zone = 9) %>% # # # interpolate missing depth values # dplyr::rename(depth = depth_m) %>% # gfplot:::interp_survey_bathymetry() # d_trawl <- all_depth$data %>% # # remove obvious sensor fails, d_trawl %>% filter(depth_max<10), fishing event ids = 481861, 2179089 # filter(depth_max>10) %>% # gfplot:::scale_survey_predictors() d_trawl2 <- d_trawl2[!is.na(d_trawl2$SST), ] d_trawl2 <- d_trawl2[!is.na(d_trawl2$meanSST), ] d_trawl2 <- mutate(d_trawl2, depth = depth_m, DOY = lubridate::yday(date), DOY_scaled = arm::rescale(DOY), DOY_scaled2 = DOY_scaled^2, depth_scaled = arm::rescale(log(depth)), SST_scaled = arm::rescale(meanSST), depth_scaled2 = depth_scaled^2, shallow = ifelse(depth>35,0,1), deep = ifelse(depth<35,0,1) ) d_trawl2 <- d_trawl2[!is.na(d_trawl2$depth_scaled), ] # d_trawl2[is.na(d_trawl2$shallow), ] # d_trawl2[is.na(d_trawl2$SST_scaled), ]
Explore SST on survey day
ggplot(d_trawl2, aes_string("X", "Y", colour = "SST")) + geom_point() + facet_wrap(~year) + coord_fixed() + scale_color_viridis_c()
Is bottom temp related to SST in shallower samples?
ggplot(dplyr::filter(d_trawl2, depth <35), aes(SST,temperature_c, colour = -depth)) + geom_point() ggplot(dplyr::filter(d_trawl2, depth <35), aes(SST,temperature_c, colour = DOY)) + geom_point()
Explore mean SST
ggplot(d_trawl2, aes_string("X", "Y", colour = "meanSST")) + geom_point() + facet_wrap(~year) + coord_fixed() + scale_color_viridis_c()
ggplot(d_trawl2, aes_string("depth_m", "meanSST", colour = "meanSST")) + geom_point() + facet_wrap(~ssid, scales = "free") + # coord_fixed() + scale_color_viridis_c()
ggplot(dplyr::filter(d_trawl2, shallow==1), aes(SST_scaled,temperature_c, colour = -depth)) + geom_point()
Try BOOSTED REGRESSION MODEL to check for non-linearity
library(gbm) d_trawl2$log_depth <- log(d_trawl2$depth_m) m <- gbm( temperature_c ~ log_depth + DOY + X + Y, data = d_trawl2, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02) m2 <- gbm( temperature_c ~ SST_scaled + log_depth + DOY + X + Y, data = d_trawl2, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02) m3 <- gbm( temperature_c ~ SST_scaled + depth_m + DOY + X + Y, data = d_trawl2, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02) plot(m,i.var=1) # exp(3.4)= ~30 m is when depth effect begins plot(m,i.var=2) # strange pattern plot(m,i.var=3) plot(m,i.var=1:2) # 50 m depth is when effect of DOY disappears entirely plot(m2,i.var=1) plot(m2,i.var=2) plot(m2,i.var=3) plot(m2,i.var=4) plot(m2,i.var=c(2,1)) plot(m2,i.var=c(2,3)) # only < 160 DOY show cooler temps at <30m depths plot(m2,i.var=c(1,3)) plot(m3,i.var=1:2) d_trawl2$r <- predict(m, n.trees = 2000) - d_trawl2$temperature_c ggplot(d_trawl2, aes(depth_scaled, r)) + geom_point(alpha=0.4) + ylim(-8,8) + geom_smooth() # complex interaction with depth could prove useful at depths shallower than 35m, but for vast majority of depths sampled, DOY and SST should have little effect so test simplier linear models below
complex interaction with depth could prove useful at depths shallower than 35m, but for vast majority of depths sampled, DOY and SST should have little effect so test simplier linear models below
Rerun main model with SST dataset for comparisons
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800) m_temp_new <- sdmTMB( temperature_c ~ 0 + as.factor(year), time_varying = ~ 0 + depth_scaled + depth_scaled2, data = d_trawl2, spde = spde, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = FALSE ) saveRDS(m_temp_new, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-sstdata.rds"))) m_temp_new AIC(m_temp_new)
try just DOY
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800) m_temp_doy_only <- sdmTMB( temperature_c ~ 0 + as.factor(year) + DOY_scaled, data = d_trawl2, spde = spde, time_varying = ~ 0 + depth_scaled + depth_scaled2, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = TRUE ) saveRDS(m_temp_doy_only, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-DOY-only.rds"))) tidy(m_temp_doy_only, conf.int = T)
try DOY and SST
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800) m_temp_sst <- sdmTMB( temperature_c ~ 0 + as.factor(year) + DOY_scaled + SST_scaled, data = d_trawl2, spde = spde, time_varying = ~ 0 + depth_scaled + depth_scaled2, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = TRUE ) saveRDS(m_temp_sst, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-DOY-only.rds"))) tidy(m_temp_sst, conf.int = T)
compare models
m_temp_new <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-sstdata.rds"))) m_temp_doy_only <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-DOY-only.rds"))) m_temp_sst <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-SST.rds"))) AIC(m_temp_new, m_temp_doy_only, m_temp_sst)
Plotting predictions from depth only temperature model
nd_all <- readRDS(here::here("analysis/VOCC/data/nd_all_synoptic.rds")) %>% filter(year<2019) unique(nd_all$year) glimpse(nd_all) p_nd <- predict(m_temp_rw, newdata = nd_all) # saveRDS(p_nd, file = "../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds") plot_map_raster(p_nd, "est") + scale_fill_viridis_c(trans = sqrt) + facet_wrap(~year) + ggtitle("Prediction (fixed effects + all random effects)")
plot_facet_map(p_nd, "est", viridis_option = "D", raster_limits = c(min(p_nd$est), quantile(p_nd$est, 0.99)), transform_col = sqrt )
p_nd1 <- p_nd %>% filter(ssid!=16) %>% filter(ssid!=4) breaks <- scales::trans_breaks(no_trans[["transform"]], no_trans[["inverse"]],n = 6) labels <- function(x) format(x, digits = 2, scientific = FALSE) pred_1 <- plot_map_raster(p_nd1, "est") + scale_fill_viridis_c(labels = labels, option = "C", limits = c(4, 11), na.value = "red") + facet_wrap(~year) + labs(fill = "°C") + ggtitle("Prediction (fixed effects + all random effects)") pred_1
Time-varying depth
# pd <- expand.grid( # depth_scaled = seq(min(m_temp_rw$data$depth_scaled), # max(m_temp_rw$data$depth_scaled), length.out = 50), # year = unique(m_temp_rw$data$year) # ) # pd$depth_scaled2 <- pd$depth_scaled^2 # # # p <- predict(m_temp_rw, newdata = pd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y")) # p <- predict(m_temp_new, newdata = pd, se_fit = TRUE, re_form = NA) # # p$Depth <- exp(p$depth_scaled*0.661 + 5.07) # p$Temperature <- p$est # # # p$lowCI <- (p$est - 1.96 * p$est_se) # p$highCI <- (p$est + 1.96 * p$est_se) # # # saveRDS(p, here::here("analysis/tmb-sensor-explore/data/depth-tv-m-temp-rw-800-kn.rds")) p <- readRDS(here::here("analysis/tmb-sensor-explore/data/depth-tv-m-temp-rw-800-kn.rds")) ggplot(filter(p, Depth >15 & year >2007), aes(Depth, Temperature, ymin = lowCI, ymax = highCI, group = as.factor(year), fill=year, colour = year)) + geom_line(size = 0.5, alpha =0.85) + # geom_smooth(method= "loess", size = 0.4, alpha =0.85, se =F ) + # geom_smooth(span = 0.5, size = 0.4, alpha =0.85, se =F ) + scale_fill_viridis_c(option = "C") + scale_colour_viridis_c(option = "C") + # ylab("DO") + geom_ribbon(alpha = 0.2, colour = NA) + # scale_y_continuous(limits = c(0,4.5)) + scale_x_reverse(limits = c(600,20)) + ylim( 2, 15) + coord_flip() + # scale_x_continuous(limits = c(18,210)) + gfplot::theme_pbs() + theme(legend.position = c(.8,.4)) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/depth-tv-m-temp-rw-800-kn.png"), width = 3, height = 2.5) # ggplot(filter(p, Depth >15 & year == 2018), aes(Depth, Temperature, # ymin = (Temperature - 1.96 * est_se), # ymax = (Temperature + 1.96 * est_se))) + # geom_line(size = 1.5, alpha =0.85) + # # scale_color_viridis_c(option = "C") + # geom_ribbon(alpha=0.2) + # scale_x_reverse(limits = c(210,20)) + # ylim(5, 15) + # coord_flip() + # gfplot::theme_pbs() + # theme(legend.position = c(.8,.4))
p_nd4 <- p_nd %>% filter(ssid==4) pred_4 <- plot_map_raster(p_nd4, "est") + scale_fill_viridis_c(labels = labels, option = "C", limits = c(5, 11), na.value = "black")+ facet_wrap(~year) + theme(legend.position = c(.85,.12)) + labs(fill = "°C") pred_4
plot_map_raster(p_nd, "est_non_rf") + ggtitle("Prediction (fixed effects only)") + facet_wrap(~year) + scale_fill_viridis_c()
plot_map_raster(p_nd, "est_rf") + ggtitle("Prediction (all random effects only)") + facet_wrap(~year) + scale_fill_viridis_c()
(omega_s_all <- plot_map_raster(p_nd, "omega_s") + ggtitle("Spatial random effects only") + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), na.value = "red") + labs(fill = "") + theme(legend.position = c(0.2, 0.3)) ) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/omegas-temp-rw-800-kn.png"), width = 3, height = 3) omega_s_all
plot_map_raster(p_nd, "epsilon_st") + ggtitle("Spatiotemporal random effects only") + facet_wrap(~year) + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), na.value = "red")
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds") p_nd1 <- p_nd %>% filter(ssid!=16) %>% filter(ssid!=4) %>% filter(year>2007) (st_1n3 <- plot_map_raster(p_nd1, "epsilon_st") + ggtitle("Spatiotemporal random effects") + facet_grid(~year) + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), #trans = sqrt, na.value = "red") + labs(fill = "epsilon") + theme( legend.position = "none")) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-1n3.png"), width =6, height = 2) # png( # file = "figs/temp-predictions-1n3.png", # res = 400, # units = "in", # width = 8, # height = 8 # ) # st_1n3 # dev.off()
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds") p_nd4 <- p_nd %>% filter(ssid==4) %>% filter(year>2007)#%>% filter(ssid!=4) (st_4 <- plot_map_raster(p_nd4, "epsilon_st") + #ggtitle(" ") + facet_grid(~year) + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), #trans = sqrt, na.value = "red") + labs(fill = " ") + theme( legend.position = "none") ) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-4.png"), width = 6, height = 2) st_4
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds") p_nd16 <- p_nd %>% filter(ssid==16) %>% filter(year>2007)#%>% filter(ssid!=4) (st_16 <- plot_map_raster(p_nd16, "epsilon_st") + #ggtitle(" ") + facet_grid(~year) + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), na.value = "red") + xlim(20,360) + ylim(5850,6050)+ # xlim(65,345) + ylim(5800,6100)+ theme(panel.border = element_blank(), legend.position = "none") ) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-16.png"), width = 6, height = 1) st_16
png( file = "figs/temp-predictions.png", res = 600, units = "in", width = 11, height = 10 ) gridExtra::grid.arrange( grobs = list(pred_1, st_1n3, omega_s_all, pred_4, st_4, st_16 ), widths = c(2, 2, 1.75), heights = c(2, 1.35), layout_matrix = rbind(c(1, 2, 3), c(4, 5, 6)), top = grid::textGrob("Bottom temperature estimates from sdmTMB model") ) # gridExtra::grid.arrange( nrow = 2, heights = c(2,1.5), # top = grid::textGrob("Dissolved oxygen estimates from sdmTMB model")) 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.