Speeding Up Nearest Neighbour Joins with a Spatial Index

Replace O(n·m) distance loops with sindex and sjoin_nearest in GeoPandas: R-tree queries, max_distance pruning, projected CRS requirements and benchmarks that show the difference.

Use gpd.sjoin_nearest(left, right, max_distance=..., distance_col="dist_m") on layers that are both in a projected CRS. Behind it sits an R-tree that prunes candidates by bounding box before any exact distance is computed, which turns a quadratic scan into something close to n log m. The .apply(lambda p: others.distance(p).min()) pattern that appears in most tutorials is the same computation without the pruning, and it is the reason feature builds take hours.

This page is about making proximity features fast. For what those features mean and how to design them, see Vector Proximity and Buffer Generation.

The index prunes candidates before any exact distance is computed On the left a single source point is connected to every one of twelve target geometries, representing a brute force scan. On the right the same point is tested against three bounding boxes, only one of which is expanded, so exact distances are computed for just two candidate geometries. brute force: n × m R-tree: prune, then measure 12 exact distance computations, every query 4 box tests + 2 exact distances only this box is expanded

Why This Fails in Geospatial ML Pipelines

The quadratic loop is usually written once, on a small development sample where it takes four seconds, and only becomes a problem at full scale — by which time it is buried in a feature module and its cost is attributed to “the data being big”. Scaling is unforgiving: ten times more sources and ten times more targets is a hundred times the work. A build that took a minute on a county takes nearly three hours on a country.

The second failure is degrees. sjoin_nearest in EPSG:4326 returns a number, and the number is in degrees. Since a degree of longitude shrinks with latitude, max_distance=0.01 is roughly 1.1 km in Kenya and 640 m in Denmark, so a “distance to nearest road” feature is silently rescaled across the study area. Models pick this up as a latitude effect, which is exactly the kind of artefact that survives a random validation split and fails a spatial one — see spatial cross-validation strategies.

Third, unbounded nearest queries on sparse targets. If a source point sits 300 km from the nearest hospital, the tree has to expand until it finds one, scanning a large fraction of the index. Worse, the resulting feature — “nearest hospital: 312 km” — is not informative; it is a proxy for remoteness that a truncated feature plus an explicit “none within 50 km” flag encodes far better.

Core Principles

  • Project first. Nearest-neighbour work belongs in a metric CRS, always.
  • Use sjoin_nearest, not .apply. It builds the index for you and runs in C.
  • Bound the search with max_distance. Faster, and the resulting nulls are honest.
  • Build the index once and reuse it. gdf.sindex is cached on the frame; do not rebuild per batch.
  • Simplify target geometries when exactness is not needed. Fewer vertices means cheaper exact tests.
  • Encode “no neighbour within range” as its own feature, rather than as a very large distance.

Production-Ready Code

from __future__ import annotations

import logging
import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import box

logger = logging.getLogger(__name__)


def require_projected(*frames: gpd.GeoDataFrame) -> None:
    """Nearest-neighbour distances are only meaningful in a metric CRS."""
    for i, gdf in enumerate(frames):
        if gdf.crs is None:
            raise ValueError(f"frame {i} has no CRS")
        if not gdf.crs.is_projected:
            raise ValueError(
                f"frame {i} is in {gdf.crs.to_string()} (geographic). Distances would be "
                "in degrees, which vary with latitude. Reproject to a metric CRS first.")


def nearest_distance(sources: gpd.GeoDataFrame, targets: gpd.GeoDataFrame,
                     max_distance: float, target_cols: list[str] | None = None,
                     prefix: str = "nearest") -> gpd.GeoDataFrame:
    """Distance from each source to its nearest target, bounded by max_distance.

    Returns `sources` with `<prefix>_dist_m`, `<prefix>_found` and any requested
    target attributes joined on. Sources with no target inside max_distance get
    NaN distance and found=False, which is more useful than a huge number.
    """
    require_projected(sources, targets)
    if sources.crs != targets.crs:
        logger.info("reprojecting targets %s -> %s",
                    targets.crs.to_string(), sources.crs.to_string())
        targets = targets.to_crs(sources.crs)

    keep = ["geometry"] + (target_cols or [])
    joined = gpd.sjoin_nearest(
        sources, targets[keep],
        how="left", max_distance=max_distance,
        distance_col=f"{prefix}_dist_m", exclusive=False,
    )
    # sjoin_nearest can emit ties; keep the first match per source index.
    joined = joined[~joined.index.duplicated(keep="first")]
    joined[f"{prefix}_found"] = joined[f"{prefix}_dist_m"].notna()

    found = float(joined[f"{prefix}_found"].mean())
    logger.info("%s: %.1f%% of sources found a target within %.0f m",
                prefix, 100 * found, max_distance)
    if found < 0.5:
        logger.warning("more than half of sources have no target within %.0f m — "
                       "is max_distance too tight, or the target layer incomplete?",
                       max_distance)
    return joined.drop(columns=[c for c in joined.columns if c.startswith("index_right")])


def candidates_in_window(gdf: gpd.GeoDataFrame, bounds: tuple) -> gpd.GeoDataFrame:
    """Subset a layer to a bounding window using its cached spatial index."""
    idx = list(gdf.sindex.query(box(*bounds), predicate="intersects"))
    return gdf.iloc[idx]

Step-by-Step Walkthrough

Step 1 — Project both layers once, at load time

import geopandas as gpd

TARGET_CRS = 32633                     # UTM 33N — metres

