Reworked Biodiversity Prediction
Alan Ivory
2026-06-25
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
Start Here
Aside from the code chunk “load_data” there is nothing else important above this point. It is how I collected nlcd, elevation etc. rasters, but those are all saved and ready to go. The code below is the actual workflow for building the env stack, running predictions, and analyzing results.
# 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.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
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(!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...
## |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= 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. Check NLCD mask.")
# 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 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.0001356290 -0.007701703 0.01886696
## 2 Diff_2_vs_0 0.0012985415 -0.030605133 0.02513047
## 3 Diff_3_vs_0 0.0003464533 -0.031955825 0.02569142
## 4 Diff_4_vs_0 0.0001623097 -0.034865115 0.02578952
## 5 Diff_5_vs_0 0.0009469812 -0.037114557 0.03552967
## 6 Diff_6_vs_0 0.0028589211 -0.048444702 0.03943584
## 7 Diff_7_vs_0 0.0030448051 -0.049212279 0.06042167
## 8 Diff_8_vs_0 0.0044482088 -0.049077183 0.07059618
## 9 Diff_9_vs_0 0.0031520803 -0.048343908 0.09670917
## 10 Diff_10_vs_0 0.0026407888 -0.049525306 0.10060042
## 11 Diff_11_vs_0 0.0034059948 -0.068899976 0.11157174
## 12 Diff_12_vs_0 0.0036137828 -0.079794582 0.12430000
## 13 Diff_13_vs_0 0.0051067339 -0.080592944 0.12430000
## 14 Diff_14_vs_0 0.0033961005 -0.092562239 0.12678151
## 15 Diff_15_vs_0 0.0027228508 -0.094338721 0.13314818
## 16 Diff_16_vs_0 0.0015774240 -0.101627664 0.13343950
## 17 Diff_17_vs_0 0.0006745478 -0.105587769 0.13735877
## 18 Diff_18_vs_0 0.0005670232 -0.106557257 0.13672746
## 19 Diff_19_vs_0 -0.0016408148 -0.109777571 0.13704183
## 20 Diff_20_vs_0 0.0002359400 -0.107769950 0.13852169
## 21 Diff_21_vs_0 -0.0014878391 -0.191695098 0.12947139
## 22 Diff_22_vs_0 0.0006770293 -0.198325316 0.13970027
## 23 Diff_23_vs_0 0.0018363881 -0.200017303 0.13878561
## 24 Diff_24_vs_0 0.0051559487 -0.200443834 0.16608199
## 25 Diff_25_vs_0 0.0050485733 -0.211993547 0.17854369
## 26 Diff_26_vs_0 0.0046767242 -0.213418565 0.18167894
## 27 Diff_27_vs_0 0.0029505598 -0.218194000 0.18875150
## 28 Diff_28_vs_0 0.0037415075 -0.217153256 0.19500866
## 29 Diff_29_vs_0 0.0057939135 -0.222023531 0.19918652
## 30 Diff_30_vs_0 0.0061298274 -0.219707626 0.19997227
## 31 Diff_31_vs_0 0.0095591858 -0.218397092 0.19984448
## 32 Diff_32_vs_0 0.0217528721 -0.217011236 0.20775359
## 33 Diff_33_vs_0 0.0212440256 -0.219742078 0.20977435
## 34 Diff_34_vs_0 0.0216392578 -0.220650905 0.21390007
## 35 Diff_35_vs_0 0.0233031734 -0.220495020 0.21366785
## 36 Diff_36_vs_0 0.0205190716 -0.237198662 0.21116477
## 37 Diff_37_vs_0 0.0222684845 -0.236160410 0.21336591
## 38 Diff_38_vs_0 0.0231988647 -0.235448398 0.21809616
## 39 Diff_39_vs_0 0.0293404946 -0.229471021 0.23044641
## 40 Diff_40_vs_0 0.0443139795 -0.221677157 0.26103981
## 41 Diff_41_vs_0 0.0527579741 -0.217161538 0.27773679
## 42 Diff_42_vs_0 0.0544097079 -0.216552544 0.28416686
## 43 Diff_43_vs_0 0.0566125262 -0.216446702 0.28924028
## 44 Diff_44_vs_0 0.0568104537 -0.215304206 0.28908263
## 45 Diff_45_vs_0 0.0622864549 -0.215366797 0.31383556
## 46 Diff_46_vs_0 0.0703279264 -0.211153479 0.33169245
## 47 Diff_47_vs_0 0.0781536020 -0.204317608 0.36582370
## 48 Diff_48_vs_0 0.0780383861 -0.204091955 0.36582370
## 49 Diff_49_vs_0 0.0775760817 -0.204296854 0.36328152
## 50 Diff_50_vs_0 0.0789341054 -0.201530620 0.36463211
## 51 Diff_51_vs_0 0.0784113122 -0.202125930 0.36346601
## 52 Diff_52_vs_0 0.0772504178 -0.203994742 0.35888934
## 53 Diff_53_vs_0 0.0771715049 -0.206847094 0.35967520
## 54 Diff_54_vs_0 0.0774628452 -0.207028951 0.35983005
## 55 Diff_55_vs_0 0.0772931183 -0.210197121 0.35921672
## 56 Diff_56_vs_0 0.0768955041 -0.212063524 0.35921672
## 57 Diff_57_vs_0 0.0773741561 -0.210871863 0.36059436
## 58 Diff_58_vs_0 0.0774650346 -0.210871863 0.36003046
## 59 Diff_59_vs_0 0.0809516528 -0.210170471 0.36853038
## 60 Diff_60_vs_0 0.0809737673 -0.210170471 0.36853038
## 61 Diff_61_vs_0 0.0811028651 -0.211346589 0.37089878
## 62 Diff_62_vs_0 0.0790629287 -0.215286409 0.37189814
## 63 Diff_63_vs_0 0.0815258992 -0.212022934 0.38080909
## 64 Diff_64_vs_0 0.0813978087 -0.213006372 0.38230923
## 65 Diff_65_vs_0 0.0780655218 -0.210710996 0.38493751
## 66 Diff_66_vs_0 0.0776210128 -0.210710996 0.38252104
## 67 Diff_67_vs_0 0.0770663241 -0.210710996 0.38252104
## 68 Diff_68_vs_0 0.0730758048 -0.218580457 0.38252104
## 69 Diff_69_vs_0 0.0734828674 -0.218151694 0.38129775
## 70 Diff_70_vs_0 0.0623700319 -0.242969883 0.37986219
## 71 Diff_71_vs_0 0.0623479910 -0.242969883 0.37986219
## 72 Diff_72_vs_0 0.0585921356 -0.247482725 0.37905681
## 73 Diff_73_vs_0 0.0335432996 -0.302407189 0.37436766
## 74 Diff_74_vs_0 0.0162075082 -0.336114407 0.36887995
## 75 Diff_75_vs_0 0.0133512096 -0.340534597 0.36881159
## 76 Diff_76_vs_0 0.0128212363 -0.340242420 0.36479731
## 77 Diff_77_vs_0 0.0125261076 -0.340242420 0.36479731
## 78 Diff_78_vs_0 0.0116188403 -0.340242420 0.36194245
## 79 Diff_79_vs_0 0.0039604238 -0.352064209 0.34795824
## 80 Diff_80_vs_0 -0.0169450402 -0.375664229 0.31747742
## 81 Diff_81_vs_0 -0.0171830120 -0.377589096 0.31747742
## 82 Diff_82_vs_0 -0.0171830120 -0.377589096 0.31747742
## 83 Diff_83_vs_0 -0.0683956110 -0.419646346 0.26749144
## 84 Diff_84_vs_0 -0.0721275413 -0.422996055 0.26576706
## 85 Diff_85_vs_0 -0.0721275413 -0.422996055 0.26576706
## 86 Diff_86_vs_0 -0.0721275413 -0.422996055 0.26576706
## 87 Diff_87_vs_0 -0.4473657701 -0.852807458 -0.17817782
## 88 Diff_88_vs_0 -0.4759827990 -0.894780528 -0.19835943
## 89 Diff_89_vs_0 -0.4759827990 -0.894780528 -0.19835943
## 90 Diff_90_vs_0 -0.4831894253 -0.899419291 -0.20447579
## 91 Diff_91_vs_0 -0.4821759852 -0.898405851 -0.20346235
## 92 Diff_92_vs_0 -0.4682266067 -0.884329236 -0.18938573
## 93 Diff_93_vs_0 -0.4682266067 -0.884329236 -0.18938573
## 94 Diff_94_vs_0 -0.4682266067 -0.884329236 -0.18938573
## 95 Diff_95_vs_0 -0.4682266067 -0.884329236 -0.18938573
## 96 Diff_96_vs_0 -0.4682266067 -0.884329236 -0.18938573
## 97 Diff_97_vs_0 -0.4682266067 -0.884329236 -0.18938573
## 98 Diff_98_vs_0 -0.4682266067 -0.884329236 -0.18938573
## 99 Diff_99_vs_0 -0.4682266067 -0.884329236 -0.18938573
## 100 Diff_100_vs_0 -0.4682266067 -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")
# 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.00171
## 2 2 arkansas Diff_1_vs_0 1 Shannon_1 0.000649
## 3 3 florida Diff_1_vs_0 1 Shannon_1 -0.000134
## 4 4 georgia Diff_1_vs_0 1 Shannon_1 0.000930
## 5 5 louisiana Diff_1_vs_0 1 Shannon_1 0.00130
## 6 6 mississippi Diff_1_vs_0 1 Shannon_1 0.00157
## 7 7 north carolina Diff_1_vs_0 1 Shannon_1 0.000228
## 8 8 south carolina Diff_1_vs_0 1 Shannon_1 0.000499
## 9 9 tennessee Diff_1_vs_0 1 Shannon_1 0.000317
## 10 10 virginia Diff_1_vs_0 1 Shannon_1 -0.000208
## # ℹ 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)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 ONCE using the first raster
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))
)
}
# Check your column names first:
# 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 = "#d73027",
mid = "white",
high = "#1a9850",
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 = "#d73027",
mid = "white",
high = "#1a9850",
midpoint = 0
) +
theme_bw() +
labs(
title = "Δ Shannon Diversity (Level 2 Ecoregions, 100% Cogongrass)",
fill = "Change"
)