Nearest-Neighbour Joins Explained: Distance, Ties and Search Radius

Problem statement

"Attach the nearest road to each property" sounds like one operation. It is four decisions, and the defaults make three of them for you.

Here are 293,729 properties — one point for every registered title in Leeds — and the 211,344 OpenStreetMap road segments of West Yorkshire, both in longitude and latitude as they are usually delivered (the setup is under The data used on this page):

joined = gpd.sjoin_nearest(properties, roads)
print(f"{len(properties):,} → {len(joined):,}")

293,729 → 295,073

More rows out than in, from a join that was supposed to be one-to-one. Those are exact distance ties — a property equally close to two road segments, usually two segments that meet at the point nearest to it — and sjoin_nearest returns all of them by default.

Then the distances look wrong:

joined = gpd.sjoin_nearest(properties, roads, distance_col="dist")
print(joined["dist"].describe().round(6).to_string())

count    295073.000000
mean          0.000187
std           0.000152
min           0.000000
25%           0.000110
50%           0.000164
75%           0.000234
max           0.010107

Metres? No — the layers are in EPSG:4326, so those are degrees, and GeoPandas says so in a warning that is easy to scroll past. Worse, as step 1 shows, degrees do not just change the units: they change which road is nearest.

Nearest-neighbour joins are genuinely useful and quietly full of decisions: what counts as distance, what happens on a tie, how far is too far, and whether "nearest" is even the right question.

Quick answer

import geopandas as gpd

CRS = 27700                                     # metres, before anything else
properties = properties.to_crs(CRS)
roads = roads.to_crs(CRS)

joined = gpd.sjoin_nearest(
    properties,
    roads[["road_name", "geometry"]],
    how="left",
    max_distance=500,          # bound the search
    distance_col="dist_m",     # keep the distance
)
joined = joined[~joined.index.duplicated(keep="first")]     # resolve ties

The one that bites hardest: sjoin_nearest measures in the CRS's units, so an unprojected join returns degrees and max_distance=500 means 500 degrees.

Project, define nearest, bound it, handle ties

1. Project first, always

A degree is not a length: 111 km of latitude, but only about 66 km of longitude at Leeds' 53.8°N. So a "nearest" computed in degrees is measured with a ruler that is 1.7 times longer north–south than east–west, and the nearest feature by that metric is not always the nearest on the ground. How often, on real data:

CRS = 27700                               # British National Grid, metres
in_degrees = gpd.sjoin_nearest(properties, roads).groupby(level=0)["osm_id"].first()

properties_m, roads_m = properties.to_crs(CRS), roads.to_crs(CRS)
joined = gpd.sjoin_nearest(properties_m, roads_m, distance_col="dist_m")
in_metres = joined.groupby(level=0)["osm_id"].first()

print(f"{joined['dist_m'].median():.1f} m median, {joined['dist_m'].max():.0f} m max")
print(f"{(in_degrees != in_metres).sum():,} of {len(properties):,} properties "
      f"get a different nearest road when measured in degrees")

15.1 m median, 766 m max
33,205 of 293,729 properties get a different nearest road when measured in degrees

That is not a rounding difference. Roughly one property in nine is attached to the wrong road by the unprojected join, and nothing in its output says which. properties.estimate_utm_crs() is a reasonable choice where there is no national grid. Full treatment in how to measure distance accurately.

2. Decide what "nearest" means for your geometry

sjoin_nearest measures the minimum distance between the geometries, which is not always what people expect. Two shapes built by hand make the difference plain:

from shapely.geometry import Point, box

pt = gpd.GeoDataFrame(geometry=[Point(5, 5)], crs=27700)
poly = gpd.GeoDataFrame({"id": ["A"]}, geometry=[box(20, 0, 40, 20)], crs=27700)

print(f"to the polygon edge:     {pt.geometry.iloc[0].distance(poly.geometry.iloc[0]):.1f} m")
print(f"to the polygon centroid: {pt.geometry.iloc[0].distance(poly.geometry.centroid.iloc[0]):.1f} m")

to the polygon edge:     15.0 m
to the polygon centroid: 25.5 m

