⚙️ How-To Geocoding & addresses Joins & relationships Geocoding & Addresses · step 11 of 20 Updated

How to Reverse Geocode Points to Addresses and Areas in Python

Problem statement

You have coordinates and you need names: which country is this sensor in, which neighbourhood does this GPS fix belong to, what is the nearest street address to this photograph.

There are two completely different ways to answer that, and choosing the wrong one is the usual mistake:

  • Call a reverse geocoding API — one HTTP request per point, rate-limited, gives you an address string.
  • Spatially join to a boundary layer you already have — one in-process operation for the whole file, gives you the area's name and code.

For 200,000 GPS fixes, the API route is a policy violation that takes days. The spatial join runs in seconds. Measured on 13,464,017 points against 4,596 province polygons, an in-process join completed in 61.8 seconds and produced 12,942,217 matches.

The rule: if the answer you need is an area, join. Only call an API when you genuinely need a street address.

Quick answer

For areas — countries, regions, postcodes, neighbourhoods — join, do not call:

import geopandas as gpd

points = gpd.GeoDataFrame(
    df, geometry=gpd.points_from_xy(df.lon, df.lat), crs="EPSG:4326")
areas = gpd.read_file("admin2.gpkg")[["name", "code", "geometry"]]

labelled = gpd.sjoin(points, areas.to_crs(points.crs),
                     how="left", predicate="within")
print(f"{labelled['name'].isna().sum():,} points fell outside every area")

For a street address, call an API — once per point, throttled, cached:

def reverse(client, lat, lon, zoom=18):
    """zoom controls the level of detail: 18 = building, 10 = city, 3 = country."""
    r = client.session.get(f"{client.base}/reverse", timeout=30, params={
        "lat": lat, "lon": lon, "format": "jsonv2", "zoom": zoom})
    r.raise_for_status()
    data = r.json()
    return None if "error" in data else data
Decision diagram routing three reverse-geocoding questions to a join, a nearest join or an API.
If the answer you need is an area, join. The API is for street addresses only.

Step-by-step solution

1. Decide which question you are actually asking

Three different needs get called "reverse geocoding":

Need Right tool Cost
Which admin area contains this point? spatial join to boundaries one pass over the file
Which of my own features is nearest? nearest-neighbour join one pass, with an index
What is the postal address here? reverse geocoding API one request per point

Only the third genuinely needs an API, and it is the least common of the three in analytical work.

2. For areas, get the boundary layer once

Administrative boundaries, postcode areas and statistical geographies are published as files. Downloading one and joining against it is faster, reproducible, offline, and free of usage policy.

It also gives you the code as well as the name — E09000007, FR-75, US-CA — which is what joins to every other statistical dataset. An API's address string gives you a name that has to be matched back to a code, and that matching is the address-matching problem all over again.

3. Use within for containment, and check the misses

A point-in-polygon join is a within (or intersects) predicate. Points that match nothing are the interesting output:

  • genuinely offshore or outside the coverage of the layer
  • in a gap between polygons, because the boundary layer has slivers
  • exactly on a boundary, where floating-point comparison decides

In the measured global join, 13,464,017 points produced 12,942,217 matches — 521,800 points, 3.9%, fell outside every province polygon. Most were marine features and small islands the generalised boundary layer does not cover. That number is a property of the layer, not a bug, and it is worth knowing before somebody asks why the totals do not add up.

4. For nearest-feature questions, use a nearest join and keep the distance

sjoin_nearest answers "which is closest" — but always keep the distance, because a nearest match with no distance limit always succeeds, however far away it is:

joined = gpd.sjoin_nearest(points.to_crs(27700), depots.to_crs(27700),
                           distance_col="dist_m", max_distance=5000)

Work in a projected CRS so the distance is in metres, and set max_distance so that "nearest" cannot silently mean "on the other side of the country".

5. When you do call an API, control the zoom and cache the result

