Nearest Neighbor & KD-Tree Search in Geospatial Python

Nearest-neighbor search answers the "which feature is closest" question that sits underneath routing, facility assignment, sensor deduplication, and spatial joins — one of the core query patterns in Spatial Analysis & Advanced Query Techniques. This guide covers how a KD-tree partitions coordinate space to turn an O(n²) distance sweep into O(n log n) lookups, when to reach for scipy.spatial.cKDTree versus sklearn.neighbors.BallTree, and how to keep every result correct by projecting to a metric CRS first. It builds on the point geometries produced by Proximity & Buffer Analysis and feeds enriched matches into Geometric Intersections & Overlays; the head-to-head benchmark against GeoPandas' convenience API lives in sjoin_nearest vs cKDTree.

KD-tree spatial partitioning A plane of points recursively split by alternating vertical and horizontal lines, so a nearest-neighbour query only searches the partition containing the query point and its neighbours. Partition the plane, prune the search query split on x → then on y → search O(log n)
A KD-tree recursively halves the plane, so a nearest-neighbour query examines a small partition rather than every point.

Architecture & Data Structures

A KD-tree is a binary tree over a k-dimensional coordinate array. At each level it splits the remaining points on one axis — x, then y, then back to x — placing the median on the node and the two halves in the child subtrees. A query descends to the leaf containing the target, then backtracks upward, only crossing a splitting plane when the hypersphere of the current best distance actually reaches past it. That pruning is why the expected cost of a nearest-neighbor lookup is logarithmic in the number of points rather than linear.

Three index structures dominate practical geospatial work, and they are not interchangeable:

The minimal construction is a single call over an (n, 2) array of projected coordinates:

import numpy as np
from scipy.spatial import cKDTree

# station_xy: (n, 2) array of projected easting/northing in metres, not degrees
station_xy = np.array([[512340.0, 5403120.0], [514902.0, 5401876.0]])
station_tree = cKDTree(station_xy, leafsize=32)

The three constructor arguments that matter are rarely tuned and occasionally decisive. leafsize (default 16) is the point count below which a node stops splitting and the query falls back to a brute-force scan inside the leaf; raising it to 32–64 shallows the tree, cuts build time, and usually speeds up large batch queries because the linear scan over a small leaf vectorizes well. balanced_tree=True (the default) splits on the median, which costs a partial sort per level but guarantees O(n log n) depth; setting it to False splits on the midpoint of the bounding box instead, building noticeably faster on tens of millions of points but degrading badly when the points are clustered — which spatial data always is. compact_nodes=True shrinks each node's hyperrectangle to the data it actually holds, paying build time to make pruning tighter later.

Budget the two costs separately. Construction is O(n log n) and single-threaded; querying is O(log n) per point and parallelizable. If you build a tree over 2 million stations and query 500 points against it, the build dominates and you should reuse the tree. If you query 50 million points against it, the build is noise. The memory footprint is the coordinate array itself — n × 2 × 8 bytes — plus roughly one node per leafsize points, so a 10-million-point tree is on the order of 200 MB, cheap relative to the GeoDataFrame it came from.

One structural limit is worth knowing before you hit it: a KD-tree partitions an infinite plane, so it has no concept of wrap-around. Points either side of the antimeridian are as far apart as the coordinate system says they are, and a tree built on longitudes will happily report a station at 179.9° E as the far side of the world from one at 179.9° W. SciPy's boxsize argument enables periodic boundaries and can model that wrap in degree space, but the distances it returns are then degrees, not metres. For anything metric near the dateline, project into a CRS whose extent does not straddle it.

Choosing a spatial index: cKDTree vs BallTree vs STRtree A comparison matrix. scipy cKDTree uses axis-aligned boxes, points only, Euclidean only, best for fast planar k-nearest-neighbour. sklearn BallTree uses nested hyperspheres, points only, supports the haversine metric, best for querying directly in latitude and longitude on the sphere. shapely STRtree is a bounding-box R-tree, indexes any geometry, matches by envelope overlap, best for sjoin and mixed geometries. Pick the index that fits the query scipy cKDTree sklearn BallTree shapely STRtree Structure Geometry Custom metric (Haversine) Best use case Axis-aligned boxes Nested hyperspheres Bounding-box R-tree Points only Points only Any geometry Euclidean only Yes — haversine Envelope overlap Fast planar k-NN Query in lat/lon sjoin, mixed geoms
The three indexes are not interchangeable: cKDTree wins on planar point speed, BallTree is the only one that measures great-circle distance, and STRtree is the one that indexes non-point geometry.