Which is right depends on the question. "How far to the nearest park?" means the edge — you can enter at the boundary. "How far to the nearest town centre?" probably means the centre, so join to a point that represents it rather than to the polygon.

A point inside a polygon has distance zero, so sjoin_nearest on a point inside two overlapping polygons returns both as ties.

And for lines, "nearest" is the nearest point along the line, which is right for a road and wrong for "the nearest road junction" — that needs the junctions as points.

3. Bound the search with max_distance

Without a bound, every feature gets a nearest neighbour however far away:

print(joined["dist_m"].describe(percentiles=[.5, .9, .99, .999]).round(1).to_string())

count    295218.0
mean         16.6
std          13.3
min           0.0
50%          15.1
90%          26.3
99%          56.5
99.9%       168.6
max         766.2

In a dense city with every footpath mapped, the tail is short — but it exists, and in a layer with gaps it is where the nonsense lives: features matched to whatever happened to be closest because the real answer was "none".

for limit in (500, 100):
    bounded = gpd.sjoin_nearest(properties_m, roads_m, how="left",
                                max_distance=limit, distance_col="dist_m")
    none = bounded["dist_m"].isna().sum()
    print(f"{none:,} of {len(properties):,} have no road within {limit} m")

6 of 293,729 have no road within 500 m
1,055 of 293,729 have no road within 100 m

max_distance also makes the query faster, because the index search is bounded rather than walking the tree until it finds something.

Choose the threshold from the question, not from the data. "Within walking distance" is 400–800 m; "served by this facility" might be 5 km; "in the same postcode" is not a distance at all.

4. Handle ties explicitly

joined = gpd.sjoin_nearest(properties_m, roads_m, how="left",
                           max_distance=500, distance_col="dist_m")
dupes = joined.index.duplicated().sum()
print(f"{len(properties_m):,} in → {len(joined):,} out, {dupes:,} extra rows from ties")

293,729 in → 295,218 out, 1,489 extra rows from ties

A tie is an exact distance equality. With real coordinates that sounds rare, but on a road network it is routine: road segments share their end vertices, and a property whose nearest point on the network is a junction is exactly as close to every segment that meets there.

Three ways to resolve, and they are not equivalent:

# (a) arbitrary — fast, and the choice depends on row order
first = joined[~joined.index.duplicated(keep="first")]

# (b) by a rule — deterministic and defensible
RANK = {"motorway": 0, "trunk": 1, "primary": 2, "secondary": 3, "tertiary": 4, "residential": 5}
ruled = (joined.assign(rank=joined["road_class"].map(RANK).fillna(9))
               .sort_values(["dist_m", "rank"]).groupby(level=0).head(1))

# (c) keep them, recorded
ties = joined.groupby(level=0)["road_name"].agg(list)
recorded = properties_m.assign(nearest_roads=ties, n_tied=ties.str.len())
print(f"(a) {len(first):,} rows  (b) {len(ruled):,} rows  "
      f"(c) {(recorded['n_tied'] > 1).sum():,} properties with more than one nearest road")

(a) 293,729 rows  (b) 293,729 rows  (c) 1,460 properties with more than one nearest road

Option (a) is what most code does implicitly, and it is fine when the tied features are interchangeable. When they are not — a trunk road and a footpath meeting at the same junction — sort by something meaningful before taking the first.

5. Keep the unmatched rows

inner = gpd.sjoin_nearest(properties_m, roads_m, max_distance=100)
left = gpd.sjoin_nearest(properties_m, roads_m, how="left", max_distance=100)
print(f"inner {len(inner):,}   left {len(left):,}")

unmatched = left[left["osm_id"].isna()]
print(f"{len(unmatched):,} unmatched, extent {[round(v) for v in unmatched.total_bounds]}")

inner 294,152   left 295,207
1,055 unmatched, extent [412626, 422944, 446968, 450290]

how="inner" — the default — silently drops the properties with no road within the limit. Those are a finding: either a gap in the road data or genuinely remote properties — farmland titles, woodland, reservoirs — and both are worth knowing. If they cluster, the road layer has a hole; if they are scattered across the rural edge of the district, they are genuinely remote.

