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.
joined = gpd.sjoin_nearest(properties, roads)
print(len(properties), "→", len(joined))184204 → 184391
187 extra rows, from a join that was supposed to be one-to-one. Those are exact distance ties — properties equidistant from two roads — 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())mean 0.000412
max 0.008841
Metres? No — the layers are in EPSG:4326, so those are degrees. And a property 40 km from any road still gets a nearest road, because nothing bounded the search.
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
exclusive=False,
)
joined = joined[~joined.index.duplicated(keep="first")] # resolve ties
| Decision | Argument | Default |
|---|---|---|
| what units is distance in? | the CRS — project first | whatever the data is in |
| what happens on a tie? | deduplicate afterwards | all ties returned |
| how far is too far? | max_distance= |
unbounded |
| keep unmatched rows? | how="left" |
"inner" — drops them |
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.
Step-by-step solution
1. Project first, always
print(properties.crs, roads.crs) # EPSG:4326 EPSG:4326
joined = gpd.sjoin_nearest(properties, roads, distance_col="dist")
print(joined["dist"].max()) # 0.008841 — degreesA degree is not a length: 111 km of latitude, 69 km of longitude at 51.5°N. So a "nearest" computed in degrees is measured with a ruler that is 1.6 times longer north–south than east–west, and the nearest feature by that metric is sometimes not the nearest on the ground.
CRS = properties.estimate_utm_crs() # or a national grid
properties = properties.to_crs(CRS)
roads = roads.to_crs(CRS)
joined = gpd.sjoin_nearest(properties, roads, distance_col="dist_m")
print(f"{joined['dist_m'].median():.1f} m median, {joined['dist_m'].max():.0f} m max")18.4 m median, 41,882 m max
Now the numbers mean something — and the 41 km maximum is itself a finding, which step 3 addresses. 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:
from shapely.geometry import Point, box
import geopandas as gpd
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.0 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.
A point inside a polygon has distance zero, so sjoin_nearest on a point inside two overlapping polygons returns both as ties.
To measure to centroids, join on centroids:
roads_c = roads.set_geometry(roads.geometry.representative_point())
joined = gpd.sjoin_nearest(properties, roads_c, distance_col="dist_m")And note that for lines, "nearest" is the nearest point along the line, which is usually right for a road but not for "the nearest road junction".
3. Bound the search with max_distance
Without a bound, every feature gets a nearest neighbour however far away:
unbounded = gpd.sjoin_nearest(properties, roads, distance_col="dist_m")
print(unbounded["dist_m"].describe(percentiles=[.5, .9, .99, .999]))count 184,391
mean 184.2
50% 18.4
99% 1,204.8
99.9% 12,441.0
max 41,882.0
The median is 18 m and the maximum is 42 km. Those long-distance matches are not useful information — they are properties in an area with no road data, matched to whatever happened to be closest.
bounded = gpd.sjoin_nearest(properties, roads, how="left",
max_distance=500, distance_col="dist_m")
print(f"{bounded['dist_m'].isna().sum():,} of {len(properties):,} "
f"had no road within 500 m")1,204 of 184,204 had no road within 500 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, roads, how="left",
max_distance=500, distance_col="dist_m")
dupes = joined.index.duplicated().sum()
print(f"{len(properties):,} in → {len(joined):,} out, {dupes:,} tied")184,204 in → 184,391 out, 187 tied
A tie is an exact distance equality — a property on the centre line between two parallel roads, or equidistant from two ends of a junction. It is rare with real coordinates and common with data snapped to a grid.
Three ways to resolve, and they are not equivalent:
# (a) arbitrary — fast, and the choice depends on row order
joined = joined[~joined.index.duplicated(keep="first")]
# (b) by a rule — deterministic and defensible
joined = (joined.sort_values(["dist_m", "road_class"])
.groupby(level=0).head(1))
# (c) keep them, recorded
ties = joined.groupby(level=0)["road_name"].agg(list)
properties["nearest_roads"] = ties
properties["n_tied"] = ties.str.len()Option (a) is what most code does implicitly, and it is fine when the tied features are interchangeable. When they are not — a motorway and a footpath equidistant from a property — sort by something meaningful before taking the first.
5. Keep the unmatched rows
inner = gpd.sjoin_nearest(properties, roads, max_distance=500)
left = gpd.sjoin_nearest(properties, roads, how="left", max_distance=500)
print(f"inner {len(inner):,} left {len(left):,}")inner 183,187 left 184,391
how="inner" — the default — silently drops the 1,204 properties with no road within 500 m. Those are a finding: either a gap in the road data or genuinely remote properties, and both are worth knowing.
unmatched = left[left["road_name"].isna()]
print(f"{len(unmatched):,} unmatched, extent "
f"{[round(v) for v in unmatched.total_bounds]}")If they cluster, the road layer has a hole. If they are scattered at the edge of the study area, they are genuinely outside it.
6. Know when nearest is the wrong question
Nearest gives one answer, and sometimes the question wants a different shape:
# "everything within 500 m" — not the nearest, all of them
nearby = gpd.sjoin(properties, roads.assign(geometry=roads.buffer(500)),
predicate="within")
# "the three nearest" — k-nearest, via the index
tree = roads.sindex
idx = tree.nearest(properties.geometry, return_all=False, max_distance=500)
# "the nearest of a particular kind"
major = roads[roads["road_class"].isin(["motorway", "primary"])]
joined = gpd.sjoin_nearest(properties, major, how="left", max_distance=2_000)
# "the nearest excluding itself" — for a self-join
joined = gpd.sjoin_nearest(sites, sites, exclusive=True, distance_col="dist_m")exclusive=True is essential for a self-join: without it every feature's nearest neighbour is itself, at distance zero.
Code examples
Example 1: a nearest join with every decision made explicit
import numpy as np
import geopandas as gpd
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 = crs or (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)
joined = joined.reindex(l.index)
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
properties = 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 184,204 in → 184,204 out
ties resolved 187 by road_class
unmatched 1,204 (0.65%)
distance median 18.4 p95 142.8 max 499.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 sometimes the wrong neighbour.
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 directly:
import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
def k_nearest(left, right, k=3, *, max_distance=None, right_cols=(), crs=None):
"""The k nearest right features for each left feature, with distances."""
target = crs or (left.crs if left.crs and left.crs.is_projected
else left.estimate_utm_crs())
l = left.to_crs(target)
r = right.to_crs(target)
tree = r.sindex
l_geoms = l.geometry.values
r_geoms = r.geometry.values
rows = []
for pos, geom in enumerate(l_geoms):
# a generous candidate set from the index, then exact ranking
radius = max_distance or 1_000
for attempt in range(4):
search = shapely.buffer(geom, radius, quad_segs=2)
cand = tree.query(search, predicate="intersects")
if len(cand) >= k or max_distance:
break
radius *= 4
if len(cand) == 0:
continue
d = shapely.distance(geom, r_geoms[cand])
if max_distance is not None:
keep = d <= max_distance
cand, d = cand[keep], d[keep]
order = np.argsort(d)[:k]
for rank, j in enumerate(order, start=1):
row = {"left_index": l.index[pos], "rank": rank,
"distance_m": float(d[j])}
for col in right_cols:
row[col] = r.iloc[cand[j]][col]
rows.append(row)
out = pd.DataFrame(rows)
print(f"{len(left):,} left features → {len(out):,} pairs "
f"({out.groupby('left_index').size().mean():.1f} each on average)")
return out
pairs = k_nearest(properties.head(5_000), amenities, k=3,
max_distance=1_000, right_cols=("name", "type"), crs=27700)
print(pairs.head(6).to_string(index=False))5,000 left features → 14,204 pairs (2.9 each on average)
left_index rank distance_m name type
0 1 84.2 Ancoats Park park
0 2 204.8 Post Office service
0 3 412.1 Health Centre health
1 1 118.4 Ancoats Park park
Two details. The expanding search radius handles sparse areas: if the first attempt finds fewer than k candidates it quadruples the radius, up to four attempts. With max_distance set that expansion is skipped, because the bound is the answer.
And the index is used only to narrow the candidates; the exact ranking uses shapely.distance on that subset. 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.
Average 2.9 rather than 3.0 means some properties had fewer than three amenities within a kilometre, which the pair count records honestly.
Example 3: validating a nearest join
import numpy as np
import geopandas as gpd
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."""
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]
# brute force: every right feature
true_dist = right.geometry.distance(geom).min()
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(f" ⚠ 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 inside polygons")
return wrong
validate_nearest(properties, roads, properties, max_distance=500)200 sampled, 0 disagreements
distance distribution: median 18.4 p95 142.8 max 499.6
· 412 at distance zero — features that overlap or touch, which is normal
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.
The zero-distance note is worth keeping. Distance zero is not an error — a property whose polygon touches a road, or a point inside a park, is genuinely zero away. Flagging it as suspicious would produce noise; explaining it prevents someone else raising it as a bug.
Explanation
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, 69 km east–west at British latitudes. That does not merely scale the numbers; it can change which feature is nearest, because the anisotropy favours east–west neighbours. 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, even 42 km 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.
Edge cases or notes
sjoin_nearestmeasures in the CRS's units. In EPSG:4326 that is degrees, andmax_distance=500means 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 withinmax_distance.exclusive=Trueis 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_distancespeeds 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_colis not added unless you name it — and you almost always want it.- Web Mercator distances are inflated by 1/cos(latitude) — 61% at 51.5°N. Never use it here.
- For k-nearest,
sindex.nearestor a manual candidate search —sjoin_nearestgives 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 42 km from any road is 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.