Nearest-Neighbour Matching with GeoPandas sjoin_nearest

Turn two datasets into a bounded set of plausible pairings, with distances measured in metres, a cap that survives dense city centres, and an output that is identical on every run.

Prerequisites

Conceptual minimum

sjoin_nearest finds, for each row of the left frame, the nearest rows of the right frame, optionally within a maximum distance. Three properties decide whether it does what you want.

Distance units follow the CRS. In a geographic CRS the “distance” is in degrees, which is not a distance: a degree of longitude is 111 kilometres at the equator and 55 at sixty degrees north. A max_distance of 0.001 therefore means different things in different parts of one dataset. Reprojecting both frames to a metric CRS first is not an optimisation, it is a correctness requirement.

It returns the nearest, not all within the radius. By default it returns the single nearest match per left row. Candidate generation wants several, which means either raising the number returned or — more controllably — doing a radius join and ranking afterwards.

Ties are resolved arbitrarily. Two features at identical distance produce an order that depends on internal index layout, so two runs over the same data can differ. Sorting explicitly by distance and then by a stable identifier makes the output reproducible.

Choosing between a nearest join and a radius join for candidate generation A decision node about how many candidates each record needs, with three outcomes. When exactly one match per record is wanted and the data is sparse, a nearest join with a maximum distance is simplest. When several candidates per record are needed for scoring, a buffer-and-join approach returns everything within the radius and lets you rank and cap afterwards. When the external side has polygons rather than points, an intersection join is more meaningful than any distance-based approach. How many candidates, and what geometry? One match, or several? Scoring needs several Cap them afterwards Nearest join One candidate per record, sparse data, simplest code Buffer and join Several candidates per record, ranked and capped after Intersection join Polygon external data; containment beats distance Scoring on a single candidate throws away the runner-up gap, which is one of the most informative signals available.
The middle branch is the default for conflation precisely because the second-best candidate matters.

Runnable solution

python 
from __future__ import annotations

import logging

import geopandas as gpd
import pandas as pd

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
logger = logging.getLogger("osm.conflate.candidates")


def to_metric(frame: gpd.GeoDataFrame, epsg: int) -> gpd.GeoDataFrame:
    """Reproject to a metric CRS. Distances in degrees are not distances."""
    if frame.crs is None:
        raise ValueError("frame has no CRS; set it before reprojecting")
    return frame.to_crs(epsg=epsg)


def candidates(external: gpd.GeoDataFrame, osm: gpd.GeoDataFrame,
               radius_m: float, epsg: int, max_per_record: int = 10,
               left_id: str = "record_id",
               right_id: str = "osm_key") -> pd.DataFrame:
    """Every OSM feature within `radius_m` of each external record, ranked.

    Returns one row per candidate PAIR, with the distance and the rank, so the
    scoring stage can see the runner-up as well as the best candidate.
    """
    left = to_metric(external, epsg)
    right = to_metric(osm, epsg)
    for frame, name, key in ((left, "external", left_id), (right, "osm", right_id)):
        if key not in frame.columns:
            raise ValueError(f"{name} frame is missing the identifier {key!r}")

    # Buffer the left side and join: this returns EVERYTHING within the radius,
    # unlike a nearest join which returns only the closest.
    probe = left.copy()
    probe["geometry"] = probe.geometry.buffer(radius_m)
    pairs = gpd.sjoin(probe[[left_id, "geometry"]], right[[right_id, "geometry"]],
                      how="inner", predicate="intersects")
    logger.info("%d record(s) produced %d raw candidate pair(s)",
                len(left), len(pairs))

    # Recover the true point-to-feature distance; the buffer was only a filter.
    left_geom = left.set_index(left_id).geometry
    right_geom = right.set_index(right_id).geometry
    pairs = pairs.reset_index(drop=True)
    pairs["distance_m"] = [
        left_geom.loc[a].distance(right_geom.loc[b])
        for a, b in zip(pairs[left_id], pairs[right_id])
    ]

    # Deterministic order: distance, then the identifier, so ties never wobble.
    pairs = pairs.sort_values([left_id, "distance_m", right_id],
                              kind="mergesort").reset_index(drop=True)
    pairs["rank"] = pairs.groupby(left_id).cumcount()

    capped = pairs[pairs["rank"] < max_per_record].copy()
    dropped = len(pairs) - len(capped)
    if dropped:
        logger.warning("capped %d pair(s) beyond rank %d — those records are "
                       "dense enough to need review regardless of score",
                       dropped, max_per_record)

    per_record = capped.groupby(left_id).size()
    logger.info("candidates per record: median %.0f, max %d, %d record(s) with none",
                per_record.median() if len(per_record) else 0,
                per_record.max() if len(per_record) else 0,
                len(left) - per_record.shape[0])
    return capped[[left_id, right_id, "distance_m", "rank"]]


if __name__ == "__main__":
    logger.info("call candidates(external, osm, radius_m=120, epsg=32633)")

