library(sf)
library(FNN)
library(dplyr)

# -------- Read inputs from the current working directory --------
footprints    <- st_read("Footprints_county_sample.shp", quiet = TRUE)
transit_stops <- st_read("National_Transit_Map_Stops.shp", quiet = TRUE)

# Align CRS: start by matching stops to footprints
transit_stops <- st_transform(transit_stops, st_crs(footprints))

# If CRS is geographic (degrees), reproject both to a projected CRS (meters)
# EPSG:5070 (NAD83 / Conus Albers) is a good CONUS-wide choice.
if (sf::st_is_longlat(footprints)) {
  target_crs   <- 5070
  footprints    <- st_transform(footprints, target_crs)
  transit_stops <- st_transform(transit_stops, target_crs)
}

# Ensure an 'index' column exists
if (!"index" %in% names(footprints)) {
  footprints$index <- seq_len(nrow(footprints))
}

# Centroid coordinates (matrices with x,y)
foot_xy   <- st_coordinates(st_centroid(footprints))
stop_xy   <- st_coordinates(st_centroid(transit_stops))

# 1-NN distance from each footprint to the nearest transit stop (meters if projected)
nn <- FNN::get.knnx(data = stop_xy, query = foot_xy, k = 1)

# Build the output table: index + distance
Foot_transit_distances <- tibble::tibble(
  index = footprints$index,
  Foot_transit_distances = as.numeric(nn$nn.dist)
)

# Save result
save(Foot_transit_distances, file = "Foot_transit_distances.RData")