6. Know when nearest is the wrong question

Nearest gives one answer, and sometimes the question wants a different shape:

# "everything within 25 m" — not the nearest, all of them
nearby = gpd.sjoin(properties_m, roads_m[["geometry"]], predicate="dwithin", distance=25)

# "the nearest of a particular kind"
major = roads_m[roads_m["road_class"].isin(["motorway", "trunk", "primary"])]
to_major = gpd.sjoin_nearest(properties_m, major, how="left", max_distance=2_000, distance_col="dist_m")

# "the nearest excluding itself" — for a self-join
neighbour = gpd.sjoin_nearest(properties_m, properties_m[["geometry"]], exclusive=True, distance_col="dist_m")

print(f"{nearby.index.value_counts().median():.0f} road segments within 25 m of the median property")
print(f"median distance to a motorway, trunk or primary road: {to_major['dist_m'].median():.0f} m")
print(f"median distance to the nearest other title: {neighbour['dist_m'].median():.1f} m")

2 road segments within 25 m of the median property
median distance to a motorway, trunk or primary road: 381 m
median distance to the nearest other title: 9.2 m

exclusive=True is essential for a self-join: without it every feature's nearest neighbour is itself, at distance zero. And "the three nearest" is not something sjoin_nearest or sindex.nearest can answer — both return the single nearest (plus ties). Example 2 builds it.

Leeds properties to roads with every decision explicit

The data used on this page

Leeds' properties are HM Land Registry's INSPIRE Index Polygons — the extent of every registered freehold title — reduced to one point inside each. The roads, and the points of interest in Example 2, are OpenStreetMap's, from Geofabrik's West Yorkshire extract of 1 January 2026. Both are converted to longitude and latitude, as a web service would deliver them, so that step 1 has something to fix. Land Registry data © Crown copyright and database right 2026, Open Government Licence v3.0; OpenStreetMap data © OpenStreetMap contributors, ODbL.

import geopandas as gpd
import numpy as np
import pandas as pd
import pyogrio
import requests

# HM Land Registry's download sets a cookie and redirects to itself, so use a session
reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
                               timeout=600)
reply.raise_for_status()
with open("leeds_inspire.zip", "wb") as f:
    f.write(reply.content)
pyogrio.set_gdal_config_options({"GML_DOWNLOAD_SCHEMA": "NO"})   # the file's schema URL is unreachable
titles = gpd.read_file("zip://leeds_inspire.zip!Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"])
properties = titles.set_geometry(titles.representative_point())[["INSPIREID", "geometry"]].to_crs(4326)

OSM = "https://download.geofabrik.de/europe/united-kingdom/england/west-yorkshire-260101-free.shp.zip"
with open("west-yorkshire.zip", "wb") as f:
    f.write(requests.get(OSM, timeout=600).content)
roads = (gpd.read_file("zip://west-yorkshire.zip!gis_osm_roads_free_1.shp", columns=["osm_id", "name", "fclass"])
            .rename(columns={"name": "road_name", "fclass": "road_class"}))
amenities = (gpd.read_file("zip://west-yorkshire.zip!gis_osm_pois_free_1.shp", columns=["name", "fclass"])
                .rename(columns={"fclass": "type"}))
print(f"{len(properties):,} properties, {len(roads):,} road segments, {len(amenities):,} points of interest, "
      f"all in EPSG:{roads.crs.to_epsg()}")

293,729 properties, 211,344 road segments, 24,448 points of interest, all in EPSG:4326

Example 1: a nearest join with every decision made explicit

