Reworked Biodiversity Prediction
Alan Ivory
2026-08-10
Load and Prepare Data
model_vars <- attr(rf_final$terms, "term.labels")
cat("Model expects these predictors:\n")## Model expects these predictors:print(model_vars)## [1] "Cogongrass_Cover" "elevation" "slope" "pr"
## [5] "srad" "TreeCanopy_NLCD" "SoilType" "ndvi"
## [9] "avg_pop_den_1km"env_vars <- model_vars
cat(sprintf("\nEnvironmental predictors : %s\n", paste(env_vars, collapse = ", ")))##
## Environmental predictors : Cogongrass_Cover, elevation, slope, pr, srad, TreeCanopy_NLCD, SoilType, ndvi, avg_pop_den_1kmNLCD
env_vars <- env_vars[env_vars != "Cogongrass_Cover"]
raster_paths <- list(
pr = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/5_Year_Average_Precipitation.tif",
sph = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/5_Year_Average_Specific_Humidity.tif",
srad = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/5_Year_Average_Radiation.tif",
TreeCanopy_NLCD = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/Tree_Cover_SE.tif",
ndvi = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/CombinedMedianNDVI.tif",
tmmn = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/5_Year_Average_Min_Temperature.tif",
elevation = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/Elevation_CONUS.tif",
SoilType = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/taxorder_CONUS_30m.tif",
avg_pop_den_1km = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/US_PopDensity_2020.tif",
nlcd = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/Annual_NLCD_LndCov_2024_CU_C1V1.tif",
slope = "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/Slope_30m_Merged.tif"
)
missing_rasts <- setdiff(env_vars, names(raster_paths))
if (length(missing_rasts) > 0) {
stop(sprintf("Missing raster paths for: %s", paste(missing_rasts, collapse = ", ")))
}
## Load only the rasters needed by the model
cat("Loading rasters...\n")## Loading rasters...rast_list <- lapply(env_vars, function(v) {
cat(sprintf(" Loading: %s\n", v))
rast(raster_paths[[v]])
})## Loading: elevation
## Loading: slope
## Loading: pr
## Loading: srad
## Loading: TreeCanopy_NLCD
## Loading: SoilType
## Loading: ndvi
## Loading: avg_pop_den_1kmnames(rast_list) <- env_vars
# Loop through and plot each raster
for (name in names(raster_paths)) {
r <- rast(raster_paths[[name]])
plot(r, main = name)
}
Align Rasters
cat("CRS of each raster:\n")## CRS of each raster:for (v in names(rast_list)) {
cat(sprintf(" %-20s %s\n", v, crs(rast_list[[v]], describe = TRUE)$code))
}## elevation 4326
## slope 5070
## pr 4326
## srad 4326
## TreeCanopy_NLCD 5070
## SoilType 5070
## ndvi 4326
## avg_pop_den_1km 4326# Use NDVI raster extent as the boundary
cat("\nExtracting extent from NDVI raster...\n")##
## Extracting extent from NDVI raster...ndvi_rast <- rast(raster_paths[["ndvi"]])
ndvi_crs <- crs(ndvi_rast)
ndvi_extent <- ext(ndvi_rast)
cat(sprintf("NDVI CRS: %s\n", crs(ndvi_rast, describe = TRUE)$code))## NDVI CRS: 4326cat("NDVI Extent:\n")## NDVI Extent:print(ndvi_extent)## SpatExtent : -131.325070711546, -59.42391537062, 19.9520316179366, 55.2019233667867 (xmin, xmax, ymin, ymax)# Load, crop each raster to NDVI extent
local_dir <- "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/local_cache"
dir.create(local_dir, showWarnings = FALSE, recursive = TRUE)
rast_list <- lapply(names(raster_paths), function(v) {
src_path <- raster_paths[[v]]
local_path <- file.path(local_dir, paste0(v, ".tif"))
if (!file.exists(local_path)) {
cat(sprintf(" Copying %-20s to local cache...\n", v))
r_src <- rast(src_path)
writeRaster(r_src, local_path, overwrite = TRUE)
}
r <- rast(local_path)
# Reproject the NDVI extent into this raster's CRS for cropping
if (!same.crs(r, ndvi_rast)) {
ndvi_ext_reproj <- ext(project(vect(ndvi_extent, crs = ndvi_crs), crs(r)))
} else {
ndvi_ext_reproj <- ndvi_extent
}
r_crop <- crop(r, ndvi_ext_reproj)
cat(sprintf(" %-20s done | dims: %d x %d | CRS: %s\n",
v, nrow(r_crop), ncol(r_crop),
crs(r_crop, describe = TRUE)$code))
r_crop
})## pr done | dims: 3924 x 8004 | CRS: 4326
## sph done | dims: 3924 x 8004 | CRS: 4326
## srad done | dims: 3924 x 8004 | CRS: 4326
## |---------|---------|---------|---------|========================================= TreeCanopy_NLCD done | dims: 93792 x 154328 | CRS: 5070
## ndvi done | dims: 3924 x 8004 | CRS: 4326
## tmmn done | dims: 3924 x 8004 | CRS: 4326
## elevation done | dims: 1252 x 2841 | CRS: 4326
## |---------|---------|---------|---------|========================================= SoilType done | dims: 93376 x 153996 | CRS: 5070
## avg_pop_den_1km done | dims: 2766 x 6433 | CRS: 4326
## |---------|---------|---------|---------|========================================= nlcd done | dims: 98117 x 160000 | CRS: NA
## |---------|---------|---------|---------|========================================= slope done | dims: 54490 x 109186 | CRS: 5070names(rast_list) <- names(raster_paths)Env Stack for Prediction
# Use the NDVI raster as the reference grid for alignment
ref_rast <- rast_list[["ndvi"]]
cat("Aligning all rasters to the NDVI 30m reference grid...\n")## Aligning all rasters to the NDVI 30m reference grid...cat("\nReprojecting and resampling TreeCanopy and SoilType to match ref grid...\n")##
## Reprojecting and resampling TreeCanopy and SoilType to match ref grid...rast_aligned <- lapply(names(rast_list), function(v) {
r <- rast_list[[v]]
# Define categorical layers that MUST use Nearest Neighbor
categorical_layers <- c("SoilType", "nlcd")
method <- if (v %in% categorical_layers) "near" else "bilinear"
if (!same.crs(r, ref_rast)) {
cat(sprintf(" Reprojecting: %s to 30m...\n", v))
r <- project(r, ref_rast, method = method)
} else if (!compareGeom(r, ref_rast, stopOnError = FALSE)) {
cat(sprintf(" Resampling: %s to 30m...\n", v))
r <- resample(r, ref_rast, method = method)
} else {
cat(sprintf(" Already Aligned: %s\n", v))
}
r
})## Already Aligned: pr
## Already Aligned: sph
## Already Aligned: srad
## Reprojecting: TreeCanopy_NLCD to 30m...
## |---------|---------|---------|---------|========================================= Already Aligned: ndvi
## Already Aligned: tmmn
## Resampling: elevation to 30m...
## Reprojecting: SoilType to 30m...
## |---------|---------|---------|---------|========================================= Resampling: avg_pop_den_1km to 30m...
## Reprojecting: nlcd to 30m...
## |---------|---------|---------|---------|========================================= Reprojecting: slope to 30m...
## |---------|---------|---------|---------|========================================= names(rast_aligned) <- names(rast_list)
env_stack <- rast(rast_aligned)
names(env_stack) <- names(rast_aligned)
cat(sprintf("\nenv_stack: %d layers, %d x %d cells\n",
nlyr(env_stack), nrow(env_stack), ncol(env_stack)))##
## env_stack: 11 layers, 3924 x 8004 cells# Check for NaN/NA in each layer
cat("\nNon-NA cell counts per layer:\n")##
## Non-NA cell counts per layer:for (v in names(env_stack)) {
n <- global(env_stack[[v]], fun = "notNA")[1,1]
cat(sprintf(" %-20s %d\n", v, n))
}## pr 10408368
## sph 10408368
## srad 10408368
## TreeCanopy_NLCD 16872041
## ndvi 6217899
## tmmn 10408368
## elevation 18998331
## SoilType 9512108
## avg_pop_den_1km 10160873
## nlcd 10514534
## slope 4124352NLCD Masking
Only focus on forest, shrubland, or heraceous areas NLCD Classes: 41, 42, 43 (Forests); 52 (Shrub); 71 (Herbaceous); 90 (Woody Wetland); 95 (Emergent Wetland)
mask_habitats <- function(env_stack, nlcd_raster) {
cat("Masking landscape to Forest, Shrubland, and Herbaceous cover...\n")
# Align NLCD to the stack
if(!compareGeom(nlcd_raster, env_stack, stopOnError = FALSE)) {
nlcd_raster <- project(nlcd_raster, env_stack, method = "near")
}
# Target classes
target_classes <- c(41, 42, 43, 52, 71, 90, 95)
mask_layer <- nlcd_raster %in% target_classes
mask_layer[mask_layer == 0] <- NA
# Apply mask to the entire stack
env_stack_masked <- mask(env_stack, mask_layer)
return(env_stack_masked)
}
env_stack <- mask_habitats(env_stack, rast_list[["nlcd"]])## Masking landscape to Forest, Shrubland, and Herbaceous cover...
## |---------|---------|---------|---------|========================================= Longleaf Pine Mask
# Path to the longleaf pine range shapefile.
# NOTE: st_read needs the companion files (.shx, .dbf, .prj) in the same folder.
longleaf_path <- "C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/pinupalu.shp"
cat("Loading longleaf pine historic range...\n")## Loading longleaf pine historic range...longleaf_range <- sf::st_read(longleaf_path, quiet = TRUE)
# The Little's range shapefiles often ship without a .prj, so the CRS comes in
# undefined. They are in geographic coordinates (lon/lat, decimal degrees), so we
# assign that here before reprojecting. Change assumed_crs if you know the file's
# true CRS (the documented original datum for these maps is NAD27 = EPSG:4267).
assumed_crs <- 4326 # WGS84 lon/lat
if (is.na(sf::st_crs(longleaf_range))) {
cat(sprintf(" Shapefile has no CRS; assigning EPSG:%d (assumed).\n", assumed_crs))
sf::st_crs(longleaf_range) <- assumed_crs
}## Shapefile has no CRS; assigning EPSG:4326 (assumed).# Match the CRS of the environmental stack, then convert to a terra vector
longleaf_range <- sf::st_transform(longleaf_range, crs(env_stack))
longleaf_vect <- terra::vect(longleaf_range)
# Crop to the range's bounding box (large speed-up on the 101 scenarios),
# then mask so that cells outside the range polygon(s) become NA.
cat("Cropping and masking env_stack to the longleaf pine range...\n")## Cropping and masking env_stack to the longleaf pine range...env_stack <- terra::crop(env_stack, longleaf_vect)
env_stack <- terra::mask(env_stack, longleaf_vect)
cat(sprintf("env_stack restricted to longleaf range: %d x %d cells\n",
nrow(env_stack), ncol(env_stack)))## env_stack restricted to longleaf range: 1139 x 2162 cells# Quick visual check of the masked extent
plot(env_stack[[1]], main = "env_stack masked to longleaf pine range")
plot(sf::st_geometry(longleaf_range), add = TRUE, border = "black")
Predict Function
# Builds predictor dataframe, inserts cogongrass scenario, predicts, returns raster
predict_diversity <- function(cogongrass_value, env_stack, rf_model) {
cat(sprintf(" Predicting at Cogongrass_Cover = %d%%...\n",
cogongrass_value))
# Convert stack to data frame
pred_df <- as.data.frame(env_stack, xy = FALSE, na.rm = FALSE)
# Only predict on valid cells
valid_cells <- complete.cases(pred_df)
if(sum(valid_cells) == 0)
stop("No valid cells found for prediction.")
# Add scenario variable
pred_df$Cogongrass_Cover <- cogongrass_value
# Match model variables
all_vars <- attr(rf_model$terms, "term.labels")
pred_df_valid <- pred_df[valid_cells,
all_vars, drop = FALSE]
# Predict
predicted_vals <- predict(rf_model,
newdata = pred_df_valid)
# Map back to raster
predicted_full <- rep(NA_real_,
nrow(pred_df))
predicted_full[valid_cells] <- predicted_vals
result_rast <- env_stack[[1]]
values(result_rast) <- predicted_full
names(result_rast) <- paste0("Shannon_",
cogongrass_value)
return(result_rast)
}Run Cogongrass Scenarios
scenario_levels <- seq(0, 100, by = 1)
prediction_list <- list()
cat("\n========== Running All Cogongrass Scenarios ==========\n") ##
## ========== Running All Cogongrass Scenarios ==========for(level in scenario_levels) {
prediction_list[[paste0("p", level)]] <- predict_diversity(level,
env_stack, rf_final)
} ## Predicting at Cogongrass_Cover = 0%...
## Predicting at Cogongrass_Cover = 1%...
## Predicting at Cogongrass_Cover = 2%...
## Predicting at Cogongrass_Cover = 3%...
## Predicting at Cogongrass_Cover = 4%...
## Predicting at Cogongrass_Cover = 5%...
## Predicting at Cogongrass_Cover = 6%...
## Predicting at Cogongrass_Cover = 7%...
## Predicting at Cogongrass_Cover = 8%...
## Predicting at Cogongrass_Cover = 9%...
## Predicting at Cogongrass_Cover = 10%...
## Predicting at Cogongrass_Cover = 11%...
## Predicting at Cogongrass_Cover = 12%...
## Predicting at Cogongrass_Cover = 13%...
## Predicting at Cogongrass_Cover = 14%...
## Predicting at Cogongrass_Cover = 15%...
## Predicting at Cogongrass_Cover = 16%...
## Predicting at Cogongrass_Cover = 17%...
## Predicting at Cogongrass_Cover = 18%...
## Predicting at Cogongrass_Cover = 19%...
## Predicting at Cogongrass_Cover = 20%...
## Predicting at Cogongrass_Cover = 21%...
## Predicting at Cogongrass_Cover = 22%...
## Predicting at Cogongrass_Cover = 23%...
## Predicting at Cogongrass_Cover = 24%...
## Predicting at Cogongrass_Cover = 25%...
## Predicting at Cogongrass_Cover = 26%...
## Predicting at Cogongrass_Cover = 27%...
## Predicting at Cogongrass_Cover = 28%...
## Predicting at Cogongrass_Cover = 29%...
## Predicting at Cogongrass_Cover = 30%...
## Predicting at Cogongrass_Cover = 31%...
## Predicting at Cogongrass_Cover = 32%...
## Predicting at Cogongrass_Cover = 33%...
## Predicting at Cogongrass_Cover = 34%...
## Predicting at Cogongrass_Cover = 35%...
## Predicting at Cogongrass_Cover = 36%...
## Predicting at Cogongrass_Cover = 37%...
## Predicting at Cogongrass_Cover = 38%...
## Predicting at Cogongrass_Cover = 39%...
## Predicting at Cogongrass_Cover = 40%...
## Predicting at Cogongrass_Cover = 41%...
## Predicting at Cogongrass_Cover = 42%...
## Predicting at Cogongrass_Cover = 43%...
## Predicting at Cogongrass_Cover = 44%...
## Predicting at Cogongrass_Cover = 45%...
## Predicting at Cogongrass_Cover = 46%...
## Predicting at Cogongrass_Cover = 47%...
## Predicting at Cogongrass_Cover = 48%...
## Predicting at Cogongrass_Cover = 49%...
## Predicting at Cogongrass_Cover = 50%...
## Predicting at Cogongrass_Cover = 51%...
## Predicting at Cogongrass_Cover = 52%...
## Predicting at Cogongrass_Cover = 53%...
## Predicting at Cogongrass_Cover = 54%...
## Predicting at Cogongrass_Cover = 55%...
## Predicting at Cogongrass_Cover = 56%...
## Predicting at Cogongrass_Cover = 57%...
## Predicting at Cogongrass_Cover = 58%...
## Predicting at Cogongrass_Cover = 59%...
## Predicting at Cogongrass_Cover = 60%...
## Predicting at Cogongrass_Cover = 61%...
## Predicting at Cogongrass_Cover = 62%...
## Predicting at Cogongrass_Cover = 63%...
## Predicting at Cogongrass_Cover = 64%...
## Predicting at Cogongrass_Cover = 65%...
## Predicting at Cogongrass_Cover = 66%...
## Predicting at Cogongrass_Cover = 67%...
## Predicting at Cogongrass_Cover = 68%...
## Predicting at Cogongrass_Cover = 69%...
## Predicting at Cogongrass_Cover = 70%...
## Predicting at Cogongrass_Cover = 71%...
## Predicting at Cogongrass_Cover = 72%...
## Predicting at Cogongrass_Cover = 73%...
## Predicting at Cogongrass_Cover = 74%...
## Predicting at Cogongrass_Cover = 75%...
## Predicting at Cogongrass_Cover = 76%...
## Predicting at Cogongrass_Cover = 77%...
## Predicting at Cogongrass_Cover = 78%...
## Predicting at Cogongrass_Cover = 79%...
## Predicting at Cogongrass_Cover = 80%...
## Predicting at Cogongrass_Cover = 81%...
## Predicting at Cogongrass_Cover = 82%...
## Predicting at Cogongrass_Cover = 83%...
## Predicting at Cogongrass_Cover = 84%...
## Predicting at Cogongrass_Cover = 85%...
## Predicting at Cogongrass_Cover = 86%...
## Predicting at Cogongrass_Cover = 87%...
## Predicting at Cogongrass_Cover = 88%...
## Predicting at Cogongrass_Cover = 89%...
## Predicting at Cogongrass_Cover = 90%...
## Predicting at Cogongrass_Cover = 91%...
## Predicting at Cogongrass_Cover = 92%...
## Predicting at Cogongrass_Cover = 93%...
## Predicting at Cogongrass_Cover = 94%...
## Predicting at Cogongrass_Cover = 95%...
## Predicting at Cogongrass_Cover = 96%...
## Predicting at Cogongrass_Cover = 97%...
## Predicting at Cogongrass_Cover = 98%...
## Predicting at Cogongrass_Cover = 99%...
## Predicting at Cogongrass_Cover = 100%...baseline <- prediction_list[["p0"]]
diff_list <- list()
cat("\nCalculating differences relative to 0% baseline...\n") ##
## Calculating differences relative to 0% baseline...for(level in scenario_levels[-1]) {
# Skip 0
diff_name <- paste0("Diff_", level,
"_vs_0")
diff_list[[diff_name]] <- prediction_list[[paste0("p", level)]] -
baseline
# Save result
writeRaster(diff_list[[diff_name]], paste0(diff_name, ".tif"), overwrite = TRUE)
}Summary Statistics Table
# Create the summary table
stats_summary <- lapply(names(diff_list), function(nm) {
v <- values(diff_list[[nm]], na.rm = TRUE)
data.frame(
Scenario = nm,
Mean_Change = mean(v),
SD_Change = sd(v),
Max_Loss = min(v),
Max_Gain = max(v)
)
}) %>%
bind_rows() %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))
# Create formatted table
ft <- flextable(stats_summary)
ft <- autofit(ft)
# Create Word document
doc <- read_docx()
doc <- body_add_par(
doc,
"Scenario Summary Statistics (Relative to Baseline)",
style = "heading 1"
)
doc <- body_add_flextable(doc, ft)
# Save to the output directory
output_file <- file.path(output_dir, "Scenario_Summary_Statistics.docx")
print(doc, target = output_file)
cat("Table saved to:\n", output_file, "\n")## Table saved to:
## C:/Users/alanivory34428/Desktop/03_Biodiversity/03_Outputs/Scenario_Summary_Statistics.docxPlot
plot_raster <- function(rast_layer, title, subtitle = NULL,
palette = "viridis", midpoint = NULL,
low = NULL, high = NULL, mid = NULL,
limits = NULL) {
# Convert SpatRaster
df <- as.data.frame(rast_layer, xy = TRUE)
colnames(df)[3] <- "value"
df <- df[!is.na(df$value), ]
p <- ggplot() +
geom_raster(data = df, aes(x = x, y = y, fill = value)) +
geom_sf(data = states_sf, fill = NA, colour = "grey30", linewidth = 0.3) +
coord_sf(xlim = c(-92, -75), ylim = c(24, 37)) +
labs(title = title, subtitle = subtitle, x = NULL, y = NULL) +
theme_bw(base_size = 12) +
theme(legend.position = "right",
axis.text = element_text(size = 8))
# Logic for diverging (diff) vs sequential (absolute) scales
if (!is.null(midpoint)) {
p <- p + scale_fill_gradient2(
low = low, mid = mid, high = high,
midpoint = midpoint,
limits = limits,
name = "ΔShannon",
na.value = "transparent"
)
} else {
p <- p + scale_fill_viridis_c(
option = palette,
limits = limits,
name = "Shannon\nDiversity",
na.value = "transparent"
)
}
return(p)
}
se_states <- c("Florida", "Georgia", "Alabama", "Mississippi", "South Carolina",
"North Carolina", "Tennessee", "Arkansas", "Louisiana", "Virginia")
states_sf <- maps::map("state", regions = tolower(se_states), fill = TRUE, plot = FALSE) %>%
st_as_sf()
p_final <- plot_raster(diff_list[["Diff_100_vs_0"]],
title = "Total Predicted Biodiversity Impact",
subtitle = "100% Cogongrass vs. 0% Baseline",
midpoint = 0,
low = "#d73027", mid = "white", high = "#1a9850")
print(p_final)## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
Thresholds Analysis
thresholds <- c(0.5, 0.6, 0.7)
threshold_rasters <- lapply(thresholds, function(thresh) {
diff_list[["Diff_100_vs_0"]] < -thresh
})
names(threshold_rasters) <- paste0("Loss > ", thresholds)
par(mfrow = c(1, 3), mar = c(3, 3, 3, 4))
for (i in seq_along(thresholds)) {
plot(
threshold_rasters[[i]],
main = paste0("Diversity Loss > ", thresholds[i]),
legend = FALSE,
col = c("grey90", "#d73027"),
axes = FALSE,
box = FALSE
)
v <- values(threshold_rasters[[i]], na.rm = TRUE)
mtext(
sprintf("%.1f%% of cells affected", mean(v, na.rm = TRUE) * 100),
side = 1, line = 1, cex = 0.8
)
}
ggsave(
plot = last_plot(),
filename = "threshold_analysis.png",
width = 10,
height = 4,
units = "in",
dpi = 300
)## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.par(mfrow = c(1, 1)) # reset layout
plot_threshold_maps <- function(diff_raster,
diff_name,
thresholds = c(0.5, 0.6, 0.7),
out_dir = ".") {
library(terra)
# Extract cover percentage from the raster name
cover <- as.numeric(sub("Diff_([0-9]+)_vs_0", "\\1", diff_name))
# Create threshold rasters
threshold_rasters <- lapply(thresholds, function(thresh) {
diff_raster < -thresh
})
letters <- LETTERS[seq_along(thresholds)]
## Save PNG
png(
file.path(out_dir,
paste0(diff_name, "_SHDI_loss_thresholds.png")),
width = 12,
height = 4,
units = "in",
res = 300
)
par(mfrow = c(1, length(thresholds)),
mar = c(3, 3, 4, 2))
for (i in seq_along(thresholds)) {
plot(
threshold_rasters[[i]],
main = sprintf(
"%s) %d%% Cover \u2014 SHDI Loss > %.2f",
letters[i],
cover,
thresholds[i]
),
legend = FALSE,
col = c("grey90", "#d73027"),
axes = FALSE,
box = FALSE
)
}
dev.off()
## Save PDF
pdf(
file.path(out_dir,
paste0(diff_name, "_SHDI_loss_thresholds.pdf")),
width = 12,
height = 4
)
par(mfrow = c(1, length(thresholds)),
mar = c(3, 3, 4, 2))
for (i in seq_along(thresholds)) {
plot(
threshold_rasters[[i]],
main = sprintf(
"%s) %d%% Cover \u2014 SHDI Loss > %.2f",
letters[i],
cover,
thresholds[i]
),
legend = FALSE,
col = c("grey90", "#d73027"),
axes = FALSE,
box = FALSE
)
}
dev.off()
}Varying Threshold and Raster
thresholds <- c(0.25, 0.5, 0.75)
diff_layers <- c("Diff_50_vs_0", "Diff_75_vs_0", "Diff_100_vs_0")
cover_labels <- c(
"Diff_50_vs_0" = "50% Cover",
"Diff_75_vs_0" = "75% Cover",
"Diff_100_vs_0" = "100% Cover"
)
# --- State outlines (same approach as your other script) ---
se_states <- c("Florida", "Georgia", "Alabama", "Mississippi", "South Carolina",
"North Carolina", "Tennessee", "Arkansas", "Louisiana", "Virginia")
states_sf <- maps::map("state", regions = tolower(se_states), fill = TRUE, plot = FALSE) %>%
st_as_sf()
# --- Build a threshold raster + subplot for every cover% x threshold combo ---
# NOTE: letters are assigned in the same nested-loop order as the original
# script (diff_layers outer, thresholds inner), so panel A-I map to:
# A) 50% Cover: D 0.25 SHDI D) 75% Cover: D 0.25 SHDI G) 100% Cover: D 0.25 SHDI
# B) 50% Cover: D 0.5 SHDI E) 75% Cover: D 0.5 SHDI H) 100% Cover: D 0.5 SHDI
# C) 50% Cover: D 0.75 SHDI F) 75% Cover: D 0.75 SHDI I) 100% Cover: D 0.75 SHDI
# Swap the order of the two loops below if you need a different pairing.
panel_letters <- LETTERS[seq_len(length(diff_layers) * length(thresholds))]
threshold_rasters <- list()
plot_list <- list()
cell_stats <- data.frame(
panel = character(),
cover = character(),
threshold = numeric(),
n_affected = integer(),
n_total = integer(),
pct_affected = numeric(),
stringsAsFactors = FALSE
)
i <- 1
for (dl in diff_layers) {
for (thresh in thresholds) {
key <- paste0(dl, "_gt_", thresh)
r <- diff_list[[dl]] < -thresh
threshold_rasters[[key]] <- r
# count + percent of cells affected -> now goes into the table, not onto the plot
v <- terra::values(r, na.rm = TRUE)
n_total <- length(v)
n_aff <- sum(v, na.rm = TRUE)
pct <- mean(v, na.rm = TRUE) * 100
letter <- panel_letters[i]
panel_title <- sprintf("%s) %s: \u0394 %.2g SHDI", letter, cover_labels[[dl]], thresh)
cell_stats <- rbind(cell_stats, data.frame(
panel = letter,
cover = cover_labels[[dl]],
threshold = thresh,
n_affected = n_aff,
n_total = n_total,
pct_affected = round(pct, 1)
))
df <- as.data.frame(r, xy = TRUE)
colnames(df)[3] <- "value"
df <- df[!is.na(df$value), ]
df$value <- factor(df$value, levels = c(FALSE, TRUE))
p <- ggplot() +
geom_raster(data = df, aes(x = x, y = y, fill = value)) +
geom_sf(data = states_sf, fill = NA, colour = "grey30", linewidth = 0.3) +
coord_sf(xlim = c(-91, -74), ylim = c(24, 37), datum = NA) +
scale_fill_manual(values = c("FALSE" = "grey90", "TRUE" = "#FA4616"), guide = "none") +
labs(title = panel_title, x = NULL, y = NULL) +
theme_bw(base_size = 10) +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 10, face = "bold")
)
plot_list[[key]] <- p
i <- i + 1
}
}
# --- Combine the 9 panels into one 3x3 figure ---
final_plot <- wrap_plots(plot_list, ncol = 3, nrow = 3)
ggsave(
filename = file.path(output_dir, "threshold_raster_analysis.png"),
plot = final_plot,
width = 12,
height = 12,
units = "in",
dpi = 300
)## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.# --- Separate table of % cells affected (no longer printed on the plot) ---
print(cell_stats)## panel cover threshold n_affected n_total pct_affected
## 1 A 50% Cover 0.25 0 411926 0.0
## 2 B 50% Cover 0.50 0 411926 0.0
## 3 C 50% Cover 0.75 0 411926 0.0
## 4 D 75% Cover 0.25 2 411926 0.0
## 5 E 75% Cover 0.50 0 411926 0.0
## 6 F 75% Cover 0.75 0 411926 0.0
## 7 G 100% Cover 0.25 407258 411926 98.9
## 8 H 100% Cover 0.50 209885 411926 51.0
## 9 I 100% Cover 0.75 158 411926 0.0write.csv(
cell_stats,
file = file.path(output_dir, "threshold_cell_stats.csv"),
row.names = FALSE
)EPA Ecoregions
eco_l2 <- st_read("C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/NA_CEC_Eco_Level2.shp") %>% st_transform(crs(diff_list[[1]]))## Reading layer `NA_CEC_Eco_Level2' from data source
## `C:\Users\alanivory34428\Desktop\03_Biodiversity\02_Data\NA_CEC_Eco_Level2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2261 features and 8 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Areaeco_l3 <- st_read("C:/Users/alanivory34428/Desktop/03_Biodiversity/02_Data/NA_CEC_Eco_Level3.shp") %>% st_transform(crs(diff_list[[1]])) ## Reading layer `NA_CEC_Eco_Level3' from data source
## `C:\Users\alanivory34428\Desktop\03_Biodiversity\02_Data\NA_CEC_Eco_Level3.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2548 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Area# Convert to SpatVector for terra::extract
eco_l2_vect <- terra::vect(eco_l2)
eco_l3_vect <- terra::vect(eco_l3)
summarize_by_ecoregion <- function(diff_list, eco_vect, name_col) {
# Rasterize polygons
template_raster <- diff_list[[1]]
eco_rast <- terra::rasterize(eco_vect, template_raster, field = name_col)
lapply(names(diff_list), function(nm) {
r <- diff_list[[nm]]
z <- terra::zonal(r, eco_rast, fun = "mean", na.rm = TRUE)
colnames(z) <- c("Ecoregion", "Mean_Change")
z$Scenario <- nm
z
}) %>%
bind_rows() %>%
filter(!is.na(Mean_Change)) %>%
mutate(
Scenario_num = as.numeric(gsub("Diff_|_vs_0", "", Scenario))
)
}
# names(as.data.frame(eco_l2_vect))
# names(as.data.frame(eco_l3_vect))
# Common EPA shapefile name columns:
# Level 2 -> "NA_L2NAME"
# Level 3 -> "US_L3NAME"
eco_l2_summary <- summarize_by_ecoregion(diff_list, eco_l2_vect, name_col = "NA_L2NAME")
eco_l3_summary <- summarize_by_ecoregion(diff_list, eco_l3_vect, name_col = "NA_L3NAME")Plot EPA
# L3 Ecoregion
eco_l3_100 <- eco_l3_summary %>%
dplyr::filter(Scenario_num == 100)
eco_l3_map <- eco_l3 %>%
dplyr::left_join(eco_l3_100, by = c("NA_L3NAME" = "Ecoregion"))
eco_l3_map_vect <- terra::vect(eco_l3_map)
template_raster <- diff_list[[1]]
eco_l3_map_vect <- terra::crop(eco_l3_map_vect, template_raster)
terra::plot(
eco_l3_map_vect,
"Mean_Change",
col = hcl.colors(50, "RdYlGn", rev = TRUE),
main = "Δ Shannon Diversity (100% Cogongrass)"
)
ggplot(eco_l3_map) +
geom_sf(aes(fill = Mean_Change), color = NA) +
scale_fill_gradient2(
low = "#1a9850",
mid = "white",
high = "#d73027",
midpoint = 0
) +
theme_bw() +
labs(
title = "Δ Shannon Diversity (100% Cogongrass)",
fill = "Change"
)
#L2 Ecoregion
eco_l2_100 <- eco_l2_summary %>%
dplyr::filter(Scenario_num == 100)
eco_l2_map <- eco_l2 %>%
dplyr::left_join(eco_l2_100, by = c("NA_L2NAME" = "Ecoregion"))
eco_l2_map_vect <- terra::vect(eco_l2_map)
template_raster <- diff_list[[1]]
eco_l2_map_vect <- terra::crop(eco_l2_map_vect, template_raster)
terra::plot(
eco_l2_map_vect,
"Mean_Change",
col = hcl.colors(50, "RdYlGn", rev = TRUE),
main = "Δ Shannon Diversity (Level 2 Ecoregions, 100% Cogongrass)"
)
ggplot(eco_l2_map) +
geom_sf(aes(fill = Mean_Change), color = "black", linewidth = 0.2) +
scale_fill_gradient2(
low = "#1a9850",
mid = "white",
high = "#d73027",
midpoint = 0
) +
theme_bw() +
labs(
title = "Δ Shannon Diversity (Level 2 Ecoregions, 100% Cogongrass)",
fill = "Change"
)
## two panel figure
# Common color limits so both maps are directly comparable
lims <- range(
c(eco_l2_map$Mean_Change, eco_l3_map$Mean_Change),
na.rm = TRUE
)
# -----------------------------
# Level III Ecoregions
# -----------------------------
p_l3 <- ggplot(eco_l3_map) +
geom_sf(aes(fill = Mean_Change), color = "black", linewidth = 0.1) +
scale_fill_gradient2(
low = "#1a9850",
mid = "white",
high = "#d73027",
midpoint = 0,
limits = lims,
name = expression(Delta*" Shannon")
) +
coord_sf(xlim = c(-92, -75), ylim = c(24, 37), expand = FALSE) +
labs(title = "Level III Ecoregions") +
theme_bw(base_size = 12) +
theme(
panel.grid = element_blank(),
axis.title = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5)
)
# -----------------------------
# Level II Ecoregions
# -----------------------------
p_l2 <- ggplot(eco_l2_map) +
geom_sf(aes(fill = Mean_Change), color = "black", linewidth = 0.2) +
scale_fill_gradient2(
low = "#1a9850",
mid = "white",
high = "#d73027",
midpoint = 0,
limits = lims,
name = expression(Delta*" Shannon")
) +
coord_sf(xlim = c(-92, -75), ylim = c(24, 37), expand = FALSE) +
labs(title = "Level II Ecoregions") +
theme_bw(base_size = 12) +
theme(
panel.grid = element_blank(),
axis.title = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5)
)
eco_panel <-
p_l3 + p_l2 +
plot_layout(guides = "collect") +
plot_annotation(
title = "Predicted Change in Shannon Diversity Under 100% Cogongrass",
tag_levels = "A"
) &
theme(
legend.position = "right",
plot.tag = element_text(face = "bold", size = 14)
)
eco_panel
# Save
ggsave(
filename = "Ecoregion_Shannon_Change_100Cogongrass.png",
plot = eco_panel,
width = 12,
height = 6,
dpi = 600,
bg = "white"
)thresholds <- c(0.5, 0.6, 0.7)
threshold_rasters <- lapply(thresholds, function(thresh) {
diff_list[["Diff_100_vs_0"]] < -thresh
})
names(threshold_rasters) <- paste0("Loss > ", thresholds)
## Quick on-screen panel + saved PNG (base graphics -> capture with png(), NOT ggsave)
png(file.path(output_dir, "threshold_analysis.png"),
width = 10, height = 4, units = "in", res = 300)
par(mfrow = c(1, 3), mar = c(3, 3, 3, 4))
for (i in seq_along(thresholds)) {
plot(
threshold_rasters[[i]],
main = paste0("Diversity Loss > ", thresholds[i]),
legend = FALSE,
col = c("grey90", "#FA4616"),
axes = FALSE,
box = FALSE
)
v <- values(threshold_rasters[[i]], na.rm = TRUE)
mtext(
sprintf("%.1f%% of cells affected", mean(v, na.rm = TRUE) * 100),
side = 1, line = 1, cex = 0.8
)
}
dev.off()## png
## 2par(mfrow = c(1, 1)) # reset layout
## Reusable function: saves a PNG + PDF panel of loss thresholds for one diff raster.
## out_dir defaults to the global output_dir defined in the setup chunk.
plot_threshold_maps <- function(diff_raster,
diff_name,
thresholds = c(0.5, 0.6, 0.7),
out_dir = output_dir) {
library(terra)
# Extract cover percentage from the raster name
cover <- as.numeric(sub("Diff_([0-9]+)_vs_0", "\\1", diff_name))
# Create threshold rasters
threshold_rasters <- lapply(thresholds, function(thresh) {
diff_raster < -thresh
})
panel_letters <- LETTERS[seq_along(thresholds)]
draw_panel <- function() {
par(mfrow = c(1, length(thresholds)), mar = c(3, 3, 4, 2))
for (i in seq_along(thresholds)) {
plot(
threshold_rasters[[i]],
main = sprintf(
"%s) %d%% Cover \u2014 SHDI Loss > %.2f",
panel_letters[i], cover, thresholds[i]
),
legend = FALSE,
col = c("grey90", "#FA4616"),
axes = FALSE,
box = FALSE
)
}
}
## Save PNG
png(file.path(out_dir, paste0(diff_name, "_SHDI_loss_thresholds.png")),
width = 12, height = 4, units = "in", res = 300)
draw_panel()
dev.off()
## Save PDF
pdf(file.path(out_dir, paste0(diff_name, "_SHDI_loss_thresholds.pdf")),
width = 12, height = 4)
draw_panel()
dev.off()
}
## Example: plot_threshold_maps(diff_list[["Diff_100_vs_0"]], "Diff_100_vs_0")New Figures
## ============================================================
## Objective 3 figures, rebuilt around what the scenarios actually show
##
## The 3x3 threshold panel has six empty panels because nothing crosses
## 0.25 SHDI until ~87% cover. The scenario series locates the change as a
## step between 86% and 87%: mean goes from -0.234 to -0.486, five times
## larger than any other one-percent increment. That is the finding, and
## the figure should show it.
##
## Requires: diff_list, states_sf, output_dir (from the prediction Rmd)
## ============================================================
pacman::p_load(terra, sf, tidyverse, patchwork, maps)
gator_orange <- "#FA4616"
gator_blue <- "#0021A5"
## ---- 0. Correct cell area -----------------------------------------------
## The reference grid is the NDVI raster: ~0.009 degrees, i.e. roughly 1 km,
## NOT 30 m as the alignment comments say. In a lon/lat CRS cell area also
## varies with latitude, so use cellSize() rather than a constant.
cell_km2 <- terra::cellSize(diff_list[[1]], unit = "km")
total_km2 <- as.numeric(terra::global(
terra::mask(cell_km2, diff_list[[1]]), "sum", na.rm = TRUE))
cat(sprintf("Study extent: %s cells, %.0f km2 (%.1f million ha)\n",
format(as.numeric(terra::global(!is.na(diff_list[[1]]), "sum",
na.rm = TRUE)), big.mark = ","),
total_km2, total_km2 * 100 / 1e6))## Study extent: 411,926 cells, 348177 km2 (34.8 million ha)## ---- 1. The scenario curve: this is the key panel ------------------------
scen <- map_dfr(names(diff_list), function(nm) {
v <- terra::values(diff_list[[nm]], na.rm = TRUE)
tibble(cover = as.numeric(gsub("Diff_|_vs_0", "", nm)),
mean = mean(v), sd = sd(v),
min = min(v), max = max(v))
}) %>% arrange(cover)
## Locate the transition: the single largest one-percent drop
steps <- scen %>% mutate(delta = mean - lag(mean)) %>% filter(!is.na(delta))
brk <- steps$cover[which.min(steps$delta)]
cat(sprintf("Largest one-percent drop at %d%% cover: %.3f units\n",
brk, min(steps$delta, na.rm = TRUE)))## Largest one-percent drop at 87% cover: -0.251 unitscat(sprintf(" next largest is %.3f, so the step is %.1f times larger\n",
sort(steps$delta)[2], min(steps$delta) / sort(steps$delta)[2]))## next largest is -0.051, so the step is 4.9 times larger## Cover above which predictions stop changing = the observed data ceiling
flat <- scen %>% filter(abs(mean - last(scen$mean)) < 0.001) %>% pull(cover) %>% min()
cat(sprintf("Predictions flat from %d%% cover onward (RF cannot extrapolate)\n", flat))## Predictions flat from 92% cover onward (RF cannot extrapolate)p_curve <- ggplot(scen, aes(cover, mean)) +
geom_hline(yintercept = 0, colour = "grey60", linewidth = 0.3) +
geom_ribbon(aes(ymin = mean - sd, ymax = mean + sd),
fill = gator_orange, alpha = 0.25) +
geom_line(colour = gator_blue, linewidth = 1.1) +
geom_vline(xintercept = brk, linetype = "dashed", colour = "grey30") +
annotate("text", x = brk - 2, y = min(scen$mean) * 0.55,
label = sprintf("%d%%", brk), hjust = 1, size = 3.5, fontface = "bold") +
annotate("rect", xmin = flat, xmax = 100, ymin = -Inf, ymax = Inf,
fill = "grey50", alpha = 0.12) +
annotate("text", x = (flat + 100) / 2, y = max(scen$mean + scen$sd),
label = "beyond\nobserved\nrange", size = 2.7, colour = "grey40", vjust = 1) +
labs(x = "Cogongrass cover (%)",
y = expression(Delta*" Shannon diversity"),
title = "A) Range-wide mean change") +
theme_classic(base_size = 11) +
theme(plot.title = element_text(face = "bold", size = 10))
## ---- 2. Maps that are not empty ------------------------------------------
## Bracket the transition instead of showing 50/75/100, where two rows are blank.
map_panel <- function(diff_name, thresh, title) {
r <- diff_list[[diff_name]] < -thresh
df <- as.data.frame(r, xy = TRUE)
colnames(df)[3] <- "value"
df <- df[!is.na(df$value), ]
df$value <- factor(df$value, levels = c(FALSE, TRUE))
pct <- 100 * mean(terra::values(r, na.rm = TRUE), na.rm = TRUE)
ggplot() +
geom_raster(data = df, aes(x, y, fill = value)) +
geom_sf(data = states_sf, fill = NA, colour = "grey30", linewidth = 0.25) +
coord_sf(xlim = c(-95, -75), ylim = c(24, 37), datum = NA) +
scale_fill_manual(values = c("FALSE" = "grey90", "TRUE" = gator_orange),
guide = "none") +
labs(title = sprintf("%s (%.1f%%)", title, pct), x = NULL, y = NULL) +
theme_bw(base_size = 10) +
theme(axis.text = element_blank(), axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 9, face = "bold"))
}
## Row 1: same threshold, cover levels that bracket the step
row1 <- list(
map_panel("Diff_80_vs_0", 0.25, "B) 80% cover, loss > 0.25"),
map_panel(paste0("Diff_", brk, "_vs_0"), 0.25, sprintf("C) %d%% cover, loss > 0.25", brk)),
map_panel(paste0("Diff_", flat, "_vs_0"), 0.25, sprintf("D) %d%% cover, loss > 0.25", flat))
)
## Row 2: at the observed ceiling, vary the threshold - this is where the
## spatial signal lives, since almost everything crosses 0.25
row2 <- list(
map_panel(paste0("Diff_", flat, "_vs_0"), 0.25, sprintf("E) %d%% cover, loss > 0.25", flat)),
map_panel(paste0("Diff_", flat, "_vs_0"), 0.50, sprintf("F) %d%% cover, loss > 0.50", flat)),
map_panel(paste0("Diff_", flat, "_vs_0"), 0.75, sprintf("G) %d%% cover, loss > 0.75", flat))
)
fig6 <- p_curve / wrap_plots(row1, nrow = 1) / wrap_plots(row2, nrow = 1) +
plot_layout(heights = c(1, 1, 1))
ggsave(file.path(output_dir, "Figure6_threshold_and_maps.png"),
fig6, width = 10, height = 11, units = "in", dpi = 300, bg = "white")## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.## ---- 3. Onset map: the cover at which each cell first crosses -0.25 ------
## This is the direct answer to Victoria's comments 111, 113, 114 and 116 -
## it shows spatial variation in a single map instead of asking the reader to
## compare panels. Cells that never cross stay NA.
cat("\nBuilding onset map (this loops over all scenarios)...\n")##
## Building onset map (this loops over all scenarios)...onset <- diff_list[[1]] * NA
for (nm in names(diff_list)) {
cov_i <- as.numeric(gsub("Diff_|_vs_0", "", nm))
crossed <- diff_list[[nm]] <= -0.25
onset <- terra::ifel(is.na(onset) & crossed, cov_i, onset)
}
names(onset) <- "onset_cover"
onset_df <- as.data.frame(onset, xy = TRUE)
colnames(onset_df)[3] <- "cover"
onset_df <- onset_df[!is.na(onset_df$cover), ]
p_onset <- ggplot() +
geom_raster(data = onset_df, aes(x, y, fill = cover)) +
geom_sf(data = states_sf, fill = NA, colour = "grey30", linewidth = 0.25) +
coord_sf(xlim = c(-95, -75), ylim = c(24, 37), datum = NA) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "Cogongrass\ncover (%)") +
labs(title = "Cover at which predicted diversity first falls by 0.25 units",
x = NULL, y = NULL) +
theme_bw(base_size = 11) +
theme(axis.text = element_blank(), axis.ticks = element_blank(),
panel.grid = element_blank())
ggsave(file.path(output_dir, "Figure_onset_cover.png"),
p_onset, width = 8, height = 6, units = "in", dpi = 300, bg = "white")## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.cat("\nOnset cover summary:\n")##
## Onset cover summary:print(summary(onset_df$cover))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 72.00 84.00 87.00 85.36 87.00 90.00## ---- 4. Ecoregion table for the Results text ----------------------------
## eco_l3_summary is built earlier in the prediction Rmd.
if (exists("eco_l3_summary")) {
eco_top <- eco_l3_summary %>%
filter(Scenario_num == flat) %>%
arrange(Mean_Change) %>%
mutate(Mean_Change = round(Mean_Change, 3))
cat(sprintf("\n===== Ecoregions ranked by mean loss at %d%% cover =====\n", flat))
print(as.data.frame(head(eco_top, 10)))
cat("\n--- least affected ---\n")
print(as.data.frame(tail(eco_top, 5)))
}##
## ===== Ecoregions ranked by mean loss at 92% cover =====
## Ecoregion Mean_Change Scenario Scenario_num
## 1 Ridge and Valley -0.559 Diff_92_vs_0 92
## 2 South Central Plains -0.532 Diff_92_vs_0 92
## 3 Southwestern Appalachians -0.526 Diff_92_vs_0 92
## 4 Mississippi Alluvial Plain -0.525 Diff_92_vs_0 92
## 5 Mississippi Valley Loess Plains -0.510 Diff_92_vs_0 92
## 6 Southern Coastal Plain -0.496 Diff_92_vs_0 92
## 7 Southeastern Plains -0.486 Diff_92_vs_0 92
## 8 Western Gulf Coastal Plain -0.482 Diff_92_vs_0 92
## 9 Piedmont -0.478 Diff_92_vs_0 92
## 10 Blue Ridge -0.468 Diff_92_vs_0 92
##
## --- least affected ---
## Ecoregion Mean_Change Scenario Scenario_num
## 7 Southeastern Plains -0.486 Diff_92_vs_0 92
## 8 Western Gulf Coastal Plain -0.482 Diff_92_vs_0 92
## 9 Piedmont -0.478 Diff_92_vs_0 92
## 10 Blue Ridge -0.468 Diff_92_vs_0 92
## 11 Middle Atlantic Coastal Plain -0.442 Diff_92_vs_0 92