K-Nearest-Neighbor Queries with ST_DWithin

Finding the $k$ closest features per probe point in DuckDB Spatial has no native operator, so the durable pattern is a LATERAL join whose subquery orders candidates by ST_Distance under an ST_DWithin radius guard that keeps the R-tree engaged — a technique this walkthrough builds up from a naive cross-join baseline, extends with a QUALIFY ROW_NUMBER() alternative and an expanding-radius fallback for sparse regions, and grounds in the units and correctness caveats that decide whether the answer is even right. It sits under Spatial Range and Nearest-Neighbor Search, which frames how bounded range search and unbounded kNN diverge; here the focus is narrow — making per-probe kNN both correct and index-driven.

Root-Cause Analysis #

kNN in DuckDB fails along four axes, each with a distinct fix, because unlike a fixed-radius range query there is no single window that returns exactly $k$ rows:

• No distance-order index descent. DuckDB has no <-> operator and no LIMIT k R-tree traversal, so a plain ORDER BY ST_Distance(...) LIMIT k computes an exact distance for every candidate before sorting. Against $M$ candidates per probe that is $O(N \times M \log M)$ — the query returns correct results but scales quadratically. The fix is to bound the candidate set with an ST_DWithin guard so the R-tree index prunes before the sort.
• Guard radius too tight. A radius that fits dense regions returns fewer than $k$ neighbours in sparse ones. The result is silently short, and the fix is an expanding retry, not a blanket-wide radius that penalizes every dense probe.
• Unit mismatch on the radius. ST_DWithin(a, b, r) measures r in coordinate units; a metric radius on EPSG:4326 data — or a degree radius on projected data — makes the guard meaningless. Both inputs must share a CRS, per CRS mapping and transformations.
• Non-deterministic ties. When candidates share the $k$-th distance, LIMIT and ROW_NUMBER pick an arbitrary subset unless the ordering carries a stable tiebreak.

Deterministic Configuration #

Load the extension, cap memory below total RAM, and ensure an R-tree exists on the candidate layer — the ST_DWithin guard is only fast if there is an index for it to drive.

INSTALL spatial;
LOAD spatial;

-- Cap resident memory below total RAM. A well-guarded kNN needs little; a low ceiling
-- makes an accidental unbounded sort spill or error early instead of exhausting RAM.
SET memory_limit = '4GB';

-- Match physical cores; R-tree probing is bandwidth-bound, so oversubscription adds
-- contention rather than throughput.
SET threads = 8;

-- Parallelize the index scan and drop order guarantees; re-impose order on the result.
SET preserve_insertion_order = false;

-- R-tree on the candidate layer so the ST_DWithin guard is index-servable, not a filter.
CREATE INDEX idx_shops_geom ON shops USING RTREE (geom);

Confirm the preconditions before trusting any result:

• Probe and candidate layers share one CRS, with the radius expressed in that CRS’s units.
• An R-tree index exists on the candidate layer’s bare GEOMETRY column.
• The guard radius is wide enough to contain the $k$-th neighbour for the sparsest probe you expect.

Optimized Execution Pattern #

The naive kNN pairs every probe with every candidate and sorts the full product per probe — correct, but it never touches the index:

-- BEFORE: cross join, exact distance on every pair, per-probe sort of all candidates.
SELECT incident_id, shop_id, metres
FROM (
    SELECT i.id AS incident_id, s.id AS shop_id,
           ST_Distance(i.geom, s.geom) AS metres,
           ROW_NUMBER() OVER (PARTITION BY i.id
                              ORDER BY ST_Distance(i.geom, s.geom)) AS rk
    FROM incidents i
    CROSS JOIN shops s          -- N x M pairs materialized before ranking
) ranked
WHERE rk <= 3;

The indexed rewrite bounds each probe’s candidate set with an ST_Expand bbox window and an ST_DWithin radius guard, so the R-tree hands the ORDER BY ... LIMIT k only a handful of rows:

-- AFTER: LATERAL per-probe subquery; the && / ST_DWithin guard keeps the R-tree engaged.
SELECT i.id AS incident_id, n.id AS shop_id, n.metres
FROM incidents i
CROSS JOIN LATERAL (
    SELECT s.id, ST_Distance(i.geom, s.geom) AS metres
    FROM shops s
    WHERE s.geom && ST_Expand(i.geom, 2000.0)  -- index-servable MBR window
      AND ST_DWithin(i.geom, s.geom, 2000.0)    -- exact radius bound on survivors
    ORDER BY metres, s.id                        -- s.id breaks distance ties deterministically
    LIMIT 3
) AS n;

The behavioural change is entirely in the candidate bound. The CROSS JOIN sorted all $M$ shops per incident; the LATERAL subquery’s s.geom && ST_Expand(i.geom, 2000.0) is index-servable, so the R-tree returns only shops whose envelope falls within 2 km, ST_DWithin trims the square window’s corners to a true circle, and the sort runs over that small set. The s.id in the ORDER BY makes ties reproducible.

The QUALIFY ROW_NUMBER() form expresses the same result without a correlated subquery — one bounded join, ranked in a single pass. QUALIFY filters on the window result directly, avoiding the wrapping subquery the naive form needs:

-- Ranking alternative: bounded join once, ROW_NUMBER per probe, QUALIFY keeps top k.
SELECT i.id AS incident_id, s.id AS shop_id,
       ST_Distance(i.geom, s.geom) AS metres
FROM incidents i
JOIN shops s
  ON s.geom && ST_Expand(i.geom, 2000.0)
