pivot_longer(
cols = starts_with("resid_"),
names_to = "Model",
values_to = "Residual"
) %>%
# Optional: Clean up the model names for the plot labels
mutate(Model = gsub("resid_", "", Model),
Model = factor(Model, levels = c("lm", "lm_krig", "gam", "gam_krig",
"ranger", "ranger_krig", "xgboost",
"xgboost_krig", "krig")))
# 3. Convert to a spatial (sf) object
# Replace 'Longitude' and 'Latitude' with your actual coordinate column names
# crs = 4326 is standard GPS coordinates (WGS84)
val_sf <- st_as_sf(val_long,
coords = c("longitude", "latitude"),
crs = 4326)
val_sf <- val_sf[val_sf$Model %in% c("lm","lm_krig", "ranger","ranger_krig"),]
model_names <- c(
"lm"           = "Linear Regression",
"lm_krig"      = "Linear Regression + Kriging",
"ranger"       = "Random Forest",
"ranger_krig"  = "Random Forest + Kriging"
)
counties <- counties(state = "WI", cb = TRUE, year = 2022) # 'cb = TRUE' for cartographic boundary (simplified) version
counties = counties[,c("NAME","geometry")]
colnames(counties) = c("County","geometry")
ggplot() +
geom_sf(data = val_sf[val_sf$Model %in% c("lm","lm_krig", "ranger","ranger_krig"),],
aes(color = pmin(pmax(Residual,-1),1)), # trim it to make to look better
# aes(color = Residual),
size = 0.01) +
geom_sf(data = counties, fill = NA, color = "grey80", size = 0.2) +
scale_color_gradient2(low = "#000080",
mid = "#ffffbf",
high = "#e31a1c",
name = expression(paste("Prediction Error\n[log10(mg/L)]"))
) +
facet_wrap(~Model, labeller = as_labeller(model_names)) +
theme_minimal() +
theme(legend.position = "right",
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
strip.text = element_text(size = 12))
ggsave("figures/prediction_error_map.png", width = 10, height = 8, dpi = 300)
#### transform back to the original scale
rmse_data <- data.frame(
Method = c("LR", "LR", "GAM", "GAM", "XGBoost", "XGBoost", "RF", "RF", "Kriging"),
Variant = c("Base only", "Base and Kriging", "Base only", "Base and Kriging",
"Base only", "Base and Kriging", "Base only", "Base and Kriging", "Kriging Only"),
RMSE = c(
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_lm)),
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_lm_krig)),
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_gam)),
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_gam_krig)),
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_xgboost)),
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_xgboost_krig)),
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_ranger)),
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_ranger_krig)),
Metrics::rmse(back_trans(data_split$data_test$logconcentration_plus_median), back_trans(data_split$data_test$pred_krig))
)
)
#### 2. Create the Barplot
ggplot(rmse_data, aes(x = reorder(Method, -RMSE), y = RMSE, fill = Variant)) +
geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7) +
geom_text(aes(label = round(RMSE, 3)),
position = position_dodge(width = 0.8),
vjust = -0.5, size = 4.5) +
scale_fill_manual(values = c("Base only" = "#A6CEE3",
"Kriging Only" = "#B2DF8A",
"Base and Kriging" = "#1F78B4"
)) +
coord_cartesian(ylim = c(2.5,4.4))+
theme_minimal() +
labs(title = "Validation RMSE Comparison",
x = "Model Type",
y = "Root Mean Square Error (mg/L)",
fill = "Method") +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.title = element_text(size = 14, face = "bold"),
axis.text = element_text(size = 12, color = "black"),
legend.title = element_text(size = 13, face = "bold"),
legend.text = element_text(size = 12),
legend.position = c(0.98, 0.98),
legend.justification = c("right", "top"),
legend.background = element_rect(fill = "white", color = "grey80"))
ggsave("figures/rmse_barplot.png", width = 8, height = 6, dpi = 300)
predictions_3D_grid <- readRDS('results/predictions_3D_grid.rds')
plss_covariates <- read.csv('data/plss_covariates.csv')
pred_matrix <- predictions_3D_grid$pred_matrix
krige_values_PLSS_ranger <- predictions_3D_grid$krige_values_PLSS_ranger
treatment_values <- predictions_3D_grid$log_well_depth_grid
depths_meters <- c(15,30,60,120)
feet_meters_scale <- 0.3048
depths_feet <- depths_meters/feet_meters_scale
grid_indice <- lapply(depths_feet, FUN = function(x){which.min(abs(treatment_values - feet))}, simplify = TRUE)
all_depths_df <- map_dfr(depths_feet, function(d) {
plss_covariates %>%
mutate(
WellDepth = d,
logWellDepth = log(d),
Depth_Label = paste(round(d*0.3048,0), "Meters"), # Useful for facet titles
pred_ranger = pred_matrix[which.min(abs(exp(treatment_values) - d)),],
krige_ranger = krige_values_PLSS_ranger
)
})
# 2. Transform to original scale and apply the floor (pmax)
all_depths_df <- all_depths_df %>%
mutate(
pred_original = exp(pred_ranger + krige_ranger),
pred_original = pmax(pred_original, 0.4274596)
)
# 3. Convert to SF object once
all_depths_sf <- st_as_sf(all_depths_df, coords = c("longitude", "latitude"), crs = 4326)
counties <- counties(state = "WI", cb = TRUE, year = 2022) # 'cb = TRUE' for cartographic boundary (simplified) version
counties = counties[,c("NAME","geometry")]
colnames(counties) = c("County","geometry")
ggplot() +
# Points Layer
geom_sf(data = all_depths_sf, aes(color =pred_original), size = 0.01) +
# Counties Layer (drawn on top for clarity)
geom_sf(data = counties, fill = NA, color = "black", size = 0.1) +
# Shared Color Scale
scale_color_viridis_c(
option = "plasma",
name = "Predicted Nitrate\n(mg/L)",
trans = "log10",
# trans = "log",
limits = c(min(all_depths_df$pred_original), max(all_depths_df$pred_original)),
breaks = c(1, 5, 10, 20, 40),
labels = c( "1", "5", "10", "20", "40")
) +
# Create the 2x2 Grid
facet_wrap(~reorder(Depth_Label, WellDepth), ncol = 2) +
labs(
title = "Nitrate Concentration Predictions by Depth",
subtitle = "Integration of Residual Kriging and Random Forest",
# x = "Longitude",
# y = "Latitude"
) +
theme_minimal() +
theme(
strip.text = element_text(size = 12, face = "plain"),
legend.position = "right",
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.key.width = unit(0.7, "cm")
)
ggsave("figures/3D_map.png", width = 10, height = 8, dpi = 300)
ggplot() +
# Points Layer
geom_sf(data = all_depths_sf, aes(color =pred_original), size = 0.01) +
# Counties Layer (drawn on top for clarity)
geom_sf(data = counties, fill = NA, color = "black", size = 0.1) +
# Shared Color Scale
scale_color_viridis_c(
option = "plasma",
name = "Predicted Nitrate\n(mg/L)",
trans = "log10",
# trans = "log",
limits = c(min(all_depths_df$pred_original), max(all_depths_df$pred_original)),
breaks = c(1, 5, 10, 20, 40),
labels = c( "1", "5", "10", "20", "40")
) +
# Create the 2x2 Grid
facet_wrap(~reorder(Depth_Label, WellDepth), ncol = 2) +
labs(
title = "Nitrate Concentration Predictions by Depth",
# subtitle = "Integration of Residual Kriging and Random Forest",
# x = "Longitude",
# y = "Latitude"
) +
theme_minimal() +
theme(
strip.text = element_text(size = 12, face = "plain"),
legend.position = "right",
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.key.width = unit(0.7, "cm")
)
ggsave("figures/3D_map.png", width = 10, height = 8, dpi = 300)
pred_matrix <- predictions_3D_grid$pred_matrix
krige_values_PLSS_ranger <- predictions_3D_grid$krige_values_PLSS_ranger
treatment_values <- predictions_3D_grid$log_well_depth_grid
plss_covariates <- read.csv(file = "/data/plss_covariates.csv")
predictions_3D_grid <- readRDS('results/predictions_3D_grid.rds')
pred_matrix <- predictions_3D_grid$pred_matrix
krige_values_PLSS_ranger <- predictions_3D_grid$krige_values_PLSS_ranger
treatment_values <- predictions_3D_grid$log_well_depth_grid
plss_covariates <- read.csv(file = "/data/plss_covariates.csv")
plss_covariates <- read.csv(file = "data/plss_covariates.csv")
adjusted_matrix <- sweep(pred_matrix, 2, krige_values_PLSS_ranger, "+")
pass_condition_2mg <- adjusted_matrix < log(2.5)
first_indices_2mg <- max.col(t(pass_condition_2mg), ties.method = "first")
has_any_pass_2mg <- colSums(pass_condition_2mg) > 0
first_indices_2mg[!has_any_pass_2mg] <- nrow(pred_matrix)
pass_condition_10mg <- adjusted_matrix < log(10.5)
first_indices_10mg <- max.col(t(pass_condition_10mg), ties.method = "first")
has_any_pass_10mg <- colSums(pass_condition_10mg) > 0
first_indices_10mg[!has_any_pass_10mg] <- nrow(pred_matrix)
indirect_policy_ranger_10mg <- treatment_values[first_indices_10mg]
indirect_policy_ranger_2mg <- treatment_values[first_indices_2mg]
plss_covariates$indirect_policy_ranger_10mg <- indirect_policy_ranger_10mg
plss_covariates$indirect_policy_ranger_2mg <- indirect_policy_ranger_2mg
counties <- counties(state = "WI", cb = TRUE, year = 2022) # 'cb = TRUE' for cartographic boundary (simplified) version
counties = counties[,c("NAME","geometry")]
colnames(counties) = c("County","geometry")
plss_long_sf <- plss_covariates %>%
dplyr::select(longitude, latitude, indirect_policy_ranger_2mg, indirect_policy_ranger_10mg) %>%
pivot_longer(cols = starts_with("indirect_policy"),
names_to = "Threshold",
values_to = "Depth") %>%
mutate(
# Apply the cap and rename for the plot titles
Threshold = recode(Threshold,
"indirect_policy_ranger_2mg" = "2mg/L Threshold",
"indirect_policy_ranger_10mg" = "10mg/L Threshold")
) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
meters <- c(5, 25, 100, 400)
ggplot() +
# 1. Background layer
geom_sf(data = counties, fill = NA, color = "grey80", size = 0.1) +
# 2. Prediction points
geom_sf(data = plss_long_sf, aes(color = Depth), size = 0.05, alpha = 0.6) +
# 3. Shared Color Scale
# we will choose 5m, 50m, 100m ,300m
scale_color_viridis_c(
option = "turbo", # "H" in newer ggplot is turbo
name = "Required Well Depth\n(Meter)",
trans = "log",
limits = c(log(10), log(1350)),
# breaks = log(c(10, 40, 200, 1000)),
breaks = log(meters/0.3048),
labels = function(x) round(exp(x)*0.3048, 0)
) +
# 4. Faceting into 1 row, 2 columns
facet_wrap(~Threshold, ncol = 2) +
# 5. Styling
labs(
title = "Minimum Safe Well Depth for Nitrate Thresholds",
x = NULL, y = NULL
) +
theme_minimal() +
theme(
strip.text = element_text(size = 14, face = "plain"), # Facet Titles
legend.position = "right",
legend.key.width = unit(0.5, "cm"),
panel.spacing = unit(2, "lines") # Adds space between the two maps
)
ggsave("figures/required_depth_map.png", width = 12, height = 6, dpi = 300)
ggsave("figures/required_depth_map.png", width = 8, height = 4, dpi = 300)
ggsave("figures/required_depth_map.png", width = 10, height = 5, dpi = 300)
dim(data_split$data_test)
data_less_two <- data[data$concentration_plus_median >= 2.5, ]
nitrate_fitted_ranger_less_two <- nitrate_prediction_ranger(
data = data_split_less_two$data,
data_test = data_split_less_two$data_test
)
data_less_two <- data[data$concentration_plus_median >= 2.5, ]
data_split_less_two <- split_nitrate_data(data_less_two)
data
dim(data)
data <- load_nitrate_data(file_path = "/data/data_Nitrate_with_covar.csv", zero_inflated = TRUE)
data <- load_nitrate_data(file_path = "data/data_Nitrate_with_covar.csv", zero_inflated = TRUE)
data_split_less_two <- split_nitrate_data(data_less_two)
data <- load_nitrate_data(file_path = "data/data_Nitrate_with_covar.csv", zero_inflated = TRUE)
sum(data$concentration_plus_median >= 2.5)
sum(data$concentration_plus_median >= 2)
sum(data$concentration_plus_median >= 0.5)
data_less_two <- data[data$concentration_plus_median >= 2.5, ]
data_split_less_two <- split_nitrate_data(data_less_two)
nitrate_fitted_ranger_less_two <- nitrate_prediction_ranger(
data = data_split_less_two$data,
data_test = data_split_less_two$data_test
)
data_split_less_two$data_test$pred_ranger      <- nitrate_fitted_ranger_less_two$pred_test
data_split_less_two$data_test$pred_ranger_krig <- nitrate_fitted_ranger_less_two$pred_test + nitrate_fitted_ranger_less_two$krige_values_test
### linear regression (trained on original scale)
nitrate_fitted_lm_less_two <- nitrate_prediction_SL(
data = data_split_less_two$data,
data_test = data_split_less_two$data_test,
SL.library = "SL.lm",
response_name = "concentration_plus_median"
)
data_split_less_two$data_test$pred_lm      <- nitrate_fitted_lm_less_two$pred_test
data_split_less_two$data_test$pred_lm_krig <- nitrate_fitted_lm_less_two$pred_test + nitrate_fitted_lm_less_two$krige_values_test
saveRDS(data_split_less_two, "results/data_split_less_two.rds")
data_split_less_two <- readRDS("results/data_split_less_two.rds")
test_data_less_two <- data_split_less_two$data_test
pred_names <- c("pred_ranger_krig", "pred_ranger", "pred_lm", "pred_lm_krig")
display_names <- c("RF + Kriging", "RF", "LR", "LR + Kriging")
name_map <- setNames(display_names, pred_names)
threshold_lt2 <- 10
all_metrics_lt2 <- data.frame()
roc_list_lt2   <- list()
conf_mat_list_lt2 <- list()
all_cm_data_lt2 <- data.frame()
for (m in pred_names) {
y_obs  <- factor(ifelse(test_data_less_two$concentration_plus_median > threshold_lt2, "High", "Low"),
levels = c("High", "Low"))
y_pred <- factor(ifelse(test_data_less_two[[m]] > threshold_lt2, "High", "Low"),
levels = c("High", "Low"))
conf_mat <- confusionMatrix(y_pred, y_obs, positive = "High")
conf_mat_list_lt2[[name_map[m]]] <- conf_mat
df <- as.data.frame(conf_mat$table)
df$Method <- name_map[m]
all_cm_data_lt2 <- rbind(all_cm_data_lt2, df)
eval_df <- data.frame(truth = y_obs, estimate = y_pred)
mcc_val <- mcc(eval_df, truth = truth, estimate = estimate)$.estimate
f1_val  <- F1_Score(y_obs, y_pred)
roc_obj <- roc(test_data_less_two$concentration_plus_median > threshold_lt2,
test_data_less_two[[m]], quiet = TRUE)
auc_val <- as.numeric(auc(roc_obj))
roc_list_lt2[[m]] <- roc_obj
model_stats <- data.frame(
Method = name_map[m],
Accuracy    = conf_mat$overall["Accuracy"],
Sensitivity = conf_mat$byClass["Sensitivity"],
Specificity = conf_mat$byClass["Specificity"],
MCC = mcc_val,
F1  = f1_val,
AUC = auc_val
)
all_metrics_lt2 <- rbind(all_metrics_lt2, model_stats)
}
plot_data_lt2 <- all_metrics_lt2 %>%
pivot_longer(cols = -Method, names_to = "Metric", values_to = "Score")
ggplot(plot_data_lt2, aes(x = Method, y = Score, fill = Method)) +
geom_col(alpha = 0.8) +
facet_wrap(~Metric, scales = "free_y") +
theme_minimal() +
labs(title = "Comparative Model Performance of Binary Classification",
subtitle = paste("Threshold:", threshold_lt2, "mg/L | Trained without <2mg/L observations")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
scale_fill_brewer(palette = "Set1")
for (m in pred_names) {
y_obs  <- factor(ifelse(test_data_less_two$concentration_plus_median > threshold_lt2, "High", "Low"),
levels = c("High", "Low"))
y_pred <- factor(ifelse(back_trans(test_data_less_two[[m]]) > threshold_lt2, "High", "Low"),
levels = c("High", "Low"))
conf_mat <- confusionMatrix(y_pred, y_obs, positive = "High")
conf_mat_list_lt2[[name_map[m]]] <- conf_mat
df <- as.data.frame(conf_mat$table)
df$Method <- name_map[m]
all_cm_data_lt2 <- rbind(all_cm_data_lt2, df)
eval_df <- data.frame(truth = y_obs, estimate = y_pred)
mcc_val <- mcc(eval_df, truth = truth, estimate = estimate)$.estimate
f1_val  <- F1_Score(y_obs, y_pred)
roc_obj <- roc(test_data_less_two$concentration_plus_median > threshold_lt2,
test_data_less_two[[m]], quiet = TRUE)
auc_val <- as.numeric(auc(roc_obj))
roc_list_lt2[[m]] <- roc_obj
model_stats <- data.frame(
Method = name_map[m],
Accuracy    = conf_mat$overall["Accuracy"],
Sensitivity = conf_mat$byClass["Sensitivity"],
Specificity = conf_mat$byClass["Specificity"],
MCC = mcc_val,
F1  = f1_val,
AUC = auc_val
)
all_metrics_lt2 <- rbind(all_metrics_lt2, model_stats)
}
plot_data_lt2 <- all_metrics_lt2 %>%
pivot_longer(cols = -Method, names_to = "Metric", values_to = "Score")
ggplot(plot_data_lt2, aes(x = Method, y = Score, fill = Method)) +
geom_col(alpha = 0.8) +
facet_wrap(~Metric, scales = "free_y") +
theme_minimal() +
labs(title = "Comparative Model Performance of Binary Classification",
subtitle = paste("Threshold:", threshold_lt2, "mg/L | Trained without <2mg/L observations")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
scale_fill_brewer(palette = "Set1")
### linear regression (trained on original scale)
nitrate_fitted_lm_less_two <- nitrate_prediction_SL(
data = data_split_less_two$data,
data_test = data_split_less_two$data_test,
SL.library = "SL.lm"
)
data_split_less_two$data_test$pred_lm      <- nitrate_fitted_lm_less_two$pred_test
data_split_less_two$data_test$pred_lm_krig <- nitrate_fitted_lm_less_two$pred_test + nitrate_fitted_lm_less_two$krige_values_test
saveRDS(data_split_less_two, "results/data_split_less_two.rds")
summary(data_split_less_two$data_test$pred_lm_krig)
summary(exp(data_split_less_two$data_test$pred_lm_krig))
summary(exp(data_split_less_two$data_test$pred_lm))
summary(exp(data_split_less_two$data_test$pred_ranger))
summary(exp(data_split_less_two$data_test$pred_ranger_krig))
data_split_less_two <- readRDS("results/data_split_less_two.rds")
test_data_less_two <- data_split_less_two$data_test
pred_names <- c("pred_ranger_krig", "pred_ranger", "pred_lm", "pred_lm_krig")
display_names <- c("RF + Kriging", "RF", "LR", "LR + Kriging")
name_map <- setNames(display_names, pred_names)
threshold_lt2 <- 10
all_metrics_lt2 <- data.frame()
roc_list_lt2   <- list()
conf_mat_list_lt2 <- list()
all_cm_data_lt2 <- data.frame()
for (m in pred_names) {
y_obs  <- factor(ifelse(test_data_less_two$concentration_plus_median > threshold_lt2, "High", "Low"),
levels = c("High", "Low"))
y_pred <- factor(ifelse(back_trans(test_data_less_two[[m]]) > threshold_lt2, "High", "Low"),
levels = c("High", "Low"))
conf_mat <- confusionMatrix(y_pred, y_obs, positive = "High")
conf_mat_list_lt2[[name_map[m]]] <- conf_mat
df <- as.data.frame(conf_mat$table)
df$Method <- name_map[m]
all_cm_data_lt2 <- rbind(all_cm_data_lt2, df)
eval_df <- data.frame(truth = y_obs, estimate = y_pred)
mcc_val <- mcc(eval_df, truth = truth, estimate = estimate)$.estimate
f1_val  <- F1_Score(y_obs, y_pred)
roc_obj <- roc(test_data_less_two$concentration_plus_median > threshold_lt2,
test_data_less_two[[m]], quiet = TRUE)
auc_val <- as.numeric(auc(roc_obj))
roc_list_lt2[[m]] <- roc_obj
model_stats <- data.frame(
Method = name_map[m],
Accuracy    = conf_mat$overall["Accuracy"],
Sensitivity = conf_mat$byClass["Sensitivity"],
Specificity = conf_mat$byClass["Specificity"],
MCC = mcc_val,
F1  = f1_val,
AUC = auc_val
)
all_metrics_lt2 <- rbind(all_metrics_lt2, model_stats)
}
plot_data_lt2 <- all_metrics_lt2 %>%
pivot_longer(cols = -Method, names_to = "Metric", values_to = "Score")
ggplot(plot_data_lt2, aes(x = Method, y = Score, fill = Method)) +
geom_col(alpha = 0.8) +
facet_wrap(~Metric, scales = "free_y") +
theme_minimal() +
labs(title = "Comparative Model Performance of Binary Classification",
subtitle = paste("Threshold:", threshold_lt2, "mg/L | Trained without <2mg/L observations")) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
scale_fill_brewer(palette = "Set1")
library(plotly)
library(htmlwidgets)
depth_labels <- c("15 m", "30 m", "60 m", "120 m")
depth_vals   <- c(15, 30, 60, 120)
# Add depth in metres and log10 nitrate to the flat data frame
all_depths_3d <- all_depths_df %>%
mutate(
depth_m      = round(WellDepth * feet_meters_scale, 0),
log10_nitrate = log10(pred_original)
)
# Build one trace per depth so each slice gets its own legend entry
fig_3d <- plot_ly()
trace_colours <- c("#0d0887", "#7e03a8", "#cc4778", "#f89540")  # plasma-ish anchors
for (i in seq_along(depth_vals)) {
slice <- all_depths_3d %>% filter(depth_m == depth_vals[i])
fig_3d <- fig_3d %>%
add_trace(
data = slice,
type = "scatter3d", mode = "markers",
x = ~longitude, y = ~latitude, z = ~depth_m,
marker = list(
size    = 2,
opacity = 0.75,
color   = ~log10_nitrate,
colorscale = "Plasma",
showscale  = (i == 1),           # show colour bar only once
cmin = log10(min(all_depths_df$pred_original)),
cmax = log10(max(all_depths_df$pred_original)),
colorbar = list(
title      = "Nitrate (mg/L)<br>log\u2081\u2080 scale",
tickvals   = log10(c(1, 5, 10, 20, 40)),
ticktext   = c("1", "5", "10", "20", "40"),
thickness  = 15,
len        = 0.6
)
),
text      = ~paste0("Depth: ", depth_m, " m<br>",
"Nitrate: ", round(pred_original, 2), " mg/L<br>",
"Lon: ", round(longitude, 3), "<br>",
"Lat: ", round(latitude, 3)),
hoverinfo = "text",
name      = depth_labels[i],
showlegend = TRUE
)
}
fig_3d <- fig_3d %>%
layout(
title = list(text = "Nitrate Concentration Predictions by Depth",
font = list(size = 16)),
scene = list(
xaxis = list(title = "Longitude"),
yaxis = list(title = "Latitude"),
zaxis = list(
title    = "Depth (m)",
tickvals = depth_vals,
ticktext = depth_labels,
autorange = "reversed",          # 15 m on top, 120 m at bottom
showgrid = TRUE,
gridcolor = "lightgrey"
),
camera = list(
eye = list(x = 1.6, y = -1.8, z = 1.2)   # default viewing angle
),
aspectmode = "manual",
aspectratio = list(x = 1.2, y = 1.5, z = 0.8)
),
legend = list(title = list(text = "Depth"), x = 0.02, y = 0.95)
)
fig_3d
# Save as interactive HTML
saveWidget(fig_3d, file = "figures/3D_nitrate_map.html", selfcontained = TRUE)
# Save as static PNG (requires the webshot2 package + Chrome/Chromium)
if (requireNamespace("webshot2", quietly = TRUE)) {
webshot2::webshot("figures/3D_nitrate_map.html",
file   = "figures/3D_nitrate_map.png",
vwidth = 1200, vheight = 900, delay = 2)
}
