Reverse Geocoding Points to Administrative Boundaries

Most reverse-geocoding jobs are not asking for a street address at all — they are asking which district, ZIP code, census tract or municipality a coordinate falls in, and that question is answered by a polygon file you already have rather than by an HTTP call you pay for. This guide is for anyone labelling a table of coordinates with administrative units at scale: it sits under Geocoding & Address Data Pipelines in Geospatial Data Ingestion & Processing Workflows, and takes the local-join path that the parent topic points to once volumes outgrow a provider.

Why This Approach / What Goes Wrong

A provider round-trip costs one network call per point. At a realistic 50–200 ms of latency that is roughly ten hours of wall clock for 250,000 points before any rate limit applies, and the public Nominatim policy caps you at one request per second, which turns the same job into days. The local join reads one boundary file, builds an R-tree over a few thousand polygons, and answers every point in microseconds. The economics are not close, and they are the least interesting part of the argument.

The decisive property is determinism. A provider's answer is a function of its database on the day you asked: boundaries get re-cut, a municipality merges, an internal version ships, and the same input silently returns a different district six months later. Nothing in your code changed, but your yearly report no longer reproduces. A pinned boundary file — municipalities_2024-01.gpkg, checksummed and stored next to the code — turns the labelling step into a pure function. It runs on an air-gapped machine, in CI, inside a container with no egress, and produces byte-identical output every time. That is the same reasoning that makes a cached geocode preferable to a live one in batch geocoding with GeoPy, taken one step further: here there is no provider to cache.

What goes wrong locally is a short, well-understood list. The two layers disagree on their coordinate reference system and every predicate quietly evaluates false. The boundary file carries several administrative levels stacked in one table, so each point matches three overlapping polygons and the row count triples. A point sits exactly on a shared border, or in the sliver of sea between a coastline polygon and a harbour, and matches nothing. None of these raise; all of them corrupt the output. The rest of this guide is a join that handles each case explicitly instead of hoping.

Prerequisites

• geopandas>=1.0 — sjoin, sjoin_nearest, and pyogrio as the default I/O engine
• shapely>=2.0 — the vectorized GEOS predicates behind within
• pyogrio>=0.8 — column pushdown when reading the boundary file
• pyproj>=3.6 — CRS objects and estimate_utm_crs()

conda install -c conda-forge "geopandas=1.0.*" "shapely=2.0.*" "pyogrio=0.8.*" "pyproj=3.6.*"

You also need boundary polygons you control. National statistical agencies publish them as GeoPackage or Shapefile; GADM and Natural Earth cover the world at coarser resolution. Whatever the source, download it once, record the vintage in the filename, and store it as an artefact — a URL you re-fetch at runtime is a provider dependency wearing a different hat.

Step-by-Step Implementation

1. Load the boundary layer once, with only the columns you need. Passing columns= pushes the projection down to OGR through pyogrio, so attribute tables with sixty fields never reach memory. Reading directly to GeoParquet is faster still if you control the format.

import geopandas as gpd

boundaries = gpd.read_file(
    "data/admin/municipalities_2024-01.gpkg",
    layer="municipalities",
    columns=["muni_code", "muni_name", "admin_level"],
)

# One level only. Stacked levels in a single table are the usual cause of
# a join that returns three rows per input point.
boundaries = boundaries[boundaries["admin_level"] == 8].copy()

# Deterministic tie-break key, computed in an equal-area CRS — EPSG:3035 over
# Europe, EPSG:5070 over the conterminous US. Never a single UTM zone nationally.
boundaries["unit_area_m2"] = boundaries.to_crs("EPSG:3035").area

print(boundaries.crs, len(boundaries))      # EPSG:4326 8131

2. Build the point layer and align the CRS explicitly. points_from_xy takes longitude first, the opposite order from the (lat, lon) tuples a geocoder hands back — the ordering trap covered in the parent topic. sjoin only warns on a CRS mismatch and then computes an all-NaN result anyway, which is easy to lose in a log, so align before you call it.

import pandas as pd

incidents = pd.read_csv("data/incidents_2024.csv")     # incident_id, lon, lat

