Reworked Biodiversity Prediction
Alan Ivory
2026-07-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" "pr" "sph"
## [5] "srad" "TreeCanopy_NLCD" "ndvi" "avg_pop_den_1km"env_vars <- model_vars
cat(sprintf("\nEnvironmental predictors : %s\n", paste(env_vars, collapse = ", ")))##
## Environmental predictors : Cogongrass_Cover, elevation, pr, sph, srad, TreeCanopy_NLCD, 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"
)
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: pr
## Loading: sph
## Loading: srad
## Loading: TreeCanopy_NLCD
## 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
## pr 4326
## sph 4326
## srad 4326
## TreeCanopy_NLCD 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: NAnames(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...
## |---------|---------|---------|---------|========================================= 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: 10 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 10514534NLCD 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
cat("\n========== Scenario Summary Statistics (Relative to Baseline) ==========\n")##
## ========== Scenario Summary Statistics (Relative to Baseline) ==========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()
print(stats_summary)## Scenario Mean_Change SD_Change Max_Loss Max_Gain
## 1 Diff_1_vs_0 0.0010213375 0.001861690 -0.007701703 0.01886696
## 2 Diff_2_vs_0 -0.0008641458 0.004617821 -0.030605133 0.02513047
## 3 Diff_3_vs_0 -0.0022970856 0.005626169 -0.031955825 0.02569142
## 4 Diff_4_vs_0 -0.0034163428 0.006911026 -0.034865115 0.02578952
## 5 Diff_5_vs_0 -0.0020452065 0.007319406 -0.037114557 0.03552967
## 6 Diff_6_vs_0 -0.0012436414 0.009460904 -0.048444702 0.03943584
## 7 Diff_7_vs_0 0.0008839535 0.011550145 -0.049212279 0.05954409
## 8 Diff_8_vs_0 0.0028398687 0.012665343 -0.049077183 0.07059618
## 9 Diff_9_vs_0 0.0047368305 0.014859366 -0.048343908 0.09670917
## 10 Diff_10_vs_0 0.0048793228 0.015247550 -0.049525306 0.09657862
## 11 Diff_11_vs_0 0.0050263784 0.019803345 -0.068899976 0.11083684
## 12 Diff_12_vs_0 0.0060173562 0.021247987 -0.079794582 0.12217427
## 13 Diff_13_vs_0 0.0070660719 0.021647881 -0.080592944 0.12375827
## 14 Diff_14_vs_0 0.0058332790 0.023405847 -0.092562239 0.12678151
## 15 Diff_15_vs_0 0.0062469432 0.024095411 -0.094338721 0.13314818
## 16 Diff_16_vs_0 0.0045913567 0.024774333 -0.101627664 0.13343950
## 17 Diff_17_vs_0 0.0043063303 0.025167435 -0.105587769 0.13735877
## 18 Diff_18_vs_0 0.0021972312 0.025720140 -0.106557257 0.13672746
## 19 Diff_19_vs_0 -0.0007159006 0.028505610 -0.109777571 0.13704183
## 20 Diff_20_vs_0 0.0007526056 0.028183657 -0.107769950 0.13852169
## 21 Diff_21_vs_0 -0.0149947215 0.038399987 -0.191695098 0.12947139
## 22 Diff_22_vs_0 -0.0123990267 0.040847306 -0.198325316 0.13970027
## 23 Diff_23_vs_0 -0.0116849146 0.040776369 -0.200017303 0.13878561
## 24 Diff_24_vs_0 -0.0061536333 0.044210978 -0.200443834 0.16457902
## 25 Diff_25_vs_0 -0.0062037334 0.046278651 -0.211993547 0.17854369
## 26 Diff_26_vs_0 -0.0063663099 0.046699047 -0.213418565 0.18167894
## 27 Diff_27_vs_0 -0.0069887478 0.047678010 -0.218194000 0.18875150
## 28 Diff_28_vs_0 -0.0061911680 0.047818486 -0.217153256 0.19500866
## 29 Diff_29_vs_0 -0.0039202039 0.049165289 -0.222023531 0.19918652
## 30 Diff_30_vs_0 -0.0026326584 0.048902794 -0.219707626 0.19997227
## 31 Diff_31_vs_0 0.0005873017 0.048529942 -0.218397092 0.19984448
## 32 Diff_32_vs_0 0.0095414800 0.051244659 -0.217011236 0.20775359
## 33 Diff_33_vs_0 0.0094939428 0.051514705 -0.219742078 0.20977435
## 34 Diff_34_vs_0 0.0108112982 0.052138185 -0.220650905 0.21390007
## 35 Diff_35_vs_0 0.0115279664 0.052094839 -0.220495020 0.21366785
## 36 Diff_36_vs_0 0.0085939829 0.055019176 -0.237198662 0.21116477
## 37 Diff_37_vs_0 0.0103821338 0.055352407 -0.236160410 0.21336591
## 38 Diff_38_vs_0 0.0115885616 0.055693576 -0.235448398 0.21809616
## 39 Diff_39_vs_0 0.0192220256 0.056653985 -0.229471021 0.23044641
## 40 Diff_40_vs_0 0.0355293999 0.061163383 -0.221677157 0.26103981
## 41 Diff_41_vs_0 0.0456793395 0.064283281 -0.217161538 0.27767963
## 42 Diff_42_vs_0 0.0472482967 0.065053767 -0.216552544 0.28092681
## 43 Diff_43_vs_0 0.0486794139 0.066066823 -0.216446702 0.28577655
## 44 Diff_44_vs_0 0.0488099987 0.066456321 -0.215304206 0.28657841
## 45 Diff_45_vs_0 0.0560255622 0.070973828 -0.215366797 0.31301367
## 46 Diff_46_vs_0 0.0645265319 0.074684731 -0.211153479 0.33169245
## 47 Diff_47_vs_0 0.0795249518 0.080329273 -0.204317608 0.36582370
## 48 Diff_48_vs_0 0.0794855412 0.080381926 -0.204091955 0.36582370
## 49 Diff_49_vs_0 0.0789039935 0.080552089 -0.204296854 0.36328152
## 50 Diff_50_vs_0 0.0807898782 0.080282282 -0.201530620 0.36463211
## 51 Diff_51_vs_0 0.0800081729 0.079787685 -0.202125930 0.36346601
## 52 Diff_52_vs_0 0.0769695531 0.079479071 -0.203994742 0.35888934
## 53 Diff_53_vs_0 0.0763164570 0.079509948 -0.206847094 0.35967520
## 54 Diff_54_vs_0 0.0766752009 0.079632955 -0.207028951 0.35983005
## 55 Diff_55_vs_0 0.0762825335 0.080057582 -0.210197121 0.35921672
## 56 Diff_56_vs_0 0.0760574783 0.080191618 -0.212063524 0.35921672
## 57 Diff_57_vs_0 0.0771187162 0.080424710 -0.210871863 0.36059436
## 58 Diff_58_vs_0 0.0770907339 0.080296647 -0.210871863 0.36003046
## 59 Diff_59_vs_0 0.0815083208 0.081262391 -0.210170471 0.36853038
## 60 Diff_60_vs_0 0.0816120578 0.081227409 -0.210170471 0.36853038
## 61 Diff_61_vs_0 0.0830372918 0.082682017 -0.211346589 0.37089878
## 62 Diff_62_vs_0 0.0825901131 0.084084481 -0.215286409 0.37189814
## 63 Diff_63_vs_0 0.0869040605 0.086045570 -0.212022934 0.38080909
## 64 Diff_64_vs_0 0.0877776089 0.086443514 -0.213006372 0.38230923
## 65 Diff_65_vs_0 0.0886523027 0.087963544 -0.210710996 0.38493751
## 66 Diff_66_vs_0 0.0879734465 0.087531124 -0.210710996 0.38252104
## 67 Diff_67_vs_0 0.0876801799 0.087603037 -0.210710996 0.38252104
## 68 Diff_68_vs_0 0.0832809590 0.088346702 -0.218580457 0.38252104
## 69 Diff_69_vs_0 0.0832590543 0.088736114 -0.218151694 0.38129775
## 70 Diff_70_vs_0 0.0758684759 0.092557015 -0.242969883 0.37986219
## 71 Diff_71_vs_0 0.0757585021 0.092488051 -0.242969883 0.37986219
## 72 Diff_72_vs_0 0.0731098397 0.093046928 -0.247482725 0.37905681
## 73 Diff_73_vs_0 0.0505145507 0.103730851 -0.302407189 0.37436766
## 74 Diff_74_vs_0 0.0341543520 0.109044768 -0.336114407 0.36887995
## 75 Diff_75_vs_0 0.0301185588 0.109598350 -0.340534597 0.36881159
## 76 Diff_76_vs_0 0.0287373267 0.108989978 -0.340242420 0.36479731
## 77 Diff_77_vs_0 0.0280713553 0.109000391 -0.340242420 0.36479731
## 78 Diff_78_vs_0 0.0272104666 0.108418702 -0.340242420 0.36194245
## 79 Diff_79_vs_0 0.0159004516 0.107078187 -0.352064209 0.34795824
## 80 Diff_80_vs_0 -0.0123158585 0.104470174 -0.375664229 0.31747742
## 81 Diff_81_vs_0 -0.0129292108 0.104220989 -0.377589096 0.31747742
## 82 Diff_82_vs_0 -0.0129292108 0.104220989 -0.377589096 0.31747742
## 83 Diff_83_vs_0 -0.0672459041 0.100566954 -0.419646346 0.26749144
## 84 Diff_84_vs_0 -0.0715929178 0.100958635 -0.422996055 0.26576706
## 85 Diff_85_vs_0 -0.0715929178 0.100958635 -0.422996055 0.26576706
## 86 Diff_86_vs_0 -0.0715929178 0.100958635 -0.422996055 0.26576706
## 87 Diff_87_vs_0 -0.4698887035 0.098566288 -0.852807458 -0.17817782
## 88 Diff_88_vs_0 -0.5019176038 0.099448898 -0.894780528 -0.19835943
## 89 Diff_89_vs_0 -0.5019176038 0.099448898 -0.894780528 -0.19835943
## 90 Diff_90_vs_0 -0.5106966681 0.101054324 -0.899419291 -0.20447579
## 91 Diff_91_vs_0 -0.5096832280 0.101054324 -0.898405851 -0.20346235
## 92 Diff_92_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573
## 93 Diff_93_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573
## 94 Diff_94_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573
## 95 Diff_95_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573
## 96 Diff_96_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573
## 97 Diff_97_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573
## 98 Diff_98_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573
## 99 Diff_99_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573
## 100 Diff_100_vs_0 -0.4960097137 0.101124466 -0.884329236 -0.18938573Plot
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")
# Threshold raster for every combination of diff layer x threshold
threshold_rasters <- list()
for (dl in diff_layers) {
for (thresh in thresholds) {
key <- paste0(dl, "_gt_", thresh)
threshold_rasters[[key]] <- diff_list[[dl]] < -thresh
}
}
par(mfrow = c(3, 3), mar = c(3, 3, 3, 4))
for (dl in diff_layers) {
for (thresh in thresholds) {
key <- paste0(dl, "_gt_", thresh)
r <- threshold_rasters[[key]]
plot(
r,
main = paste0(dl, ": Loss > ", thresh),
legend = FALSE,
col = c("grey90", "#d73027"),
axes = FALSE,
box = FALSE
)
v <- values(r, 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_raster_analysis.png",
width = 10,
height = 10,
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))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", "#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
)
}
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", "#d73027"),
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")