library(sf)
library(dplyr)
library(osmextract)


setwd("")

# -------- Inputs (files expected in your working directory) --------
pbf_files <- c(
  "us-south-latest.osm.pbf",
  "us-west-latest.osm.pbf",
  "us-northeast-latest.osm.pbf",
  "us-midwest-latest.osm.pbf"
)

wanted_hw <- c("residential","motorway","trunk","primary",
               "secondary","tertiary","unclassified","living_street")  # fixed 'tertiary'

sql <- paste0(
  "SELECT geometry, highway FROM 'lines' WHERE highway IN ('",
  paste(wanted_hw, collapse = "', '"),
  "')"
)

# -------- Read & combine regional roads --------
road_list <- vector("list", length(pbf_files))
for (i in seq_along(pbf_files)) {
  message("Reading: ", pbf_files[i])
  road_list[[i]] <- oe_read(pbf_files[i], query = sql)
}
roads <- do.call(rbind, road_list)

# -------- Read footprints and find nearest road + distance --------
selection <- st_read("Footprints_county_sample.shp", quiet = TRUE) |>
  st_centroid()

# Align CRS
roads <- st_transform(roads, st_crs(selection))

# Indices of nearest road feature per footprint
nearest_idx <- st_nearest_feature(selection, roads)

# Join nearest road attributes (highway class) to each footprint
c <- cbind(selection, roads[nearest_idx, "highway"])

# Distance from each footprint centroid to its nearest road (same CRS units)
distances <- st_distance(selection, roads[nearest_idx, ], by_element = TRUE)

c$Road_dist <- as.numeric(distances)

# Keep only attributes; save result
c <- st_drop_geometry(c)
save(c, file = "Road_type_distance.RData")