def nearest_join(left, right, right_cols, *, max_distance=None, crs=None,
                 tie_breaker=None, measure_to="geometry", how="left",
                 distance_col="distance_m", verbose=True):
    """Nearest-neighbour join with explicit units, ties and bounds.

    measure_to: 'geometry' (closest point) or 'centroid'.
    tie_breaker: column in `right` to sort by when distances tie exactly.
    """
    target = gpd.GeoSeries(crs=crs).crs if crs else (
        left.crs if left.crs and left.crs.is_projected else left.estimate_utm_crs())
    if target.is_geographic:
        raise ValueError(f"{target.name} is geographic — distances would be in degrees")

    l = left.to_crs(target)
    r = right.to_crs(target)
    if measure_to == "centroid":
        r = r.set_geometry(r.geometry.representative_point())

    keep = [*right_cols, r.geometry.name]
    if tie_breaker and tie_breaker not in right_cols:
        keep.insert(0, tie_breaker)

    joined = gpd.sjoin_nearest(
        l, r[keep], how=how,
        max_distance=max_distance, distance_col=distance_col)

    n_raw = len(joined)
    if tie_breaker:
        joined = (joined.sort_values([distance_col, tie_breaker])
                        .groupby(level=0).head(1))
    else:
        joined = joined[~joined.index.duplicated(keep="first")]
    ties = n_raw - len(joined)

    out = left.copy()
    for col in [*right_cols, distance_col]:
        out[col] = joined[col].reindex(out.index).values

    if verbose:
        matched = out[distance_col].notna()
        print(f"  CRS            {target.name} ({target.axis_info[0].unit_name})")
        print(f"  measured to    the {measure_to}")
        print(f"  max_distance   {max_distance if max_distance else 'unbounded'}")
        print(f"  rows           {len(left):,} in → {len(out):,} out")
        print(f"  ties resolved  {ties:,}"
              + (f" by {tie_breaker}" if tie_breaker else " arbitrarily"))
        print(f"  unmatched      {(~matched).sum():,} "
              f"({100 * (~matched).mean():.2f}%)")
        if matched.any():
            d = out.loc[matched, distance_col]
            print(f"  distance       median {d.median():,.1f}  "
                  f"p95 {d.quantile(0.95):,.1f}  max {d.max():,.1f}")
    assert len(out) == len(left), "cardinality guarantee violated"
    return out


result = nearest_join(
    properties, roads, ["road_name", "road_class"],
    max_distance=500, crs=27700, tie_breaker="road_class")

  CRS            OSGB36 / British National Grid (metre)
  measured to    the geometry
  max_distance   500
  rows           293,729 in → 293,729 out
  ties resolved  1,489 by road_class
  unmatched      6 (0.00%)
  distance       median 15.1  p95 31.5  max 469.6

Every line of that report corresponds to a decision the default would have made silently. Raising on a geographic CRS is the most valuable one — a join in degrees produces plausible numbers and, as step 1 measured, the wrong neighbour for a tenth of the rows.

tie_breaker="road_class" makes the resolution deterministic: on a tie the alphabetically first class wins, which is arbitrary but stable, so re-running gives the same answer. Arbitrary-but-stable beats arbitrary-and-varying, since the latter makes a pipeline non-reproducible.

Example 2: k-nearest neighbours with the index

sjoin_nearest gives one match. For several, use the spatial index to gather candidates within a radius, then rank them exactly:

import shapely


def k_nearest(left, right, k=3, *, max_distance, right_cols=(), crs=27700):
    """The k nearest right features for each left feature within max_distance, with distances."""
    l = left.to_crs(crs)
    r = right.to_crs(crs)

    # every pair within the radius — dwithin is exact, unlike a buffer polygon
    li, ri = r.sindex.query(l.geometry.values, predicate="dwithin", distance=max_distance)
    pairs = pd.DataFrame({"left": li, "right": ri,
                          "distance_m": shapely.distance(l.geometry.values[li], r.geometry.values[ri])})
    pairs = pairs.sort_values(["left", "distance_m"])
    pairs["rank"] = pairs.groupby("left").cumcount() + 1
    pairs = pairs[pairs["rank"] <= k]

    out = pd.DataFrame({"left_index": l.index[pairs["left"]], "rank": pairs["rank"].to_numpy(),
                        "distance_m": pairs["distance_m"].round(1).to_numpy()})
    for col in right_cols:
        out[col] = r[col].to_numpy()[pairs["right"]]
    print(f"{len(left):,} left features → {len(out):,} pairs "
          f"({out.groupby('left_index').size().mean():.2f} each on average, "
          f"{len(left) - out['left_index'].nunique():,} with none)")
    return out