Nominatim's /reverse endpoint takes a zoom parameter that decides the level of detail returned — roughly 18 for a building, 16 for a street, 10 for a city, 3 for a country. Asking for building detail when you want the city wastes precision and returns a longer, more volatile string.

The same rate limit applies as for forward geocoding: about one request per second on the public instance. Reverse geocoding a track of GPS points is bulk geocoding, and it is exactly the case the policy exists to prevent.

6. Cache on a rounded coordinate

Reverse geocoding is a natural fit for coordinate rounding: points 10 m apart return the same answer. Round the key to four decimal places (about 11 m) — or coarser for city-level questions — and a track with a thousand fixes becomes a handful of lookups.

Bar chart comparing a DuckDB join, a GeoPandas join and an API loop for reverse geocoding.
The API bar is clipped: 13.5 million requests at one per second is 156 days.

Code examples

Example 1 — labelling points with every level at once 

import geopandas as gpd


def label_with_areas(points: gpd.GeoDataFrame, layers: dict) -> gpd.GeoDataFrame:
    """layers: {'country': gdf, 'region': gdf, 'postcode': gdf}, each with name/code."""
    out = points.copy()
    for level, areas in layers.items():
        areas = areas.to_crs(points.crs)
        joined = gpd.sjoin(out[["geometry"]], areas[["name", "code", "geometry"]],
                           how="left", predicate="within")
        # a point on a boundary can match two polygons; keep the first
        joined = joined[~joined.index.duplicated(keep="first")]
        out[f"{level}_name"] = joined["name"]
        out[f"{level}_code"] = joined["code"]
        missing = out[f"{level}_name"].isna().sum()
        print(f"{level:9s}: {missing:,} of {len(out):,} unmatched "
              f"({100 * missing / len(out):.1f}%)")
    return out

The duplicate-index guard matters more than it looks. A point exactly on a shared boundary matches both polygons, and without the guard the output silently grows longer than the input — the most common way a reverse-geocoding join corrupts a row count.

Example 2 — a cached, throttled reverse geocoder 

import json, sqlite3, time


class ReverseGeocoder:
    def __init__(self, client, cache_path="reverse.sqlite", precision=4):
        self.client, self.precision = client, precision
        self.db = sqlite3.connect(cache_path)
        self.db.execute("create table if not exists rev "
                        "(key text primary key, response text, fetched_at text)")
        self.db.commit()

    def _key(self, lat, lon, zoom):
        return f"{round(lat, self.precision)},{round(lon, self.precision)},{zoom}"

    def lookup(self, lat, lon, zoom=18):
        key = self._key(lat, lon, zoom)
        row = self.db.execute("select response from rev where key = ?", (key,)).fetchone()
        if row:
            return json.loads(row[0]) if row[0] else None

        self.client._wait()
        r = self.client.session.get(
            f"{self.client.base}/reverse", timeout=30,
            params={"lat": lat, "lon": lon, "format": "jsonv2", "zoom": zoom})
        r.raise_for_status()
        data = r.json()
        data = None if "error" in data else data
        self.db.execute("insert or replace into rev values (?, ?, datetime('now'))",
                        (key, json.dumps(data) if data else None))
        self.db.commit()
        return data

At four decimal places the cache key covers about 11 m. For a GPS track sampled every second in a city, that collapses most of a journey's stationary periods into single lookups.

Example 3 — choosing between the join and the API, in code 

def reverse_geocode(points, need, areas=None, client=None):
    """One entry point that picks the cheap route whenever it exists."""
    if need in ("country", "region", "county", "postcode_area", "neighbourhood"):
        if areas is None:
            raise ValueError(f"{need} is an area question — pass a boundary layer")
        return label_with_areas(points, {need: areas})

    if need == "nearest_feature":
        raise ValueError("use sjoin_nearest with a max_distance and keep the distance")

    if need == "street_address":
        if client is None:
            raise ValueError("a street address needs a geocoding service")
        if len(points) > 1000:
            print(f"WARNING: {len(points):,} API calls at ~1/s is "
                  f"{len(points) / 3600:.1f} hours — deduplicate on rounded "
                  f"coordinates first, or use your own instance")
        return [client.lookup(p.y, p.x) for p in points.geometry]

    raise ValueError(f"unknown need {need!r}")