Environment Configuration & Dependency Resolution

Nearest-neighbor search draws on the numerical stack (scipy, numpy, scikit-learn) and the geospatial stack (geopandas, shapely, pyproj) at once. The version floors that matter are on the geometry side: GeoPandas DataFrames 1.0 and Shapely 2.0 vectorized the .x / .y coordinate accessors and the spatial index, so extracting a coordinate array from a large layer is now a single vectorized pass rather than a Python loop. Keep GEOS, GDAL, and PROJ aligned by installing the whole stack from conda-forge — mixing a pip Shapely wheel against a conda GEOS is the usual cause of silent predicate mismatches.

# Aligned environment for nearest-neighbor pipelines
conda create -n nn-search -c conda-forge \
  "python=3.12" "geopandas>=1.0" "shapely>=2.0" "pyproj>=3.6" \
  "scipy>=1.11" "scikit-learn>=1.5" "numpy>=1.26" "joblib>=1.4"
conda activate nn-search

The workers parameter used below to parallelize queries was stabilized in SciPy 1.9, so treat scipy>=1.11 as a comfortable floor. If you plan to fall back to approximate search on very large data, add pynndescent>=0.5 or faiss-cpu in the same environment rather than a separate one, to keep the numpy ABI consistent.

Vectorized Operations & Core Workflow

The end-to-end pattern is: read a layer, project it to a metric CRS, extract a clean coordinate array, build the tree once, then run batch queries against it. Building the index once and querying many points amortizes construction cost — never rebuild the tree per query.

import geopandas as gpd
import numpy as np
from scipy.spatial import cKDTree

# 1. Load a sensor network and validate its CRS before anything else
sensors = gpd.read_file("sensor_network.geojson")
if sensors.crs is None:
    raise ValueError("Dataset has no CRS. Assign one before projecting.")

# 2. Project to a local metric CRS (UTM 33N here) — degrees are not distances
sensors_utm = sensors.to_crs(epsg=32633)

# 3. Extract a clean coordinate array; drop invalid/missing geometries first
sensors_utm = sensors_utm[sensors_utm.geometry.is_valid & sensors_utm.geometry.notna()]
sensor_xy = np.column_stack((sensors_utm.geometry.x, sensors_utm.geometry.y))
finite_mask = ~np.isnan(sensor_xy).any(axis=1)
sensor_xy = sensor_xy[finite_mask]
sensors_clean = sensors_utm.iloc[finite_mask].copy()

# 4. Build the index once (leafsize 16-32 trades traversal speed vs memory)
sensor_tree = cKDTree(sensor_xy, leafsize=32)

# 5a. k-nearest: 4 closest sensors to every sensor (col 0 is the point itself)
#     workers=-1 fans the batch across all CPU cores (SciPy >= 1.9)
knn_dist_m, knn_idx = sensor_tree.query(sensor_xy, k=4, workers=-1)

# 5b. radius: every sensor within a 500 m service radius, one array per point
within_500m = sensor_tree.query_ball_point(sensor_xy, r=500.0, workers=-1)

Two behaviours trip people up. cKDTree.query returns distances and positional indices as parallel numpy arrays — the indices are row positions in sensor_xy, not DataFrame labels, so map them back with .iloc, never .loc. And query on a point that is in the tree returns that point as its own nearest neighbor at distance zero; ask for k+1 and discard column zero when you want distinct neighbors.

Two layers, and the sentinel that means "no match"

The self-query above is the deduplication case. The far more common production case has two layers — incidents and hospitals, addresses and depots, readings and weather stations — where you build the tree over the target layer and query it with the source layer. Here the tree and the query array have different lengths, and the returned indices point into the target.