pairs = k_nearest(properties.head(5_000), amenities, k=3, max_distance=1_000, right_cols=("name", "type"))
print(pairs.head(6).to_string(index=False))

5,000 left features → 14,831 pairs (2.97 each on average, 5 with none)
 left_index  rank  distance_m name     type
          0     1       525.5  NaN post_box
          0     2       896.6  NaN  shelter
          0     3       976.7  NaN post_box
          1     1       297.2  NaN post_box
          1     2       685.9  NaN  shelter
          1     3       789.8  NaN post_box

The index is used only to gather candidates; the ranking uses shapely.distance on them. The index tests bounding boxes, so its ordering is not the true distance ordering — ranking on it directly would return the wrong neighbours. The mechanism is in how to use the spatial index directly, including why the candidate query must be dwithin rather than a buffer: a buffer polygon sits inside the true circle and silently drops neighbours near the limit.

An average under 3.0 means some properties had fewer than three points of interest within a kilometre, which the pair count records honestly.

Example 3: validating a nearest join

def validate_nearest(left, right, joined, *, distance_col="distance_m",
                     sample=200, seed=0, max_distance=None):
    """Check a nearest join by brute force on a random sample (both layers projected)."""
    rng = np.random.default_rng(seed)
    idx = rng.choice(len(left), min(sample, len(left)), replace=False)

    wrong = 0
    for pos in idx:
        geom = left.geometry.iloc[pos]
        label = left.index[pos]
        reported = joined.loc[label, distance_col]
        true_dist = right.geometry.distance(geom).min()         # brute force: every right feature
        if max_distance is not None and true_dist > max_distance:
            if not np.isnan(reported):
                wrong += 1
                print(f"  ✗ {label}: reported {reported:.1f} m but the true nearest "
                      f"is {true_dist:.1f} m, beyond max_distance")
            continue
        if np.isnan(reported) or abs(reported - true_dist) > 1e-6:
            wrong += 1
            print(f"  ✗ {label}: reported {reported}, brute force {true_dist:.3f} m")

    print(f"{len(idx)} sampled, {wrong} disagreements")

    d = joined[distance_col].dropna()
    print(f"distance distribution: median {d.median():,.1f}  "
          f"p95 {d.quantile(.95):,.1f}  max {d.max():,.1f}")
    if max_distance and d.max() > max_distance + 1e-6:
        print("  ⚠ a distance exceeds max_distance — check the CRS")
    if (d == 0).any():
        print(f"  · {(d == 0).sum():,} at distance zero — features that overlap "
              f"or touch, which is normal for points on a line")
    return wrong


validate_nearest(properties.to_crs(27700), roads.to_crs(27700), result, max_distance=500)

200 sampled, 0 disagreements
distance distribution: median 15.1  p95 31.5  max 469.6

Brute force on 200 samples is cheap and settles the question that matters: is the indexed result actually the nearest? A disagreement usually means a CRS mismatch between the two layers, or a max_distance interpreted in the wrong units — which is why both layers are projected before they are passed in, even though result keeps the input's own coordinates.

Every word of 'which is closest' hides a decision

A nearest-neighbour join answers "which feature is closest", and every part of that sentence hides a decision.

"Closest" depends on the metric, and the metric comes from the CRS. In a projected CRS, distance is Euclidean in metres and behaves the way intuition expects. In a geographic CRS, the coordinates are angles and Euclidean distance over them measures with a ruler whose length depends on direction — 111 km per unit north–south, about 66 km east–west in Leeds. That does not merely scale the numbers; it changes which feature is nearest, because a metre east–west counts as more degrees than a metre north–south, so the degree metric favours neighbours to the north and south. On the Leeds properties it picked a different road for 33,205 of 293,729. This is the one error in this article that changes the answer rather than the units.

"Closest" also depends on what is measured to. Shapely's distance returns the minimum distance between two geometries — the gap between their closest points. For a point and a polygon that is the distance to the boundary, and zero if the point is inside. That is usually right for accessibility questions and wrong for "how far to the town centre", which wants a centroid. The function cannot know which you mean, so it picks the geometrically natural one.

