library(sf)
library(raster)

setwd("")


# --- Inputs (files should be in your working directory) ---
# Vector footprints
Selection <- st_read("Footprints_county_sample.shp", quiet = TRUE)

# Land cover raster (NLCD 2019 CONUS IMG)
nlcd <- raster("nlcd_2019_land_cover_l48_20210604.img")

# --- 500 m buffer around centroids ---
Footprints_c      <- st_centroid(Selection)
Footprints_500_m  <- st_buffer(Footprints_c, 500)
Footprints_500_m$buffer <- 500

# Ensure an 'index' id exists (create one if missing)
if (!"index" %in% names(Footprints_500_m)) {
  Footprints_500_m$index <- seq_len(nrow(Footprints_500_m))
}

# Match CRS to raster
Footprints_500_m <- st_transform(Footprints_500_m, crs(nlcd))

# --- Recode NLCD classes to binary masks ---
# Lookup tables (NLCD classes: 11,21,22,23,24,31,41,42,43,52,71,81,82,90,95)
lut_water <- data.frame(
  id = c(11,21,22,23,24,31,41,42,43,52,71,81,82,90,95),
  v  = c(1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0)
)

lut_green <- data.frame(
  id = c(11,21,22,23,24,31,41,42,43,52,71,81,82,90,95),
  v  = c(0,  0,  0,  0,  0,  0,  1,  1,  1,  1,  1,  1,  1,  1,  1)
)

lut_trees <- data.frame(
  id = c(11,21,22,23,24,31,41,42,43,52,71,81,82,90,95),
  v  = c(0,  0,  0,  0,  0,  0,  1,  1,  1,  0,  0,  0,  0,  0,  0)
)

Water <- subs(nlcd, lut_water)
Green <- subs(nlcd, lut_green)
Trees <- subs(nlcd, lut_trees)

# --- % cover within 500 m buffers (mean of binary mask * 100) ---
Water_Pct <- raster::extract(Water, Footprints_500_m, fun = mean, na.rm = TRUE)
Green_Pct <- raster::extract(Green, Footprints_500_m, fun = mean, na.rm = TRUE)
Trees_Pct <- raster::extract(Trees, Footprints_500_m, fun = mean, na.rm = TRUE)

Land_cover_500m <- data.frame(
  index     = Footprints_500_m$index,
  Water_Pct = 100 * Water_Pct,
  Green_Pct = 100 * Green_Pct,
  Trees_Pct = 100 * Trees_Pct
)

write.csv(Land_cover_500m, "Land_cover_500m.csv", row.names = FALSE)