points = gpd.GeoDataFrame(
    incidents,
    geometry=gpd.points_from_xy(incidents["lon"], incidents["lat"]),
    crs="EPSG:4326",                       # declare it; never let it default to None
)

if boundaries.crs != points.crs:
    boundaries = boundaries.to_crs(points.crs)

Containment is projection-independent — a point either falls inside a polygon or it does not, and reprojecting both layers consistently cannot change the answer. So the containment join is legitimate in EPSG:4326. The moment a distance enters, in step 5, that stops being true and a metric CRS becomes mandatory.

3. Run the join with predicate="within", not the default. The default predicate is intersects, which counts boundary contact as a match and duplicates any point sitting on a shared edge into both neighbouring units. within demands strict interior containment. The how="left" mechanics — every left row survives, misses become NaN — are worked through in performing left joins with GeoPandas sjoin.

joined = gpd.sjoin(
    points,
    boundaries[["muni_code", "muni_name", "unit_area_m2", "geometry"]],
    how="left",
    predicate="within",
)
joined["match_kind"] = "within"

4. Collapse multi-matches with a deterministic tie-break. A point can still match two polygons when units genuinely overlap — enclaves, disputed strips, unclean data with slivers. GeoPandas signals this by repeating the left index, never by raising, so detect it before anything downstream counts rows. Sorting by unit area picks the most specific unit; the code column breaks remaining ties so the result never depends on row order in the source file.

multi_ids = joined.index[joined.index.duplicated(keep=False)].unique()
print(f"{len(multi_ids)} points matched more than one polygon")

joined = joined.sort_values(["unit_area_m2", "muni_code"], kind="stable")
joined = joined[~joined.index.duplicated(keep="first")].sort_index()
joined.loc[multi_ids, "match_kind"] = "tiebreak"

5. Rescue the unmatched with a nearest-boundary fallback in a projected CRS. Points on an exact border, coastal points a few metres offshore, and coordinates rounded to four decimal places all land outside every polygon. sjoin_nearest assigns the closest unit within a tolerance you choose — and its max_distance is expressed in CRS units, so a geographic CRS makes 250 mean 250 degrees. Reproject both sides to a metric CRS first; estimate_utm_crs() derives the right zone from the data extent, as in choosing a UTM zone automatically. GeoPandas emits a warning if you forget, but the numbers still come back.

metric_crs = points.estimate_utm_crs()          # e.g. EPSG:32633
missing = joined["muni_code"].isna()
print(f"{int(missing.sum())} points matched no polygon")

if missing.any():
    orphans = points.loc[missing].to_crs(metric_crs)
    nearby = boundaries.to_crs(metric_crs)[["muni_code", "muni_name", "geometry"]]

    rescued = gpd.sjoin_nearest(
        orphans, nearby, how="left", max_distance=250, distance_col="dist_m"
    )
    # Equidistant polygons duplicate the row here too — same rule, keep one.
    rescued = rescued[~rescued.index.duplicated(keep="first")]
    rescued = rescued[rescued["muni_code"].notna()]

    joined.loc[rescued.index, "muni_code"] = rescued["muni_code"]
    joined.loc[rescued.index, "muni_name"] = rescued["muni_name"]
    joined.loc[rescued.index, "match_kind"] = "nearest"

joined.loc[joined["muni_code"].isna(), "match_kind"] = "unmatched"

Choose the tolerance from the geometry of the problem, not from how many nulls you want gone. Two hundred and fifty metres recovers border and shoreline artefacts; two kilometres starts inventing answers. If your fallback needs ranked candidates rather than a single winner, the index structures in Nearest Neighbor & KD-Tree Search give you k results and their distances.

6. Keep the label column, not just the label. Downstream consumers need to know whether a district came from strict containment or from a 200 m guess. A choropleth of counts per district can happily use both; an audit trail cannot. Carrying match_kind costs one string column and removes every future argument about provenance. Write it out with the result — joined.to_parquet("incidents_labelled.parquet") keeps the geometry, the CRS and the label in one file — and make the unmatched subset its own artefact so somebody can actually look at it:

joined.to_parquet("out/incidents_labelled_2024.parquet")
joined.loc[joined["match_kind"] == "unmatched"].to_file("out/review_unmatched.gpkg")

