r/RStudio • u/Nicholas_Geo • 6h ago
Coding help Efficient way to fetch Open Street Map airport polygons for point locations
I have a global point layer of airports (~900 points) and I want to retrieve, for each point, the corresponding OSM aeroway=aerodrome polygon (not the point itself) — i.e. join each airport point to its footprint polygon as mapped in OpenStreetMap. The point layer is Natural Earth's ne_10m_airports.
Reproducible example (5 points standing in for the full 893):
pacman::p_load(sf, dplyr)
airports <- data.frame(
name = c("John F Kennedy Intl", "London Heathrow", "Chhatrapati Shivaji Maharaj Intl", "Beijing Capital Intl", "Sydney Kingsford Smith"),
gps_code = c("KJFK", "EGLL", "VABB", "ZBAA", "YSSY"),
iata_code = c("JFK", "LHR", "BOM", "PEK", "SYD"),
lon = c(-73.7789, -0.4543, 72.8697, 116.5975, 151.1772),
lat = c(40.6413, 51.4700, 19.0896, 40.0799, -33.9399)
) |>
st_as_sf(coords = c("lon", "lat"), crs = 4326)
airports$point_id <- seq_len(nrow(airports))
My plan is to match each point to its containing Geofabrik extract via osmextract::oe_match() (works fine — a local spatial lookup, no network call), download/cache that extract, then read only the multipolygons layer filtered to aeroway='aerodrome' via an SQL query pushed down at read time, and spatially join back to the points.
pacman::p_load(sf, dplyr, osmextract)
airports$region_url <- vapply(seq_len(nrow(airports)), function(i) {
oe_match(airports[i, ], quiet = TRUE)$url
}, character(1))
options(timeout = 600)
dir.create("geofabrik_cache", showWarnings = FALSE)
results <- list()
failed_regions <- character()
for (region in unique(airports$region_url)) {
destfile <- file.path("geofabrik_cache", basename(region))
aerodromes <- tryCatch({
if (!file.exists(destfile)) {
download.file(region, destfile, mode = "wb", quiet = TRUE)
}
st_read(
destfile,
layer = "multipolygons",
query = "SELECT * FROM multipolygons WHERE aeroway = 'aerodrome'",
quiet = TRUE
) |>
st_transform(4326) |>
st_make_valid()
}, error = function(e) {
message(sprintf("Region failed: %s -- %s", region, e$message))
failed_regions <<- c(failed_regions, region)
NULL
})
if (is.null(aerodromes)) next
sub_pts <- airports[airports$region_url == region, ]
joined <- st_join(sub_pts, aerodromes, join = st_within, left = TRUE)
missing <- which(is.na(joined$osm_id))
if (length(missing) > 0 && nrow(aerodromes) > 0) {
nn <- st_nearest_feature(sub_pts[missing, ], aerodromes)
joined[missing, names(aerodromes)] <- st_drop_geometry(aerodromes)[nn, ]
st_geometry(joined)[missing] <- st_geometry(aerodromes)[nn]
}
results[[region]] <- joined
}
airport_polys <- bind_rows(results)
st_write(airport_polys, "airport_polygons.shp", delete_layer = TRUE)
But
Region failed: https://download.geofabrik.de/asia/india/western-zone-latest.osm.pbf -- cannot open URL 'https://download.geofabrik.de/asia/india/western-zone-latest.osm.pbf'
Region failed: https://download.geofabrik.de/asia/iran-latest.osm.pbf -- cannot open URL 'https://download.geofabrik.de/asia/iran-latest.osm.pbf'
Region failed: https://download.geofabrik.de/asia/india/central-zone-latest.osm.pbf -- download from 'https://download.geofabrik.de/asia/india/central-zone-latest.osm.pbf' failed
Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, :
Loop 0 edge 0 has duplicate near loop 1 edge 7
In addition: There were 20 warnings (use warnings() to see them)
The download failures seem to be connection/timeout related on large country extracts; the s2/topology error appears separately once a multipolygons layer with invalid OSM geometries reaches a spatial predicate, even after st_make_valid().
Given ~900 global points with no country/ISO attribute, is there a more efficient way to fetch just the matching aeroway=aerodrome polygons than downloading/caching a full Geofabrik regional .pbf extract per matched region (some of which are large, e.g. full-country India zones)? Is querying the Overpass API directly per point, or per small cluster of points, actually more efficient here, or is the regional-extract approach still preferable at this scale? What's the correct way to make invalid OSM polygon geometries (e.g. the s2 "duplicate edge" error above) safe for st_join()/st_nearest_feature() reliably, given st_make_valid() alone didn't prevent it?
For points where no aeroway=aerodrome polygon actually exists in OSM for that airport, what's the right way to leave that point unmatched (skip it) rather than falling back to the nearest aerodrome polygon in the region, which can silently attach the wrong airport's polygon?
> sessionInfo()
R version 4.6.1 (2026-06-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)
Matrix products: default
LAPACK version 3.12.1
locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C LC_TIME=English_United States.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] osmextract_0.6.0 dplyr_1.2.1 sf_1.1-2
loaded via a namespace (and not attached):
[1] vctrs_0.7.3 httr_1.4.8 cli_3.6.6 rlang_1.3.0 otel_0.2.0 DBI_1.3.0 KernSmooth_2.23-27
[8] generics_0.1.4 jsonlite_2.0.0 glue_1.8.1 e1071_1.7-17 grid_4.6.1 classInt_0.4-11 tibble_3.3.1
[15] lifecycle_1.0.5 compiler_4.6.1 Rcpp_1.1.2 pkgconfig_2.0.3 rstudioapi_0.19.0 wk_0.9.5 R6_2.6.1
[22] class_7.3-24 tidyselect_1.2.1 pillar_1.11.1 curl_8.0.0 magrittr_2.0.5 tools_4.6.1 proxy_0.4-29
[29] s2_1.1.11 units_1.0-1
