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

I’m fitting a Beta GAM in `mgcv` to model the proportion of daytime outdoor-use time in laying hens. The response is the proportion of minutes spent outdoors during an 08:00–22:00 observation window, and the model contains repeated observations from individual hens.

My current model is approximately:

prop_outside \~
  coop_id +
  s(study_day, bs = "cr", k = 40) +
  s(max_temp, bs = "cr", k = 8) +
  s(mean_dew_point, bs = "cr", k = 7) +
  rain_any +
  s(log_total_rain, by = rain_status, bs = "cr", k = 7) +
  s(eid, bs = "re")

using:

family = betar(link = "logit")
method = "REML"

The dataset contains roughly 8,800 positive hen-day observations from 124 hens across about 137 study dates. The model explains about 60% of the deviance.

The issue is the `study_day` smooth. On a same-row model comparison, I get:

s(study_day)
k'       = 39
edf      = 34.58
k-index  = 0.949
p-value  < 0.001

I also tried an alternative weather specification using mean temperature, mean humidity and mean wind instead of maximum temperature and dew point:

prop_outside \~
  coop_id +
  s(study_day, bs = "cr", k = 40) +
  s(mean_temp, bs = "cr", k = 8) +
  s(mean_humidity, bs = "cr", k = 8) +
  s(mean_wind, bs = "cr", k = 8) +
  rain_any +
  s(log_total_rain, by = rain_status, bs = "cr", k = 7) +
  s(eid, bs = "re")

The same issue persists:

s(study_day)
k'       = 39
edf      = 34.84
k-index  = 0.966
p-value  = 0.0175

The study-day smooth is visually quite wiggly, particularly early in the study, and the EDF is already fairly close to the available basis dimension. Changing the weather specification does not materially improve overall fit: both models explain about 60% of deviance and have essentially identical AIC.

There is also substantial temporal/weather dependence. For example, in the alternative model, observed concurvity is about 0.90 for mean temperature and 0.48 for study day.

My questions are:

1. Is the low `k-index` plus EDF ≈ 35/39 sufficient reason to increase `k` for `study_day`, for example from 40 to 60?
2. If increasing `k` gives essentially the same fitted curve/predictions but the k-check remains significant, would you consider the current smooth adequate?
3. Could this be indicating residual temporal autocorrelation rather than simply an insufficient basis dimension?
4. Would you model study date differently in this setting—for example with an autocorrelation structure, a different smooth basis, or another temporal term?
5. Since weather variables themselves follow study date seasonally, how would you distinguish genuine temporal structure from weather-related temporal confounding?

My current plan is to compare `k = 40` and `k = 60` on the same observations and assess whether fitted values, the study-day effect, held-out-date prediction, and conclusions materially change, rather than selecting the model based only on the k-check p-value.

I’d appreciate advice on how to interpret this persistent `study_day` diagnostic and what additional diagnostic/model comparison would be most appropriate.