sjoin_nearest vs cKDTree: Which Nearest-Neighbor Join to Use

GeoPandas' sjoin_nearest and SciPy's cKDTree both answer "which feature is closest," but they make different trade-offs in convenience, geometry support, and speed. This guide benchmarks them on the same data so you can pick deliberately. It is for anyone matching points to their nearest feature at scale. It sits under Nearest Neighbor & KD-Tree Search in Spatial Analysis & Advanced Query Techniques, and it is the head-to-head companion to that guide's broader survey of KD-tree, BallTree, and STRtree indexes.

sjoin_nearest versus cKDTree comparison A side-by-side matrix: sjoin_nearest handles any geometry, keeps attributes and distance, and is the convenient default but slower; cKDTree is 5 to 10 times faster on points only, returns array indices, and needs a metric CRS. gpd.sjoin_nearest the convenient default scipy cKDTree raw point speed Any geometry (lines, polys) Points / centroids only Keeps attributes + distance Returns array indices R-tree under the hood Balanced KD-tree k = 1 nearest by default k-nearest is trivial ~1x baseline speed ~5-10x faster, points only Rule of thumb: sjoin_nearest unless points-only speed dominates — then cKDTree
Both solve nearest-neighbor matching; the choice is convenience and geometry support versus raw point speed.

Why This Approach / What Goes Wrong

sjoin_nearest is the same spatial-join machinery in nearest-neighbor mode: it works on any geometry type, keeps both layers' attributes, and returns the distance — it is the right default. But it builds an R-tree and carries DataFrame overhead, so for millions of points it is slower than necessary. cKDTree is a tight C implementation that finds nearest neighbors among coordinate arrays in a fraction of the time, but it only understands points (you must reduce polygons/lines to representative points), and it returns positional indices you map back yourself. The common error with cKDTree is feeding it geographic degrees: KD-tree distances are Euclidean, so coordinates must sit in a projected, metric CRS or the "nearest" answer is wrong near the poles.

Decision path from data shape to nearest-neighbour tool Three questions asked in order down the left column, each with a yes branch to an answer on the right. If the targets are lines, polygons or mixed geometry, the answer is sjoin_nearest. If you need both layers' attributes and a distance column for free, the answer is again sjoin_nearest. If instead the join is millions of points, or you need more than one neighbour, the answer is cKDTree on a metric CRS. Answering no to all three leaves sjoin_nearest as the default. Pick from the data shape, not from the benchmark Are the targets non-point geometry? lines, polygons, mixed layers sjoin_nearest the R-tree indexes any geometry Do you need both sets of attributes? plus the measured distance sjoin_nearest distance_col fills itself in Millions of points, or k > 1? the time budget dominates scipy cKDTree reproject to a metric CRS first sjoin_nearest stays the default yes yes yes no no no
Speed is the last question, not the first: geometry type and the attributes you need decide the tool long before the wall-clock difference does.

Where the time actually goes

The wall-clock gap is not mostly the tree. Both libraries descend a balanced structure in logarithmic time, and on identical point data the raw search cost is within a small factor. What sjoin_nearest adds is everything around the search: it builds an STRtree over envelopes rather than raw coordinates, calls Shapely 2.0's query_nearest, and then does DataFrame work — aligning two indexes, concatenating the right frame's columns onto the left, allocating the distance column, and materialising a result frame that carries every column of both inputs. On a left frame with sixty attribute columns, that concatenation can cost more than the neighbour search it followed.

That is useful, because it tells you where the cheap wins are. Selecting only the columns you need from the right layer — lockers[["locker_id", "geometry"]] rather than lockers — routinely halves the runtime without changing the tool. The benchmark below deliberately passes the whole frame in the timed call to show the untuned number.

The behaviours that differ, not just the timings

Three semantic differences decide correctness rather than speed, and each has caught production code.

Ties multiply rows. sjoin_nearest returns all equidistant matches, so a point exactly between two lockers produces two output rows. Any left frame with grid-snapped coordinates or duplicated targets will therefore come back longer than it went in. cKDTree.query picks one — deterministically for a given tree, but the choice is an artefact of construction order, not a rule you should rely on.

Unmatched rows are represented differently. With how="left" and max_distance set, sjoin_nearest leaves NaN in the right-hand columns for points with no target in range. cKDTree with distance_upper_bound signals the same condition with an infinite distance and an index equal to the tree's length — a sentinel that raises IndexError if you feed it straight into a lookup, as Nearest Neighbor & KD-Tree Search covers in detail.