On cost: the join itself is linear in the number of points once the R-tree exists, and the tree is built lazily on first access of boundaries.sindex. If you join the same boundary layer repeatedly in one process, hold a single boundaries frame and reuse it rather than re-reading the file, which is where most of the runtime actually goes. Tens of millions of points still fit comfortably on one machine; beyond that, partition the point layer and push the operation to parallel spatial joins with Dask-GeoPandas, where the boundary layer is small enough to broadcast to every worker.

Verification

The join is correct when the row count is unchanged, the index is unique, and the same inputs hash to the same output — the property that made the local path worth taking.

import hashlib

assert len(joined) == len(points), "row count changed — a multi-match survived dedup"
assert joined.index.is_unique, "duplicate index — tie-break did not run"
assert points.crs == boundaries.crs, "layers drifted apart in CRS"

print(joined["match_kind"].value_counts())
# within       48213
# tiebreak       117
# nearest         62
# unmatched        9
# Name: count, dtype: int64

# Reproducibility gate: pin this digest in a test and it will catch a swapped
# boundary file, a reordered source table, or a changed tie-break rule.
labels = joined.sort_index()["muni_code"].astype("string").fillna("").str.cat()
print(hashlib.sha256(labels.encode()).hexdigest()[:16])
# 9f3c1ab7e0d54c62

Edge Cases & Debugging

• Every row comes back NaN. The layers are in different CRSs; sjoin warns rather than raising, so check points.crs == boundaries.crs before the call and .to_crs() one side.
• The row count multiplied by three. The boundary file stacks administrative levels in one table. Filter to a single admin_level before joining, as in step 1 — the tie-break should be a rare correction, not the main path.
• A TopologyException, or a polygon that matches nothing. Self-intersecting rings make GEOS predicates unreliable. Run boundaries.geometry = boundaries.geometry.make_valid() at load time; the failure modes are catalogued in topology validation & repair.
• sjoin_nearest warns about a geographic CRS. max_distance was interpreted in degrees. Reproject both frames to a metric CRS — never EPSG:3857, whose scale error grows with latitude and makes the tolerance mean different distances in different rows.
• The unmatched count jumps after a data refresh. A new boundary vintage changed the geometry, the column names, or the code scheme. Diff the code sets — set(new["muni_code"]) - set(old["muni_code"]) — before trusting any comparison across vintages.
• Coastal points all land in the unmatched bucket. Administrative coastlines are generalised inland of the real shore. Either raise the fallback tolerance for maritime records only, or join against a version of the layer buffered seaward by a few hundred metres.

Frequently Asked Questions

When is the API still the right call? When you need something the polygons do not contain. A boundary file returns exactly the attributes in its table — a code, a name, maybe a population. If the answer has to include a street name or house number, only a geocoder can produce it. The two also compose well: resolve administrative units locally for all rows, then spend provider quota on the small subset that genuinely needs a street-level address. The trade-offs between hosted and self-run providers are covered in Nominatim vs Pelias for self-hosted geocoding.

At what volume does the local join start paying off? Almost immediately, but the crossover that matters is not volume — it is repetition. A one-off lookup of two hundred points is fine over HTTP. The moment the job runs on a schedule, or its output feeds a report someone will re-run next quarter, the pinned file wins regardless of size, because it is the only version that returns the same answer twice. Above roughly ten thousand points the runtime argument becomes decisive on its own.

Should I use within or contains? They are the same predicate read from opposite ends: points.within(polygon) and polygon.contains(point) evaluate the same relationship. In sjoin the predicate is applied left-to-right, so with points on the left, within is the one that expresses "each point is inside a unit". Passing contains with a point layer on the left asks whether each point contains a polygon, which is never true and yields an empty join.

How do I hold the boundary file so results stay reproducible? Treat it as a versioned artefact, not a download. Keep the vintage in the filename, store a SHA-256 of the file next to the code, and assert the digest at load time. When a new vintage arrives, bump it deliberately, re-run the whole history, and record that outputs changed — which is exactly the audit trail a provider API can never give you, since its database changes under you without notice.