Reworked Biodiversity Prediction
Alan Ivory
2026-04-09
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" "srad"
## [5] "TreeCanopy_NLCD" "avg_pop_den_1km" "MEM2" "MEM5"
## [9] "MEM23" "MEM118" "MEM203"# Then separate env vs MEM vars as before
env_vars <- model_vars[!grepl("^MEM", model_vars)]
mem_vars <- model_vars[ grepl("^MEM", model_vars)]
cat(sprintf("\nEnvironmental predictors : %s\n", paste(env_vars, collapse = ", ")))##
## Environmental predictors : Cogongrass_Cover, elevation, pr, srad, TreeCanopy_NLCD, avg_pop_den_1kmcat(sprintf("MEM predictors (→ set 0) : %s\n\n", paste(mem_vars, collapse = ", ")))## MEM predictors (→ set 0) : MEM2, MEM5, MEM23, MEM118, MEM203NLCD
# After building env_vars, exclude Cogongrass_Cover since it's not a raster predictor but a fixed value in scenarios
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 1.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: srad
## Loading: TreeCanopy_NLCD
## 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
# --- Check what CRS your rasters are in ---
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
## srad 4326
## TreeCanopy_NLCD 5070
## 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 (reprojecting extent if needed) ---
local_dir <- "C:/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
I will change ref_rast to pr for now, to make predictions at the 4km scale. I can switch back to ndvi if we decide 30 m resolution is needed. Because weather data is at 4km, weather pixels become roughly 17,000 identical pixels.
# 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 necessary
if(!compareGeom(nlcd_raster, env_stack, stopOnError = FALSE)) {
nlcd_raster <- project(nlcd_raster, env_stack, method = "near")
}
# Define 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...
## |---------|---------|---------|---------|========================================= Predict Function
# Takes a cogongrass cover value, builds predictor dataframe,
# sets MEMs to 0, predicts, returns a raster
predict_diversity <- function(cogongrass_value, env_stack, rf_model, mem_vars) {
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)
# Logic: Only predict for cells that aren't NA (masked out by NLCD or missing data)
valid_cells <- complete.cases(pred_df)
if(sum(valid_cells) == 0) stop("No valid cells found for prediction. Check NLCD mask.")
# Build scenario data
pred_df$Cogongrass_Cover <- cogongrass_value
for (mem in mem_vars) pred_df[[mem]] <- 0 # Methods: MEMs held constant at 0
# Reorder to match model requirements
all_vars <- attr(rf_model$terms, "term.labels")
pred_df_valid <- pred_df[valid_cells, all_vars, drop = FALSE]
# Generate predictions
predicted_vals <- predict(rf_model, newdata = pred_df_valid)
# Map values back to raster space
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, mem_vars)
}## 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 Statisitics 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),
Max_Loss = min(v),
Max_Gain = max(v)
)
}) %>% bind_rows()
print(stats_summary)## Scenario Mean_Change Max_Loss Max_Gain
## 1 Diff_1_vs_0 0.0000000000 0.00000000 0.000000000
## 2 Diff_2_vs_0 -0.0021540131 -0.01394934 0.002456030
## 3 Diff_3_vs_0 -0.0042723674 -0.01680358 0.001551865
## 4 Diff_4_vs_0 -0.0031893542 -0.01614842 0.003525222
## 5 Diff_5_vs_0 -0.0030533206 -0.01546306 0.003627290
## 6 Diff_6_vs_0 -0.0059809641 -0.01815757 0.003074285
## 7 Diff_7_vs_0 -0.0061097043 -0.01806224 0.002744493
## 8 Diff_8_vs_0 -0.0054600432 -0.01930045 0.003588835
## 9 Diff_9_vs_0 -0.0064533995 -0.02049284 0.003906713
## 10 Diff_10_vs_0 -0.0059528291 -0.02006224 0.004664293
## 11 Diff_11_vs_0 -0.0077585724 -0.02355668 0.004115986
## 12 Diff_12_vs_0 -0.0083093447 -0.02778413 0.005370331
## 13 Diff_13_vs_0 -0.0076678153 -0.02830676 0.005504541
## 14 Diff_14_vs_0 -0.0093519640 -0.03213204 0.005096504
## 15 Diff_15_vs_0 -0.0092119642 -0.03213204 0.005624360
## 16 Diff_16_vs_0 -0.0098959369 -0.03346820 0.005228866
## 17 Diff_17_vs_0 -0.0106298301 -0.03488670 0.004815168
## 18 Diff_18_vs_0 -0.0141680850 -0.04545048 0.006023090
## 19 Diff_19_vs_0 -0.0162319467 -0.05172267 0.005773319
## 20 Diff_20_vs_0 -0.0127924791 -0.05158535 0.015512752
## 21 Diff_21_vs_0 -0.0104323249 -0.05441639 0.018929875
## 22 Diff_22_vs_0 -0.0097424302 -0.05719260 0.020692801
## 23 Diff_23_vs_0 -0.0070556732 -0.05626770 0.025451575
## 24 Diff_24_vs_0 -0.0070456318 -0.05652031 0.025451575
## 25 Diff_25_vs_0 -0.0064397541 -0.05567266 0.025767808
## 26 Diff_26_vs_0 -0.0072063180 -0.05754866 0.026291933
## 27 Diff_27_vs_0 -0.0082399993 -0.06340502 0.026581197
## 28 Diff_28_vs_0 -0.0062874268 -0.06016956 0.027171052
## 29 Diff_29_vs_0 -0.0049062673 -0.05961793 0.027624910
## 30 Diff_30_vs_0 -0.0051644401 -0.06274534 0.027874919
## 31 Diff_31_vs_0 -0.0054857286 -0.06274534 0.027316503
## 32 Diff_32_vs_0 0.0011483770 -0.06088260 0.032315468
## 33 Diff_33_vs_0 -0.0004380289 -0.06695150 0.032896355
## 34 Diff_34_vs_0 0.0011782388 -0.06624781 0.035408420
## 35 Diff_35_vs_0 0.0072456335 -0.06665967 0.042139969
## 36 Diff_36_vs_0 0.0083568628 -0.06600022 0.043397560
## 37 Diff_37_vs_0 0.0124433856 -0.06900259 0.057737536
## 38 Diff_38_vs_0 0.0137675130 -0.06900259 0.058829244
## 39 Diff_39_vs_0 0.0134999143 -0.06908232 0.058482580
## 40 Diff_40_vs_0 0.0144412892 -0.06908232 0.059403210
## 41 Diff_41_vs_0 0.0184937047 -0.06623721 0.069487902
## 42 Diff_42_vs_0 0.0181115545 -0.06821772 0.071552340
## 43 Diff_43_vs_0 0.0187098493 -0.06821772 0.072084156
## 44 Diff_44_vs_0 0.0183470254 -0.06786004 0.070371292
## 45 Diff_45_vs_0 0.0190174970 -0.06876698 0.072117867
## 46 Diff_46_vs_0 0.0200931985 -0.06655336 0.071392016
## 47 Diff_47_vs_0 0.0206258782 -0.06655336 0.071294718
## 48 Diff_48_vs_0 0.0191064566 -0.06757407 0.070575105
## 49 Diff_49_vs_0 0.0191852137 -0.06743069 0.070976002
## 50 Diff_50_vs_0 0.0189447628 -0.06733882 0.068259534
## 51 Diff_51_vs_0 0.0186182779 -0.06733882 0.067487132
## 52 Diff_52_vs_0 0.0186360666 -0.06733882 0.067757581
## 53 Diff_53_vs_0 0.0186360666 -0.06733882 0.067757581
## 54 Diff_54_vs_0 0.0188030811 -0.06733882 0.069717367
## 55 Diff_55_vs_0 0.0146487042 -0.07372641 0.060543340
## 56 Diff_56_vs_0 0.0134787738 -0.07372641 0.058959844
## 57 Diff_57_vs_0 0.0080625522 -0.08382356 0.057403494
## 58 Diff_58_vs_0 0.0009203502 -0.08556228 0.049371938
## 59 Diff_59_vs_0 -0.0081000160 -0.10097455 0.046659279
## 60 Diff_60_vs_0 -0.0092560231 -0.10153448 0.045956722
## 61 Diff_61_vs_0 -0.0084564573 -0.10032120 0.046224504
## 62 Diff_62_vs_0 -0.0085034323 -0.10100991 0.046224504
## 63 Diff_63_vs_0 -0.0121587579 -0.10745498 0.046224504
## 64 Diff_64_vs_0 -0.0227089131 -0.12904099 0.046224504
## 65 Diff_65_vs_0 -0.0236224508 -0.12904099 0.045465173
## 66 Diff_66_vs_0 -0.0250066436 -0.12904099 0.042640401
## 67 Diff_67_vs_0 -0.0270765664 -0.13099822 0.041228519
## 68 Diff_68_vs_0 -0.0287132448 -0.13589608 0.039720320
## 69 Diff_69_vs_0 -0.0303893284 -0.13892651 0.040668821
## 70 Diff_70_vs_0 -0.0460456701 -0.18011276 0.045188279
## 71 Diff_71_vs_0 -0.0464116649 -0.18011276 0.045188279
## 72 Diff_72_vs_0 -0.0466464828 -0.18082731 0.045369487
## 73 Diff_73_vs_0 -0.0444535865 -0.18082731 0.050197245
## 74 Diff_74_vs_0 -0.0447191905 -0.18082731 0.049392118
## 75 Diff_75_vs_0 -0.0447191905 -0.18082731 0.049392118
## 76 Diff_76_vs_0 -0.0447191905 -0.18082731 0.049392118
## 77 Diff_77_vs_0 -0.0465898367 -0.18126790 0.045148113
## 78 Diff_78_vs_0 -0.0465898367 -0.18126790 0.045148113
## 79 Diff_79_vs_0 -0.0591954274 -0.19946714 0.036588753
## 80 Diff_80_vs_0 -0.0775670203 -0.21280704 0.018447754
## 81 Diff_81_vs_0 -0.0878868200 -0.22314340 0.009208218
## 82 Diff_82_vs_0 -0.0887865590 -0.22396412 0.009208218
## 83 Diff_83_vs_0 -0.2995711595 -0.45656224 -0.171893089
## 84 Diff_84_vs_0 -0.3577491988 -0.50707599 -0.218858391
## 85 Diff_85_vs_0 -0.3577936592 -0.50707599 -0.218858391
## 86 Diff_86_vs_0 -0.3577936592 -0.50707599 -0.218858391
## 87 Diff_87_vs_0 -0.6922712985 -0.91755522 -0.461275359
## 88 Diff_88_vs_0 -0.7077630935 -0.95183815 -0.463285660
## 89 Diff_89_vs_0 -0.7077630935 -0.95183815 -0.463285660
## 90 Diff_90_vs_0 -0.7027803186 -0.94475127 -0.464129708
## 91 Diff_91_vs_0 -0.6991085152 -0.94198597 -0.459767567
## 92 Diff_92_vs_0 -0.6965744487 -0.93950273 -0.456285078
## 93 Diff_93_vs_0 -0.6965744487 -0.93950273 -0.456285078
## 94 Diff_94_vs_0 -0.6965744487 -0.93950273 -0.456285078
## 95 Diff_95_vs_0 -0.6965744487 -0.93950273 -0.456285078
## 96 Diff_96_vs_0 -0.6965744487 -0.93950273 -0.456285078
## 97 Diff_97_vs_0 -0.6965744487 -0.93950273 -0.456285078
## 98 Diff_98_vs_0 -0.6965744487 -0.93950273 -0.456285078
## 99 Diff_99_vs_0 -0.6965744487 -0.93950273 -0.456285078
## 100 Diff_100_vs_0 -0.6965744487 -0.93950273 -0.456285078Plot
plot_raster <- function(rast_layer, title, subtitle = NULL,
palette = "viridis", midpoint = NULL,
low = NULL, high = NULL, mid = NULL,
limits = NULL) {
# Convert SpatRaster to data frame for ggplot
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)) +
# Ensure states_sf is already loaded in your environment
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")
# To see it:
print(p_final)## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
Additional Analysis
## Area-Based Impact
# Example for one scenario
v <- values(diff_list[["Diff_100_vs_0"]])
prop_loss <- mean(v < 0, na.rm = TRUE)
prop_gain <- mean(v > 0, na.rm = TRUE)
## Threshold Mapping
severe_loss <- diff_list[["Diff_100_vs_0"]] < -0.7
plot(severe_loss, main = "Areas with Severe Predicted Loss (>0.7 decrease)")
# Regional Summaries
states_vect <- terra::vect(states_sf)
regional_summary <- lapply(names(diff_list), function(nm) {
r <- diff_list[[nm]]
vals <- terra::extract(r, states_vect, fun = mean, na.rm = TRUE)
# Add state names
vals$State <- states_sf$ID[vals$ID]
vals$Scenario <- nm
vals
}) %>% bind_rows()## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the raster
## Warning: [extract] transforming vector data to the CRS of the rasterregional_summary$Scenario_num <- as.numeric(gsub("Diff_|_vs_0", "", regional_summary$Scenario))
regional_summary$Scenario_num <- as.numeric(gsub("Diff_|_vs_0", "", regional_summary$Scenario))
regional_long <- regional_summary %>%
pivot_longer(
cols = starts_with("Shannon_"),
names_to = "Replicate",
values_to = "Mean_Change"
) %>%
filter(!is.na(Mean_Change))
print(regional_long)## # A tibble: 1,000 × 6
## ID State Scenario Scenario_num Replicate Mean_Change
## <int> <chr> <chr> <dbl> <chr> <dbl>
## 1 1 alabama Diff_1_vs_0 1 Shannon_1 0
## 2 2 arkansas Diff_1_vs_0 1 Shannon_1 0
## 3 3 florida Diff_1_vs_0 1 Shannon_1 0
## 4 4 georgia Diff_1_vs_0 1 Shannon_1 0
## 5 5 louisiana Diff_1_vs_0 1 Shannon_1 0
## 6 6 mississippi Diff_1_vs_0 1 Shannon_1 0
## 7 7 north carolina Diff_1_vs_0 1 Shannon_1 0
## 8 8 south carolina Diff_1_vs_0 1 Shannon_1 0
## 9 9 tennessee Diff_1_vs_0 1 Shannon_1 0
## 10 10 virginia Diff_1_vs_0 1 Shannon_1 0
## # ℹ 990 more rowsggplot(regional_long, aes(x = Scenario_num, y = Mean_Change, color = State)) +
geom_line() +
theme_bw() +
labs(x = "Cogongrass Cover (%)",
y = "Δ Shannon Diversity",
color = "State")
## Subset regional_long for 100% scenario
regional_100 <- regional_long %>%
filter(Scenario_num == 100)