Self-matches need opting out. When both layers are the same frame — deduplicating sensors, finding each parcel's nearest neighbouring parcel — sjoin_nearest matches every feature to itself at distance zero unless you pass exclusive=True. The KD-tree equivalent is asking for k=2 and discarding the first column.

Prerequisites

conda install -c conda-forge "geopandas=0.14.*" "scipy=1.11.*" "numpy=1.26.*" "shapely=2.0.*"

Install both libraries from the same channel. A pip SciPy wheel beside a conda NumPy is the usual cause of an ABI error at from scipy.spatial import cKDTree, which reads as a missing package but is a binary-compatibility problem.

Step-by-Step Implementation

1. Load points and candidate targets in a projected CRS.

import geopandas as gpd

# 500k delivery addresses, 8k pickup lockers
addresses = gpd.read_file("addresses.gpkg").to_crs(epsg=25832)
lockers = gpd.read_file("lockers.gpkg").to_crs(epsg=25832)
assert addresses.crs.is_projected and lockers.crs.is_projected
assert addresses.crs.equals(lockers.crs), "Both layers must share one CRS"

is_projected alone is not enough. Two layers can each be projected and still be in different projections — one in ETRS89 / UTM 32N, the other in a state or national grid — in which case sjoin_nearest raises a CRS-mismatch warning and cKDTree says nothing at all, because it never sees a CRS. It simply treats two unrelated coordinate systems as one plane and returns distances in the hundreds of kilometres. Assert equality, not just projection.

2. The convenient path — sjoin_nearest.

nearest_gpd = gpd.sjoin_nearest(
    addresses, lockers[["locker_id", "geometry"]],
    how="left", distance_col="dist_m",
)

3. The fast path — cKDTree on coordinate arrays.

import numpy as np
from scipy.spatial import cKDTree

locker_xy = np.column_stack([lockers.geometry.x, lockers.geometry.y])
address_xy = np.column_stack([addresses.geometry.x, addresses.geometry.y])

tree = cKDTree(locker_xy)
dist_m, idx = tree.query(address_xy, k=1)          # nearest locker per address

addresses["locker_id"] = lockers["locker_id"].to_numpy()[idx]
addresses["dist_m"] = dist_m

4. Benchmark both on the same inputs.

import time

def timed(label, fn):
    t0 = time.perf_counter()
    fn()
    print(f"{label}: {time.perf_counter() - t0:.2f}s")

timed("sjoin_nearest", lambda: gpd.sjoin_nearest(addresses, lockers, distance_col="d"))
timed("cKDTree", lambda: cKDTree(locker_xy).query(address_xy, k=1))
# sjoin_nearest: 3.41s
# cKDTree: 0.42s

5. Trim the frames and re-time before concluding anything. Most of the reported gap is DataFrame overhead, so measure the tuned version too — otherwise you will migrate to cKDTree for a speed-up you could have had for free.

lean_lockers = lockers[["locker_id", "geometry"]]

timed("sjoin_nearest (lean)", lambda: gpd.sjoin_nearest(addresses[["geometry"]], lean_lockers, distance_col="d"))
timed("cKDTree (threaded)", lambda: cKDTree(locker_xy).query(address_xy, k=1, workers=-1))
# sjoin_nearest (lean): 1.63s
# cKDTree (threaded): 0.14s

The ratio survives, but the absolute numbers move enough to change the decision. If the join runs once in a nightly batch, 1.6 seconds is not worth a different code path; if it runs per request behind an API, 0.14 seconds is.

6. Enforce one row per input feature. Ties are the reason a 500,000-row left frame comes back with 500,006 rows. Decide the rule rather than discovering it.

# Keep the first match per address, breaking ties on the lower locker id
nearest_one = (
    nearest_gpd.sort_values(["dist_m", "locker_id"])
    .groupby(level=0, sort=False)
    .first()
)
assert len(nearest_one) == len(addresses)

7. Ask for more than one neighbour. This is where the two tools genuinely diverge: cKDTree returns a rectangular (n, k) array from a single call, which reshapes cleanly into tidy long form for ranking or fallback logic. sjoin_nearest has no k parameter at all.

import pandas as pd

dist_k, idx_k = tree.query(address_xy, k=3, workers=-1)