The moment a real service radius enters the picture, distance_upper_bound becomes the parameter to reach for, and it has a return convention that catches almost everyone the first time. Points with no neighbour inside the bound do not raise and do not return -1: their distance comes back as inf and their index as n, one past the last valid row of the tree data. Index that straight into a target array and NumPy raises IndexError; index it into a Python list and you get a wrap-around match to the wrong feature. Mask on the infinite distance before you use the indices.

import geopandas as gpd
import numpy as np
import pandas as pd
from scipy.spatial import cKDTree

incidents = gpd.read_file("incidents.gpkg").to_crs(epsg=25832)
hospitals = gpd.read_file("hospitals.gpkg").to_crs(epsg=25832)

hospital_xy = np.column_stack((hospitals.geometry.x, hospitals.geometry.y))
incident_xy = np.column_stack((incidents.geometry.x, incidents.geometry.y))

hospital_tree = cKDTree(hospital_xy)

# Nearest hospital within a 15 km response envelope, or nothing
dist_m, idx = hospital_tree.query(
    incident_xy, k=1, distance_upper_bound=15_000.0, workers=-1
)

matched = np.isfinite(dist_m)                 # idx == len(hospital_xy) where unmatched
print(f"unmatched incidents: {(~matched).sum():,}")   # unmatched incidents: 1,884

incidents["hospital_id"] = pd.NA
incidents.loc[matched, "hospital_id"] = hospitals["hospital_id"].to_numpy()[idx[matched]]
incidents["response_m"] = np.where(matched, dist_m, np.nan)

Radius questions have their own methods. query_ball_point answers "everything within r of each query point" and returns a ragged list of arrays. query_ball_tree answers the same question between two trees in one call, which is much faster than looping. And sparse_distance_matrix returns the pairwise distances under a cutoff as a scipy.sparse matrix — the right structure when you need the values, not just the membership, and the dense matrix would not fit.

sensor_tree = cKDTree(sensor_xy)

# All hospital-incident pairs within 2 km, as a tree-to-tree query
pairs = hospital_tree.query_ball_tree(cKDTree(incident_xy), r=2_000.0)

# Same neighbourhood, but keeping the distances, in COO sparse form
gaps = sensor_tree.sparse_distance_matrix(sensor_tree, max_distance=500.0, output_type="coo_matrix")
print(f"stored pairs: {gaps.nnz:,} of {len(sensor_xy) ** 2:,} possible")

The scaling caveat is the same for all three: output size is quadratic in local density, not in point count. Doubling the radius on a clustered urban layer can quadruple the memory the result occupies, and a radius large enough to reach most of the dataset from most points will exhaust RAM long before the query is slow enough to notice. Bound the radius to a distance that means something operationally, and check nnz before materializing anything dense.

The two parallel arrays returned by cKDTree.query with k equal to four Two aligned matrices, each three rows by four columns. The left one holds distances in metres and the right one holds positional indices, cell for cell. Column zero is shaded as the failure case: every distance there is 0.0 because the query point matched itself, and every index there is the row's own position. Columns one through three, shaded as the useful result, hold the distinct neighbours, which is why the working slice is bracket colon comma one colon. What query(sensor_xy, k=4) hands back knn_dist_m — metres knn_idx — row positions k=0 k=1 k=2 k=3 k=0 k=1 k=2 k=3 row 0 row 1 row 2 row 0 row 1 row 2 0.0 0.0 0.0 0 1 2 128.4 210.9 388.2 96.7 145.3 302.5 128.4 176.0 251.8 7 3 12 4 9 2 0 5 8 column 0 is the query point matching itself at 0.0 m — slice it off with [:, 1:] columns 1 to k-1 are the distinct neighbours, already sorted by distance Both arrays share one shape, so knn_idx[r, c] names the sensor whose gap is knn_dist_m[r, c].
The two arrays are cell-for-cell aligned, and column zero is always the self-match — the source of both the "everything is zero metres away" surprise and the off-by-one in downstream joins.

Geometry & Data Processing Details

cKDTree only understands points. Lines and polygons have to be reduced to a representative coordinate before indexing, and which coordinate you pick changes the answer. A centroid is the usual choice, but the centroid of a C-shaped or multipart polygon can fall outside the geometry entirely — use representative_point() when you need a coordinate guaranteed to lie on the feature.

import numpy as np