WHERE ST_DWithin(i.geom, s.geom, 2000.0)
QUALIFY ROW_NUMBER() OVER (PARTITION BY i.id
                           ORDER BY ST_Distance(i.geom, s.geom), s.id) <= 3;

Prefer LATERAL when the guard radius is tight and LIMIT k lets each probe stop early; prefer QUALIFY when probes are numerous and share a candidate window the planner can vectorize as one ranked scan.

Expanding-radius fallback #

When a probe in a sparse region returns fewer than $k$ neighbours, widen only that probe’s window rather than inflating the radius globally. Detect the shortfall and re-run the short probes at a larger radius:

-- Identify probes that came up short under the base guard radius.
WITH knn AS (
    SELECT i.id AS incident_id, s.id AS shop_id,
           ST_Distance(i.geom, s.geom) AS metres
    FROM incidents i
    JOIN shops s ON s.geom && ST_Expand(i.geom, 2000.0)
    WHERE ST_DWithin(i.geom, s.geom, 2000.0)
    QUALIFY ROW_NUMBER() OVER (PARTITION BY i.id
                               ORDER BY ST_Distance(i.geom, s.geom), s.id) <= 3
)
SELECT incident_id, COUNT(*) AS found
FROM knn
GROUP BY incident_id
HAVING COUNT(*) < 3;   -- these probes need a wider radius

Re-probe the short list at, say, 8000 m, then union the widened results with the probes that already had their full $k$. Doubling the radius each round bounds the retries to $O(\log)$ passes while keeping every dense probe on the tight, fast window.

Correctness and units caveats #

Two subtleties decide whether an indexed kNN is actually correct, not just fast. First, the guard radius and ST_Distance share the geometry’s coordinate units, and both must be the units you think they are. On a metre-based projection ST_DWithin(i.geom, s.geom, 2000.0) means 2 km and metres is a true ground distance; on EPSG:4326 the same call means 2000 degrees and the returned distance is a degree magnitude that stretches roughly $\cos(\text{latitude})$ east–west, so a “nearest” ranking computed in degrees can reorder neighbours near the poles. Project both layers to a metric CRS before ranking whenever the ordering must reflect real distance:

-- Rank in metres regardless of source CRS: reproject both inputs into the guard's units.
SELECT DISTINCT ST_SRID(geom) AS srid FROM shops;  -- expect the projected SRID, not 4326

Second, the ST_Expand window is a square and ST_DWithin is a circle, so the && stage over-selects the corners; that is correct — it only ever admits extra candidates the exact ST_DWithin then rejects, never drops a true neighbour. The danger is the reverse: an ST_Expand window narrower than the guard radius would silently discard valid neighbours before the exact test runs, so always expand by at least the radius you pass to ST_DWithin.

Diagnostic Queries & Plan Validation #

Confirm the R-tree drives the guard rather than a sequential scan:

EXPLAIN ANALYZE
SELECT i.id, s.id
FROM incidents i
JOIN shops s ON s.geom && ST_Expand(i.geom, 2000.0)
WHERE ST_DWithin(i.geom, s.geom, 2000.0)
QUALIFY ROW_NUMBER() OVER (PARTITION BY i.id
                           ORDER BY ST_Distance(i.geom, s.geom), s.id) <= 3;

Read the plan and check three things:

• Scan type. An RTREE_INDEX_SCAN of shops is correct; a SEQ_SCAN under a NESTED_LOOP_JOIN means the && predicate never reached the index — verify the indexed column is bare and the units match.
• Candidate depth. The row count entering the ORDER BY/window stage should be a low multiple of $k$ per probe, not the full candidate layer. A depth near $M$ per probe means the guard radius is too wide or missing.
• Completeness. Assert every probe returns exactly $k$ rows; a partition returning fewer is a sparse-region shortfall, not an error.

-- One-liner completeness monitor: any probe below k needs the expanding fallback.
SELECT incident_id, COUNT(*) AS n
FROM knn_result GROUP BY incident_id HAVING COUNT(*) < 3;

-- Confirm the index is registered before trusting the plan.
SELECT index_name, index_type FROM duckdb_indexes() WHERE table_name = 'shops';

The distance computation itself is the same kernel used to build full pairwise grids; the bounded-block strategies that keep those from OOM are in calculating distance matrices with SQL.

Geometry Validation & Fallback Routing #

Invalid candidate geometry makes ST_Distance return misleading values or throw mid-sort, so gate the layer before indexing it:

-- Quarantine invalid geometry before it poisons the distance ranking.
SELECT id, geom FROM shops WHERE NOT ST_IsValid(geom);

-- Repair, then rebuild the R-tree so its envelopes match the corrected geometry.
UPDATE shops SET geom = ST_MakeValid(geom) WHERE NOT ST_IsValid(geom);
DROP INDEX IF EXISTS idx_shops_geom;
CREATE INDEX idx_shops_geom ON shops USING RTREE (geom);

For a probe count large enough to strain memory, run the LATERAL kNN in chunks — partition the probe layer by a coarse grid cell or an id range and process one chunk per pass, so each LATERAL invocation holds only its slice of candidate sets. Chunked execution keeps peak memory flat regardless of probe count, the same bounded-pass discipline the spatial joins and proximity filters reference applies to large joins.

Related #

See also

• Spatial joins and proximity filters — the && two-stage pruning and chunked-join discipline the kNN guard reuses.
• R-tree spatial indexing internals — how the index the ST_DWithin guard drives is built and probed.
• Calculating distance matrices with SQL — the bounded pairwise-distance kernel behind per-probe ranking.

Up: Spatial Range and Nearest-Neighbor Search · Modern Spatial SQL Query Patterns