Making the expensive path announce its cost is a cheap guard against somebody discovering it at 3 am.

Explanation

Why the spatial join is not just faster but better

Speed is the obvious argument: 13.5 million points joined to 4,596 polygons in 61.8 seconds, versus 13.5 million HTTP requests. But the join is also more correct for area questions.

An API returns the address components it believes apply to the point, derived from its own reference data, which may use different boundaries and different vintages from the statistical geography your analysis is built on. A join to the exact boundary file your report uses is consistent by construction. When the totals have to reconcile with a published table, that consistency is the whole point.

Why unmatched points are informative rather than annoying

The measured 3.9% of GeoNames points falling outside every province polygon is a real property of the pair (points, boundaries). Some are marine; some are on small islands a generalised layer omits; some sit in slivers between polygons.

Knowing the number lets you decide. Under a percent, ignore it. Approaching four, use a coarser layer for the residue or a nearest-boundary fallback with a distance cap — and report how many rows took the fallback.

Why round the cache key but not the data

Rounding coordinates for the cache key is a deliberate statement that the answer is constant over that distance, which for a street address at 11 m is nearly always true.

Rounding the stored coordinates is a different act with a different consequence: it quantises the data permanently. Two decimal places is a 1.1 km grid; four is about 11 m. Keep the full precision in the data and round only the key.

Why zoom matters on reverse geocoding

Asking for the highest detail and then keeping only the city name looks harmless, but it makes the answer less stable: building-level results change whenever the reference data gains a shop or splits a unit, so the same coordinate can return a different string next month.

Asking at the level you actually need makes the result stable and the string shorter. It also lightens the load on a shared service, which matters when the service is free.

Checklist of five reasons a point matches no polygon in a reverse geocoding join.
Four of the five are properties of the boundary layer, not errors in the points.

Edge cases or notes

  • A point exactly on a boundary matches two polygons. Deduplicate the join result or the row count grows.
  • Use a projected CRS for anything with a distance, including max_distance in a nearest join.
  • Antimeridian and polar points break naive bounding-box filters; test with a point at longitude 179.9.
  • Marine coordinates have no address. Expect error responses and handle them as a normal outcome.
  • Different providers use different administrative vintages. If your report cites official boundaries, join to those.
  • sjoin_nearest with no max_distance always matches something, however far away.
  • Reverse geocoding a whole GPS track is bulk geocoding and will get you blocked; deduplicate on rounded coordinates first.
  • Keep both the name and the code. Names are for humans; codes join to data.

FAQ

What is reverse geocoding?

Turning a coordinate into a description of the place: an address, or the name of the area that contains it. The two need completely different tools.

Should I use an API or a spatial join?

Join, unless you need a street address. For area questions a join against the boundary file your report already uses is faster, offline, and consistent with your other tables.

How fast is a spatial join for reverse geocoding?

Measured: 13,464,017 points against 4,596 province polygons completed in 61.8 seconds in DuckDB and 117.5 seconds in GeoPandas, producing identical results.

Why do some points match no area at all?

Because they are offshore, on an island the generalised layer omits, or in a sliver between polygons. In the measured case 3.9% of points fell outside every province.

Can I reverse geocode a whole GPS track with an API?

Not against a public service — that is bulk geocoding. Deduplicate on rounded coordinates first, or join to boundaries instead.

What does the zoom parameter do?

It sets the level of detail returned: roughly 18 for a building, 10 for a city, 3 for a country. Ask at the level you need; higher detail makes the answer less stable over time.