# Reduce parcels (polygons) to on-surface points for point-based NN search
parcels_utm = parcels.to_crs(epsg=32633)
rep_points = parcels_utm.geometry.representative_point()
parcel_xy = np.column_stack((rep_points.x, rep_points.y))

Mapping query output back to attributes is where most bugs hide. Keep the index array positional and join through .iloc:

import pandas as pd

# For each sensor, attach its single nearest *other* sensor and the gap in metres
nearest_other = knn_idx[:, 1]          # column 0 is the point itself
gap_m = knn_dist_m[:, 1]

matched = sensors_clean.reset_index(drop=True).copy()
matched["nearest_sensor_id"] = sensors_clean["sensor_id"].to_numpy()[nearest_other]
matched["gap_m"] = gap_m

The representative-point approximation, and how to undo it

Reducing a polygon to one point is not a small simplification — it changes the question from "which parcel is nearest" to "which parcel's chosen point is nearest", and those differ whenever features vary in size. A 40-hectare industrial parcel whose representative point sits 900 m away may have a boundary 60 m from the query point, while a small parcel whose point is 400 m away is 380 m away at its closest edge. Point-based search ranks them in the wrong order, and the error scales with the size disparity in the layer.

The fix does not require abandoning the KD-tree. Use it as a shortlist generator, then compute exact geometry distance on the handful of candidates it returns. Shapely 2.0's distance is vectorized, so scoring twenty candidates per query point is trivial, and the result is exact while the search stays logarithmic.

import numpy as np
import shapely

K_SHORTLIST = 20   # generous enough that the true nearest is almost always inside

# Shortlist by representative point, then rank exactly by true polygon distance
_, cand_idx = parcel_tree.query(outfall_xy, k=K_SHORTLIST, workers=-1)

parcel_geoms = parcels_utm.geometry.to_numpy()
exact_m = np.vstack([
    shapely.distance(outfall_geoms, parcel_geoms[col]) for col in cand_idx.T
]).T                                   # (n_outfalls, K_SHORTLIST) true distances

best_col = exact_m.argmin(axis=1)
nearest_parcel_pos = cand_idx[np.arange(len(cand_idx)), best_col]
nearest_m = exact_m[np.arange(len(cand_idx)), best_col]

Two checks make this safe rather than merely plausible. Confirm the shortlist is deep enough by verifying that the winning column is rarely the last one — if best_col regularly equals K_SHORTLIST - 1, the true nearest feature is falling outside the shortlist and you need a larger k. And bound the error: the maximum possible discrepancy is the largest feature's radius, so compute parcels_utm.geometry.apply(lambda g: g.hausdorff_distance(g.representative_point())).max() once on the layer and decide whether that magnitude matters for the decision the numbers feed. When it does not — sensor networks, address points, utility poles, anything effectively point-like at the analysis scale — skip the refinement entirely and keep the plain KD-tree result.

Degenerate inputs deserve one guard. Perfectly collinear or exactly-duplicated coordinates can drive cKDTree construction into pathological depth. If a dataset stacks many points on identical coordinates (co-located sensors, snapped grid nodes), either deduplicate first or add sub-millimetre jitter (np.random.uniform(-1e-6, 1e-6, size=xy.shape)) so the partitioning stays balanced.

CRS Alignment & Projection Pipeline

This is the single rule that determines whether nearest-neighbor results are correct: a KD-tree measures Euclidean distance, so its coordinates must be in a projected, metric CRS. Feed it EPSG:4326 degrees and it computes a straight-line distance in degree units — which stretches badly with latitude and is meaningless as a metre value. Pick a projection whose distances are locally accurate for your extent, not Web Mercator (EPSG:3857), whose scale error grows with latitude and corrupts any metric query. Handle the projection with Coordinate Systems with PyProj, and remember that GeoPandas' to_crs already emits x/y (easting/northing) order, so no always_xy juggling is needed at the GeoDataFrame layer.

# Assert a projected CRS before building the tree — fail loud, not silently wrong
assert sensors_utm.crs is not None and sensors_utm.crs.is_projected, (
    "KD-tree distances require a projected metric CRS (e.g. a local UTM zone)."
)

When a dataset spans several UTM zones — say a continental point layer — no single flat projection stays accurate, and forcing one skews the "nearest" answer near the zone edges. Two robust options:

import numpy as np
from sklearn.neighbors import BallTree

# Continental layer kept in geographic coords, queried on the sphere
depots_deg = depots.to_crs(epsg=4326)
lat_lon_rad = np.radians(
    np.column_stack((depots_deg.geometry.y, depots_deg.geometry.x))  # (lat, lon)
)
depot_tree = BallTree(lat_lon_rad, metric="haversine")

ang_dist, depot_idx = depot_tree.query(lat_lon_rad, k=2)
depot_gap_m = ang_dist[:, 1] * 6_371_000.0  # radians -> metres

It helps to know how much error you are trading away, because the answer is often "not enough to matter". Inside a single UTM zone the scale factor runs from 0.9996 on the central meridian to about 1.0004 at the zone edge, so a 10 km measurement is off by at most a few metres — irrelevant for ranking neighbours that are hundreds of metres apart, and material only when two candidates are nearly tied. Force points from a neighbouring zone into the same projection and the error grows quickly with distance from the meridian; two zones out it is percent-level, and the ranking starts to flip, which is the failure that matters. Haversine on a sphere carries its own residual — the Earth is an ellipsoid, so great-circle distances computed against a mean radius are off by up to roughly 0.3%, systematically by latitude.

If every query shares one origin — a single depot, one emergency dispatch centre — a third option beats both: reproject into an azimuthal equidistant projection centred on that origin. Distances measured from the centre point are then correct by construction, at any range, with no zone boundaries to straddle. It is the right frame for one-to-many distance work and the wrong one for many-to-many, since only distances from the centre are preserved.

from pyproj import CRS, Geod

# Distances from one dispatch centre, exact at any range
dispatch_lon, dispatch_lat = 7.6869, 45.0703
aeqd = CRS.from_proj4(
    f"+proj=aeqd +lat_0={dispatch_lat} +lon_0={dispatch_lon} +datum=WGS84 +units=m +no_defs"
)
incidents_aeqd = incidents.to_crs(aeqd)

# Spot-check any planar result against the true ellipsoidal distance
geod = Geod(ellps="WGS84")
_, _, true_m = geod.inv(dispatch_lon, dispatch_lat, 7.7412, 45.1188)
print(f"geodesic: {true_m:,.1f} m")   # geodesic: 7,468.2 m

Production Export & Integration

Serialize matches back to a spatial format so downstream services can consume them. GeoJSON is the lingua franca for web maps; strip everything but the fields the client needs to keep the payload small, and always emit WGS84 (EPSG:4326) for browser mapping libraries even though the analysis ran in a metric CRS.

import geopandas as gpd

# Build a nearest-neighbor edge layer and export for the web in EPSG:4326
edges = gpd.GeoDataFrame(
    {
        "sensor_id": sensors_clean["sensor_id"].to_numpy(),
        "nearest_sensor_id": sensors_clean["sensor_id"].to_numpy()[nearest_other],
        "gap_m": gap_m.round(2),
    },
    geometry=sensors_clean.geometry.to_numpy(),
    crs=sensors_clean.crs,
)
edges.to_crs(epsg=4326).to_file("nearest_edges.geojson", driver="GeoJSON")

For repeated queries against a stable dataset, persist the built index instead of rebuilding it on every request. joblib.dump serializes a cKDTree cleanly; a FastAPI worker can load it once at startup and serve k-NN lookups with sub-millisecond query latency. When the point set lives in a database, push the work server-side — PostGIS answers ordered nearest-neighbor queries with the <-> KNN-GiST operator, covered in PostGIS Integration with Python, which avoids shipping the whole table into memory at all.

import joblib

# Cache the index for a long-lived service; reload with joblib.load at startup
joblib.dump(sensor_tree, "sensor_tree.joblib")

A cached tree brings its own constraint: a cKDTree is immutable. There is no insert and no delete, so a layer that gains rows — new sensors commissioned, deliveries added through the day — needs a rebuild, and rebuilding a multi-million-point tree inside a request handler is not an option. The standard workaround is a two-tree pattern: keep the large, stable base tree in memory, accumulate additions in a small delta tree, query both, and take the closer of the two results. Rebuild the base tree on a schedule when the delta grows past a few percent of it.