nearest_three = pd.DataFrame({
    "address_id": np.repeat(addresses["address_id"].to_numpy(), 3),
    "rank": np.tile([1, 2, 3], len(addresses)),
    "locker_id": lockers["locker_id"].to_numpy()[idx_k].ravel(),
    "dist_m": dist_k.ravel().round(1),
})
print(nearest_three.head(3))
#    address_id  rank  locker_id  dist_m
# 0     A-10041     1     LK-2287   214.6
# 1     A-10041     2     LK-2290   688.3
# 2     A-10041     3     LK-1904  1042.7
Nearest-neighbor join wall-clock time on 500k points A horizontal bar chart of the benchmark: sjoin_nearest takes 3.41 seconds, while cKDTree takes 0.42 seconds — roughly eight times faster on the same 500,000 point to 8,000 target join. 500k addresses → 8k lockers, metric CRS gpd.sjoin_nearest scipy cKDTree 0s 1s 2s 3s 3.41s 0.42s ≈ 8x faster — points only, indices returned
The same join, timed on identical inputs: cKDTree finishes in a fraction of the time, but only because the targets reduce cleanly to points.

Verification

Both methods must agree on the matches — if they don't, a CRS or coordinate-ordering bug is present.

import numpy as np

# Compare the two results' assigned locker ids
agree = (nearest_gpd.sort_index()["locker_id"].to_numpy() == addresses["locker_id"].to_numpy())
print(f"Agreement: {agree.mean():.4%}")            # Agreement: 100.0000%
assert agree.mean() > 0.999, "Methods disagree — check CRS and (x, y) order"

# Distances should match closely (sub-metre)
assert np.allclose(nearest_gpd.sort_index()["dist_m"], addresses["dist_m"], atol=0.01)

Agreement between two indexed methods is reassuring but not conclusive — both could share a CRS mistake. Add one brute-force check on a small random sample, where the answer is computed with no index at all:

rng = np.random.default_rng(42)
sample = rng.choice(len(address_xy), size=200, replace=False)

brute = np.linalg.norm(address_xy[sample][:, None, :] - locker_xy[None, :, :], axis=2)
brute_idx = brute.argmin(axis=1)
brute_dist = brute.min(axis=1)

assert (brute_idx == idx[sample]).all(), "KD-tree disagrees with exhaustive search"
assert np.allclose(brute_dist, dist_m[sample], atol=1e-6)
print("brute-force sample agrees:", len(sample))   # brute-force sample agrees: 200

Two sanity checks are worth adding permanently. Confirm the row count survived the join — len(nearest_gpd) == len(addresses) fails loudly when ties multiplied rows — and confirm the distance distribution is physically plausible for the layer, because a median nearest-locker distance of 0.0043 is a metric result computed in degrees rather than metres.

print(f"rows in : {len(addresses):,}   rows out: {len(nearest_gpd):,}")
print(addresses["dist_m"].describe()[["min", "50%", "max"]].round(1).to_dict())
# {'min': 3.9, '50%': 486.2, 'max': 9184.5}

Edge Cases & Debugging

Frequently Asked Questions

Which should I reach for by default? sjoin_nearest. It handles any geometry, keeps both layers' attributes, fills in the distance column, and produces a GeoDataFrame you can keep working with. Reach for cKDTree when a specific constraint pushes you there: both layers are genuinely point-like, you need more than one neighbour, or the join sits on a latency budget rather than in a batch.

Is the 8× speed-up real, or an artefact of the benchmark? Real in shape, exaggerated in magnitude by an untuned call. Passing whole frames makes sjoin_nearest copy every column of both layers into the result; trimming to the columns you need closes roughly half the gap, as step 5 shows. Benchmark your own frames — the ratio depends far more on column count and row width than on point count.

Can I keep the attributes from both layers with cKDTree? Yes, but you write the join yourself: the tree returns positional indices into the target array, so you attach columns with target["col"].to_numpy()[idx]. That is one line for one column and an increasingly error-prone one as columns multiply, which is precisely the convenience sjoin_nearest is selling.

Does either method need a projected CRS? Both do, for the same reason and with different consequences. cKDTree computes planar Euclidean distance and will happily treat degrees as a length unit. sjoin_nearest measures with GEOS, which is also planar, so its distance_col in EPSG:4326 is in degrees too — the ranking merely happens to survive over small extents, which makes the error easier to miss. Project both layers to a suitable metric CRS with Coordinate Systems with PyProj before either call.

What if the targets are polygons and I still want KD-tree speed? Use the tree as a shortlist and refine: query the k nearest representative points, then compute exact shapely.distance against those few candidate polygons and take the minimum. You keep logarithmic search with an exact answer, at the cost of one small vectorized distance pass per candidate column.

When does neither of them scale, and what replaces them? When the point set no longer fits comfortably in one process. At that point the question stops being which in-memory index to use: push the join into PostGIS, where the <-> operator answers ordered nearest-neighbour queries against a GiST index server-side, or partition the work with Dask-GeoPandas so each worker holds only nearby geometry. Both avoid materialising the whole layer in RAM, which is the actual constraint.

Related