WSI Prediction Map Plot

# Load necessary libraries
library(terra)
library(sf)
library(ggplot2)
library(dplyr)
library(purrr)
library(USA.state.boundaries)

# Define file paths and target CRS
effect_dir <- "C:/Users/Todd/Desktop/Eastern WSI/WSI_Predicted_Effects_EasternOnly"
years <- 2004:2024
effect_files <- file.path(effect_dir, paste0("predicted_effect_", years, ".tif"))
target_crs <- "EPSG:5070"  # Albers Equal Area (NAD83)

# Define crop box in WGS84 and transform to target CRS
crop_box_ll <- st_as_sfc(st_bbox(c(xmin = -102, xmax = -66, ymin = 39, ymax = 52), crs = 4326))
crop_box_proj <- st_transform(crop_box_ll, crs = target_crs)

# Load and process state boundaries
usa_states <- st_transform(st_make_valid(state_boundaries_wgs84), crs = target_crs)
usa_cropped <- st_crop(usa_states, crop_box_proj)

effect_colors <- c(
  "Strong Negative"   = "#67001f",  # Deep Maroon
  "Moderate Negative" = "#b2182b",  # Dark Red
  "Mild Negative"     = "#ef8a62",  # Light Red
  "Marginal Negative" = "#fddbc7",  # Pale Pink
  "Nuetral/Average"   = "gray75",  # Light Gray
  "Marginal Positive" = "#d9f0d3",  # Pale Green
  "Mild Positive"     = "#a6dba0",  # Light Green
  "Moderate Positive" = "#4dac26",  # Green
  "Strong Positive"   = "#00441b"   # Dark Green
)

# Process raster files
effect_df_list <- map2(effect_files, years, ~ {
  r <- rast(.x)
  if (is.null(crs(r))) {
    crs(r) <- "EPSG:4326"  # Assign CRS if missing
  }
  r_proj <- project(r, target_crs)
  r_crop <- crop(r_proj, crop_box_proj)
  r_down <- aggregate(r_crop, fact = 2)
  if (all(is.na(values(r_down)))) return(NULL)
  df <- as.data.frame(r_down, xy = TRUE, na.rm = TRUE)
  names(df)[3] <- "effect"
  df$year <- .y
  df$class <- cut(df$effect,
    breaks = c(-Inf, -0.20, -0.1, -0.05, -0.02, 0.02, 0.05, 0.1, 0.20, Inf),
    labels = c(
      "Strong Negative", "Moderate Negative", "Mild Negative", "Marginal Negative",
      "Nuetral/Average", "Marginal Positive", "Mild Positive", "Moderate Positive", "Strong Positive"
    ),
    include.lowest = TRUE
  )
  df
})

# Combine all valid years into a single data frame
effect_all <- bind_rows(compact(effect_df_list))

# Plot the results
ggplot() +
  geom_tile(data = effect_all, aes(x = x, y = y, fill = class)) +
  geom_sf(data = usa_cropped, fill = NA, color = "black", linewidth = 0.3) +
  scale_fill_manual(values = effect_colors, name = "Effect Class", drop = FALSE) +
  coord_sf(crs = target_crs, expand = FALSE) +
  facet_wrap(~year, ncol = 1) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "top",           # Position the legend at the top
    legend.direction = "vertical",     # Arrange legend items vertically
    legend.box = "vertical",           # Stack multiple legends vertically if present
    legend.title.align = 0.5,          # Center-align the legend title
    strip.text = element_text(face = "bold"),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank()
  ) +
  labs(
    title = "Predicted WSI Effects on Deer Harvest (2004–2024)",
    x = NULL, y = NULL
  )