import numpy as np

def nearest_across(base_tree, base_xy, delta_xy, query_xy):
    """Query a large static tree plus a small delta tree, keeping the better hit."""
    base_d, base_i = base_tree.query(query_xy, k=1, workers=-1)
    if len(delta_xy) == 0:
        return base_d, base_i
    delta_d, delta_i = cKDTree(delta_xy).query(query_xy, k=1, workers=-1)
    use_delta = delta_d < base_d
    # Offset delta positions so a single index space addresses both arrays
    return (
        np.where(use_delta, delta_d, base_d),
        np.where(use_delta, delta_i + len(base_xy), base_i),
    )

Two caveats on the persisted artefact. Pickle-based caches are not a stable format across SciPy versions — treat sensor_tree.joblib as a cache to be regenerated, never as an archive, and store the projected coordinate array beside it so a rebuild is one call away. And the cached tree only means anything alongside the CRS it was built in; write the EPSG code into the filename or a companion JSON, because nothing in the pickle records that the numbers are ETRS89 / UTM 32N metres rather than degrees.

Production checklist

Platform Edge Cases & Debugging

Most nearest-neighbor failures are silent — wrong answers, not exceptions — so the debugging discipline is to assert the invariants rather than eyeball the output.

Frequently Asked Questions

Should I use cKDTree or GeoPandas sjoin_nearest? Use sjoin_nearest by default — it handles any geometry, keeps both layers' attributes, and returns the distance. Drop to cKDTree when the data is millions of points and raw speed dominates. The sjoin_nearest vs cKDTree benchmark times both on identical data.

When do I need a BallTree instead of a KD-tree? When you must query directly in latitude/longitude across a large extent where no single projection stays accurate. BallTree with metric="haversine" measures great-circle distance on the sphere; a KD-tree cannot, because it only does Euclidean distance on flat coordinates.

How do I run nearest neighbor on polygons or lines? Reduce each geometry to a representative point with representative_point() (guaranteed on-surface, unlike centroid), index those points, then map the result indices back to the original features.

Is cKDTree different from KDTree in SciPy? They share an API; cKDTree is the compiled implementation and is what you should use. In current SciPy the two are effectively unified, but reach for cKDTree by name in performance-sensitive code.

How large a dataset before I switch to approximate search? Exact cKDTree stays fast into the low tens of millions of points on a workstation. Past that, or when query latency budgets are tight, move to pynndescent or faiss and accept a small recall trade-off.

Straight-line nearest or nearest along the road network? They answer different questions and often disagree by a factor of two — a depot across a river or a rail cutting is metres away in Euclidean terms and kilometres away by road. Use Euclidean nearest for deduplication, sensor pairing, and any screening step where you will refine later; use network distance whenever the result implies travel, which is what Street Network Analysis with OSMnx exists for. A common production compromise is a KD-tree shortlist of the ten closest candidates, then network routing on just those ten.

How do I get k nearest neighbours out of a spatial join instead? You mostly cannot, cleanly. sjoin_nearest returns the single nearest match (plus every tie), and coaxing k out of it means repeated joins with growing max_distance filters. That asymmetry is the strongest practical argument for the KD-tree route whenever k is greater than one — query(points, k=5) is one call and returns a rectangular array you can reshape into tidy long form.

Should the tree be built over the larger layer or the smaller one? Build over the layer you are matching to, regardless of size, because that is the layer whose rows the returned indices refer to. If both layers are large and you need the relationship in both directions, build both trees once and query each against the other rather than rebuilding per direction — construction is the expensive half.

Does deduplicating co-located points change my results? Yes, and usually for the better. Exactly duplicated coordinates unbalance the tree and make "nearest" ambiguous at zero distance, so a downstream join can attach the wrong attributes. Collapse duplicates to one row with an aggregated attribute set before building, and keep a count column if the multiplicity matters.

Why is my k-nearest result slower than a single nearest query by more than k times? Because the pruning radius grows with k: the search must keep the k best candidates so far, and the hypersphere it cannot prune past is the distance to the k-th best, not the first. Larger k crosses more splitting planes and touches more leaves. In practice the cost rises sub-linearly but noticeably, so ask for the k you need rather than a generous round number.

Related