parcels = gpd.read_file("parcels.gpkg").to_crs(TARGET_CRS)
roads = gpd.read_file("roads.gpkg").to_crs(TARGET_CRS)
hospitals = gpd.read_file("hospitals.gpkg").to_crs(TARGET_CRS)
require_projected(parcels, roads, hospitals)

Reprojecting once at the boundary, rather than inside a feature function, means the index is built in the right space and never rebuilt.

Step 2 — Replace the loop

# Before: quadratic, no index
# parcels["road_dist"] = parcels.geometry.apply(lambda g: roads.distance(g).min())

# After: indexed, bounded
parcels = nearest_distance(parcels, roads, max_distance=5_000,
                           target_cols=["road_class"], prefix="road")
parcels = nearest_distance(parcels, hospitals, max_distance=50_000,
                           target_cols=["beds"], prefix="hospital")

Step 3 — Turn “not found” into a usable feature

import numpy as np

parcels["road_dist_m"] = parcels["road_dist_m"].fillna(5_000)     # censored at the bound
parcels["road_censored"] = ~parcels["road_found"]
parcels["hospital_dist_km"] = parcels["hospital_dist_m"] / 1000.0
parcels["hospital_none_50km"] = ~parcels["hospital_found"]

Censoring at the search bound plus an explicit flag gives the model two clean signals instead of one contaminated one. It also keeps the feature’s distribution bounded, which matters for the scaling choices in feature scaling for geospatial inputs.

Step 4 — Benchmark on your own data

import time

sample = parcels.sample(2_000, random_state=0)

t0 = time.perf_counter()
_ = sample.geometry.apply(lambda g: roads.distance(g).min())
brute = time.perf_counter() - t0

t0 = time.perf_counter()
_ = gpd.sjoin_nearest(sample, roads[["geometry"]], max_distance=5_000,
                      distance_col="d")
indexed = time.perf_counter() - t0

print(f"brute {brute:.2f}s  indexed {indexed:.2f}s  speedup {brute / indexed:.0f}x")

Run this on a 2,000-row sample rather than the full frame — the brute-force branch on a full national layer is exactly the thing you are trying to avoid.

The gap widens with every row you add Runtime is plotted against the number of source geometries for two methods. The brute force loop rises steeply and quadratically. The indexed nearest join rises almost linearly and stays near the bottom of the chart across the whole range. Runtime against source count, 100k target geometries 10k50k100k200k time .apply loop sjoin_nearest ~2.8 h ~9 s Measured on your own geometry complexity — vertex count matters as much as row count.

Verification

import geopandas as gpd
import numpy as np
import pytest
from shapely.geometry import Point


def _pts(coords, crs=32633):
    return gpd.GeoDataFrame(geometry=[Point(*c) for c in coords], crs=crs)


def test_indexed_matches_brute_force():
    """The fast path must return the same distances as the slow one."""
    rng = np.random.default_rng(0)
    src = _pts(rng.uniform(0, 10_000, (300, 2)))
    tgt = _pts(rng.uniform(0, 10_000, (500, 2)))

    brute = src.geometry.apply(lambda g: tgt.distance(g).min()).to_numpy()
    fast = nearest_distance(src, tgt, max_distance=1e9)["nearest_dist_m"].to_numpy()
    assert np.allclose(np.sort(brute), np.sort(fast), atol=1e-6)


def test_geographic_crs_is_rejected():
    src = _pts([(13.4, 52.5)], crs=4326)
    with pytest.raises(ValueError, match="geographic"):
        require_projected(src)


def test_max_distance_yields_nan_not_a_huge_number():
    src = _pts([(0, 0)])
    tgt = _pts([(100_000, 0)])
    out = nearest_distance(src, tgt, max_distance=1_000)
    assert np.isnan(out["nearest_dist_m"].iloc[0])
    assert out["nearest_found"].iloc[0] is np.False_ or not out["nearest_found"].iloc[0]


def test_one_row_out_per_row_in():
    """sjoin_nearest can duplicate on ties; the wrapper must not change row count."""
    src = _pts([(0, 0)])
    tgt = _pts([(10, 0), (-10, 0)])          # exact tie
    assert len(nearest_distance(src, tgt, max_distance=100)) == 1

The tie test is the one people skip and then get bitten by: sjoin_nearest returns all equidistant matches, so a frame with tied geometries silently gains rows and every downstream assert len(X) == len(y) fails somewhere far from the cause.

Ties duplicate rows unless you drop them A single source point sits exactly between two target geometries. The nearest join returns two rows for that one source, so the output feature table has more rows than the input, which breaks the alignment between features and labels. One source, two equidistant targets, two output rows parcel 42 90 m 90 m indexroad_iddist_m 42R-11890 42R-20490 index 42 appears twice Deduplicate on the source index immediately after the join, and assert the row count.

FAQ

Why is my .apply distance loop so slow?

Because it is quadratic: every source is compared against every target. A hundred thousand against a hundred thousand is ten billion distance computations. An R-tree tests bounding boxes first and computes exact distances for a handful of candidates per query.

Does sjoin_nearest work in EPSG:4326?

It runs, but the numbers are degrees, and a degree is 111 km at the equator and 71 km at 50° N. A single max_distance then means different things in different parts of the same dataset. Reproject first; see fixing projection mismatches in pandas GeoDataFrames.

How much does max_distance help?

Substantially, because it bounds how far the tree must expand for isolated sources — typically another two to five times on sparse target layers. It also produces nulls that are more honest than a 400 km “nearest hospital”.

Part of: Vector Proximity and Buffer Generation Part of: Spatial Feature Engineering for Machine Learning