Ties are a genuine phenomenon, not a bug. When two features are at exactly equal distance, there is no principled way to prefer one, so sjoin_nearest returns both — which breaks the one-to-one cardinality people assume. Exact equality in floating-point arithmetic is rare with real coordinates and common with snapped or gridded data, which is why this surprises people intermittently rather than never. The fix is a deliberate tie-breaker, and "stable but arbitrary" is a real improvement over "arbitrary and dependent on row order", because it makes re-runs reproducible.

The unbounded search is the decision people forget. Without max_distance, "nearest" always succeeds — every feature has a closest feature somewhere, however far away. Those long-distance matches are not information; they are the algorithm doing what it was told in a region where the answer should have been "none". Bounding the search converts them into nulls, which is both more honest and faster, since the index search terminates instead of walking the tree.

And how="inner" compounds it. The default drops features with no match, so bounding the search and leaving the join type at its default silently removes exactly the features you just learned something about. how="left" keeps them as nulls, where they are visible and countable.

Finally, the question behind the question: nearest is a specific shape of answer, and often not the useful one. "The nearest road" tells you little if the second-nearest is a motorway two metres further. "Everything within 500 m", "the three nearest", "the nearest of a particular class" and "the nearest excluding itself" are all different queries, and the last one — exclusive=True — is essential for any self-join, since without it every feature's nearest neighbour is itself.

CRS units, self-joins and distance to polygon edges

• sjoin_nearest measures in the CRS's units. In EPSG:4326 that is degrees, and max_distance=500 means 500 degrees.
• Exact ties return several rows. Deduplicate on the index, ideally with a deterministic sort.
• how="inner" is the default and drops features with no match within max_distance.
• exclusive=True is mandatory for a self-join, or every feature matches itself at distance zero.
• Distance to a polygon is to its edge, and zero for a point inside it.
• max_distance speeds up the query as well as bounding the answer.
• The spatial index orders by bounding box, not true distance. Use it to narrow candidates, then rank exactly.
• distance_col is not added unless you name it — and you almost always want it.
• Web Mercator distances are inflated by 1/cos(latitude) — 69% at Leeds' 53.8°N. Never use it for distance.
• For k-nearest, sindex.nearest or a manual candidate search — sjoin_nearest gives one.

Internal links

• How to find the nearest point in GeoPandas — the practical guide
• Spatial join cardinality explained — ties as a cardinality problem
• How to use the spatial index directly in GeoPandas (sindex) — k-nearest and radius queries
• How to measure distance accurately in Python — what the units mean
• How to do a spatial join in PostGIS with SQL — the KNN operator
• Spatial predicates explained — when a predicate beats a distance
• How to perform a spatial join in Python (GeoPandas) — the predicate-based join
• Attribute join or spatial join? — when neither is needed

FAQ

Why did my nearest join return more rows than I put in?

Exact distance ties. sjoin_nearest returns every equally-near feature. Deduplicate on the index, preferably with a deterministic sort so re-runs agree.

Why are my distances tiny numbers like 0.0004?

The layers are in a geographic CRS, so distances are in degrees. Reproject both to a projected CRS before joining — otherwise max_distance is in degrees too.

Should I set max_distance?

Almost always. Without it every feature gets a nearest neighbour however far away, so a property far from any mapped road is still matched to one. It also makes the query faster.

Why do features disappear from my result?

how="inner" is the default and drops anything with no match within max_distance. Use how="left" so they survive as nulls.

Is distance measured to a polygon's edge or its centre?

To the closest point of the geometry, which for a polygon is its edge — and zero if you are inside it. Join on centroids if you want centre-to-centre.

How do I find the three nearest, not just one?

sjoin_nearest returns one. Use sindex to get a candidate set within a radius, compute exact distances, and take the top k.

How do I do a nearest join within one layer?

gpd.sjoin_nearest(gdf, gdf, exclusive=True). Without exclusive=True every feature's nearest neighbour is itself at distance zero.