Step-by-step walkthrough

  1. Refuse a missing CRS. A frame without a declared CRS cannot be reprojected correctly, and silently treating its coordinates as degrees or metres is the single most damaging assumption available here.
  2. Reproject both sides to the same metric CRS. Distances then mean metres everywhere in the dataset rather than varying with latitude.
  3. Require identifiers up front. Checking for them before the expensive join turns a confusing merge error into a clear message.
  4. Buffer and intersect rather than joining nearest. A nearest join returns one match; buffering the left side and intersecting returns everything within the radius, which is what a scoring stage needs.
  5. Recompute the true distance. The buffer was only a spatial filter; the distance that matters is from the original point or geometry, not from the buffer’s edge.
  6. Sort deterministically before ranking. Distance first, then the identifier as a tie-break, with a stable sort. Without this, two runs over identical input can rank tied candidates differently.
  7. Cap per record and report it. A record with more candidates than the cap is, by definition, in a dense area where the matcher is least reliable — the warning marks it for review rather than silently truncating.
  8. Report the distribution. Median and maximum candidates per record, plus the count of records with none, are the three numbers that tell you whether the radius is sensible before any scoring runs.
How candidate count per record varies with the radius on one urban dataset Five radii with the median number of candidates each produces per external record in an urban area. At twenty five metres the median is about one candidate and many records have none. At fifty metres the median is about two. At one hundred metres the median is about five. At two hundred metres the median is about fourteen. At four hundred metres the median is about forty, at which point the scoring stage is doing most of the work the spatial filter should have done. Candidates per record grow roughly with the square of the radius radius 25 m median 1 radius 50 m median 2 radius 100 m median 5 radius 200 m median 14 radius 400 m median 40 Area grows with the square of the radius, so doubling it roughly quadruples the pairs the scoring stage has to evaluate.
Calibrate the radius from the source's positional error, then let the cap handle the dense outliers.
The three numbers that tell you whether a candidate radius is right Three diagnostic measurements taken before any scoring runs. The median candidates per record should sit in the low single digits: one means the radius is too tight to see rivals, thirty means the scoring stage is being asked to do the spatial filter's job. The share of records with no candidate at all should be small and explainable by genuine absence from the map rather than by a radius problem. The share of records hitting the per-record cap should be small and concentrated in dense urban areas. Three numbers, read before any scoring Median candidates Two to six is healthy one is too tight Records with none Small and explainable most means a CRS bug Records at the cap Small and urban scattered means trouble All three come from the candidate stage alone, which means a mis-set radius is caught before any expensive scoring has run.
Reporting these three every run makes a radius regression obvious the first time somebody changes the projection.

Verification

  • Distances are plausible metres. Spot-check a pair whose separation you can estimate; a value in the thousands where you expected tens means the CRS is wrong.
  • Two runs produce identical output. Run twice and compare; any difference means the sort is not fully deterministic.
  • The median candidate count is small. Two to six is a healthy range; a median of thirty means the radius is far too large.
  • Records with no candidates are plausible. Some is normal; most means the radius is too small or the projections disagree.
  • The capped records are genuinely dense. Sample a few and confirm they are in city centres rather than scattered arbitrarily.

Common errors and fixes

Symptom Root cause One-line fix
Distances in the thousandths Join run in a geographic CRS Reproject both frames to a metric CRS first
Radius behaves differently by latitude Degrees used as a distance Same fix: project before measuring
Only one candidate per record Nearest join used instead of a radius join Buffer the left side and intersect
Output differs between runs Ties broken by index order Sort by distance then identifier with a stable sort
Memory exhausted on a large join Radius far too large for the density Calibrate the radius; cap candidates per record
Merge fails after the join Identifier columns missing or renamed Assert both identifiers exist before joining
Every record has zero candidates The two frames are in different places Check both CRS declarations and the bounding boxes

Specification reference

GeoDataFrame.sjoin_nearest joins each geometry in the left frame to the nearest geometries in the right frame, optionally limited by max_distance, and reports the separation in a distance column. Distances are computed in the units of the frames’ coordinate reference system, so a geographic CRS yields degrees rather than a metric distance. See the GeoPandas spatial joins documentation for the join predicates, the distance column and the CRS requirements.

Frequently Asked Questions

Why reproject before a spatial join?

Because distance in a geographic CRS is measured in degrees, and a degree is not a fixed distance. A degree of longitude spans about 111 kilometres at the equator and about 55 at sixty degrees north, so a single maximum-distance value silently means two different radii in two parts of the same dataset. Reprojecting both frames to a metric CRS makes the radius mean the same thing everywhere.

Should I use a nearest join or a radius join?

A radius join, for conflation. A nearest join returns the single closest feature, which discards the runner-up — and the gap between the best and second-best candidate is one of the most informative signals available for deciding whether a match is confident or ambiguous. Buffer the external side, intersect, then rank and cap, which gives the scoring stage everything it needs.

How many candidates per record should I keep?

Enough to include the true match and its nearest rivals, which in practice is a handful. Ten is a generous cap that costs little. The more useful observation is that a record legitimately exceeding the cap is in a dense area where distance-based matching is least reliable, so those records should be flagged for review regardless of what the scoring eventually says about them.

Why do two runs produce different candidate sets?

Because ties are being broken by whatever order the spatial index happened to return. Two features at exactly the same distance have no inherent ordering, so the result depends on internal layout that can change between runs or library versions. Sorting explicitly by distance and then by a stable identifier, with a stable sort algorithm, removes the nondeterminism entirely.

Up one level: Matching OSM Features to External Datasets.