Spatial Joins & Proximity Filters

Spatial joins and proximity filters are the primary computational bottleneck in analytical geospatial pipelines: a join predicate that reads cleanly as ST_Intersects or ST_DWithin quietly expands into an unbounded cross-product, forcing $O(N \times M)$ topology evaluations across two relations. This page, part of the wider Modern Spatial SQL Query Patterns reference, isolates the operator-level workflow: how to phrase a spatial join so DuckDB routes it through an R-tree instead of a nested loop, how to apply distance predicates without silent unit errors, and how to read the execution plan so a regression surfaces before it reaches production. Modern columnar engines mitigate the worst case through vectorized execution, in-memory bounding-box indexing, and predicate pushdown — but only when the query syntax aligns with the optimizer’s spatial routing heuristics. Everything below assumes you want deterministic, sub-second joins at scale rather than a query that merely returns the right answer eventually.

The lower-level mechanics that make these patterns work — how the engine stores geometry and builds its index — are covered in the R-tree index internals and ST_Geometry vs WKB references. For the single most common join shape, see the dedicated walkthrough on optimizing point-in-polygon queries in DuckDB.

Runtime Configuration & Memory Guardrails #

Spatial indexing and topology evaluation are highly sensitive to memory allocation and thread contention. DuckDB constructs its bounding-box index in memory and evaluates geometry predicates in vectorized batches, so a misconfigured session silently spills to disk or thrashes threads, degrading throughput by 3–5× without raising any error. Apply these session-level guardrails before executing any join workload:

INSTALL spatial;
LOAD spatial;

-- Cap resident memory below total RAM: a spatial join holds both inputs plus
-- the index plus the result hash region simultaneously, so leave headroom.
SET memory_limit = '8GB';

-- Match physical cores; hyperthreading degrades index-build parallelism because
-- the R-tree construction phase is memory-bandwidth bound, not compute bound.
SET threads = 4;

-- Drop order preservation to unlock parallel index construction. This trades
-- away ORDER BY guarantees on unindexed columns for a faster build.
SET preserve_insertion_order = false;

-- Fast local NVMe for spill: if the join overflows memory_limit, a slow temp
-- dir turns graceful spill into I/O thrash.
SET temp_directory = '/var/lib/duckdb/spill';

Trade-off Analysis:

• preserve_insertion_order = false unlocks parallel index builds but invalidates ORDER BY guarantees on unindexed columns. Compensate by materializing results into an explicitly sorted output table when downstream consumers depend on order.
• memory_limit must exceed the combined uncompressed geometry footprint of both join inputs plus the index. Monitor live allocation with SELECT * FROM duckdb_memory();. If EXPLAIN ANALYZE reports External Merge or a spill during index construction, raise the limit or partition the input by bounding-box bounds.
• threads should track physical cores. Whether the working set fits in memory at all depends on the in-memory vs disk storage trade-offs for your dataset; oversubscribing threads only raises peak memory once index build saturates bandwidth.

Primary Execution Patterns #

The canonical production join is point-in-polygon assignment, and the rule that governs every spatial join derives from it: never let the engine evaluate exact topology on every candidate pair. A naive JOIN ... ON ST_Within(point, polygon) bypasses the index and forces a full cross-product scan. Instead, implement a two-stage filter — a cheap bounding-box pre-filter that the index can serve, followed by precise topology evaluation only on the survivors.

The cheap bounding-box test eliminates the majority of pairs before the expensive exact-geometry check runs.

-- Two-stage spatial join: bbox pre-filter, then exact topology.
SELECT p.id, z.zone_id
FROM points p
JOIN zones z
  ON p.geom && z.geom            -- Stage 1: cheap bounding-box overlap (index-served)
  AND ST_Within(p.geom, z.geom); -- Stage 2: exact containment on survivors only

Performance Impact: Bounding-box pre-filtering typically eliminates 60–90% of candidate pairs before any topology check runs. The optimizer recognizes && in the join condition and leverages an R-tree index for the bounding-box stage when one is present — the structure and rebuild cost of that index is detailed in the R-tree index internals reference.

Materializing inputs to force index use #

DuckDB will not build an index over a geometry expression embedded in a join, and it cannot index a column that arrives mid-pipeline from a CTE that the planner chose not to materialize. When a join input is itself the output of a scan or transform, stage it into a temp table first so the index has a stable column to attach to:

-- Stage the larger relation so the planner can build a bbox index over it.
CREATE TEMP TABLE zones_mat AS SELECT zone_id, geom FROM zones;
CREATE INDEX zones_rtree ON zones_mat USING RTREE (geom);

SELECT p.id, z.zone_id
FROM points p
JOIN zones_mat z
  ON p.geom && z.geom
  AND ST_Within(p.geom, z.geom);

This pattern is also the foundation for the deeper tuning in optimizing point-in-polygon queries in DuckDB, where envelope strategies and index-utilization thresholds are tuned per workload.

Proximity Filters & Distance Predicates #

Proximity filters built on ST_DWithin require strict distance-unit handling and CRS awareness. DuckDB evaluates distances in the geometry’s native coordinate reference system, so an ST_DWithin(a, b, 500) against unprojected WGS84 (EPSG:4326) compares degrees, not metres — a query that runs without error and returns wrong matches. Always project to a metric CRS (a local UTM zone, or EPSG:3857 for coarse work) before applying a distance predicate; the mechanics of doing this safely are covered in the CRS transformation reference.

-- Pre-project the static reference table ONCE so ST_Transform does not re-run
-- per candidate pair inside the join.
CREATE TEMP TABLE facilities_m AS
SELECT facility_id, ST_Transform(geom, 'EPSG:4326', 'EPSG:32633') AS geom_m
FROM facilities;

WITH sensors_m AS (
    SELECT id, ST_Transform(geom, 'EPSG:4326', 'EPSG:32633') AS geom_m
    FROM sensors
)
SELECT s.id, f.facility_id
FROM sensors_m s
JOIN facilities_m f
  ON ST_DWithin(s.geom_m, f.geom_m, 500.0); -- 500 metres, now that the CRS is metric

Trade-off Analysis:

• Reprojection adds roughly 15–30% compute overhead per row but guarantees metric accuracy. For static reference tables, pre-project once into the target CRS (as above) to eliminate repeated ST_Transform calls inside the hot loop.
• ST_DWithin uses the same bounding-box routing as ST_Intersects when placed in ON, but the distance buffer expands each envelope by the search radius, inflating candidate-pair volume. For high-density inputs, pre-filter with a coarse grid (floor(x / cell), floor(y / cell)) so neighbours land in adjacent cells and the topology kernel sees far fewer pairs.

Vectorized distance evaluation batches calculations across SIMD lanes. When a proximity filter feeds an aggregation — counting facilities within range, summing weighted distances — keep the ST_DWithin predicate in the join condition and defer GROUP BY so intermediate geometry is never materialized twice. This is the same discipline described in the vectorized aggregations reference, and it underpins efficient distance-matrix construction.

Execution Plan Validation #

Run EXPLAIN ANALYZE on every join before it ships. A correctly routed plan shows the spatial predicate attached to an index scan, with actual row counts close to the estimates:

Diagnostic Boundaries:

• If the plan shows Cross Product or a plain Hash Join with no spatial predicate on the operator, the predicate is not being routed. Verify that both && and the topology check live in the ON clause, not WHERE — the optimizer only triggers spatial routing when the predicate is part of the join condition.
• If timing on the join operator exceeds ~500 ms for fewer than 500k rows, the index failed to build, usually due to memory pressure or an un-materialized input. Stage inputs into temp tables (see above) and re-check.
• If Actual Rows deviates more than ~20% from Estimated Rows on the join, statistics are stale; run ANALYZE <table> so the planner sizes the build side correctly.

Performance Trade-offs #

The patterns above are not universally optimal — each has a regime where it wins and a regime where it costs more than it saves. Choose deliberately:

The decisive variable is reuse. A single ad-hoc join rarely justifies an index build or a pre-projection pass; a reference table queried thousands of times almost always does. For workloads driven from Python where the same reference geometry is reused across requests, the staging cost amortizes especially well — see async execution patterns for keeping that warmed state alive across queries.

Edge Cases & Anti-Patterns #

Most spatial-join performance incidents trace back to a small set of repeatable mistakes. Each below pairs the symptom with the minimal fix.

Predicate in WHERE instead of ON. The optimizer routes spatial predicates only from the join condition. Moving the filter to WHERE forces a full cross-product first:

-- ANTI-PATTERN: cross product materializes before the filter runs.
SELECT p.id, z.zone_id
FROM points p, zones z
WHERE ST_Within(p.geom, z.geom);

-- FIX: predicate in ON, with a bbox pre-filter the index can serve.
SELECT p.id, z.zone_id
FROM points p
JOIN zones z ON p.geom && z.geom AND ST_Within(p.geom, z.geom);

Distance predicate against unprojected coordinates. ST_DWithin(a, b, 500) on EPSG:4326 compares degrees; 500 degrees is the whole globe. Project both sides to a metric CRS first (see the proximity section above), or the join returns plausible-looking but wrong matches.

Invalid geometry poisoning the topology kernel. A self-intersecting polygon or an unclosed ring makes ST_Within return inconsistent results or error mid-scan. Guard inputs and repair before joining:

SELECT zone_id, ST_MakeValid(geom) AS geom
FROM zones
WHERE NOT ST_IsValid(geom);

Function-wrapped join column. Wrapping the indexed column in any function (ST_Within(ST_Buffer(p.geom, 0), z.geom)) discards index eligibility, because the index is on p.geom, not on the buffered expression. Pre-compute the buffer into a materialized column instead.

Mixed-CRS inputs joined without checking. Two tables in different reference systems join without error and return near-empty or nonsensical results. Confirm both sides share a CRS — DuckDB stores no inline SRID, so this must be tracked out-of-band per the CRS transformation reference.

When clustering points before a join (for example, grouping sensors into service areas), prefer a density-aware grouping such as ST_ClusterDBSCAN over a raw self-join — it avoids the quadratic blowup a self-join would otherwise incur.

Which side becomes the build side, and why it matters here #

For an ordinary equi-join, the smaller relation becomes the hash build side and row count is a fine proxy for “smaller”. For a spatial join it is not: the cost that matters is total vertex count, and the two measures routinely disagree by an order of magnitude. Two million points are a lighter build side than fifty thousand municipal boundaries, and a planner working from row counts alone will pick the wrong one.

The two columns point at opposite answers, and only one of them is about memory.

When a plan puts the heavy side on the build, the visible symptom is a join that spills during its build phase rather than its probe phase — duckdb_temporary_files() grows before any output row appears. The fix is usually to give the optimizer better information rather than to fight it: reduce the polygon side first with the most selective attribute filter available, or simplify the geometry to the tolerance the answer actually needs.

-- Measure before you theorise. Row counts are free; vertex counts are what
-- the join will actually carry.
SELECT 'incidents'  AS side, count(*) AS rows, sum(ST_NPoints(geom)) AS vertices FROM incidents
UNION ALL
SELECT 'boundaries',        count(*),          sum(ST_NPoints(geom))            FROM boundaries;

Query Regression Analysis #

Spatial join performance degrades silently when data distributions shift or statistics go stale. Capture a baseline plan and diff against it in CI so a regression is caught at merge time, not in production. The cheapest harness wraps EXPLAIN ANALYZE through the Python API:

import duckdb

def capture_spatial_plan(con: duckdb.DuckDBPyConnection, query: str) -> str:
    """Execute EXPLAIN ANALYZE and return the plan text for diffing."""
    plan_df = con.sql(f"EXPLAIN ANALYZE {query}").fetchdf()
    return "\n".join(plan_df.iloc[:, 0].tolist())

con = duckdb.connect(":memory:")
con.execute("INSTALL spatial; LOAD spatial;")
con.execute("SET memory_limit = '8GB'; SET threads = 4;")

# Capture once; commit the plan as a fixture, then re-run and diff on every change.
baseline_plan = capture_spatial_plan(con, """
    SELECT p.id, z.zone_id
    FROM points p
    JOIN zones z ON p.geom && z.geom AND ST_Within(p.geom, z.geom)
""")

Diagnostic Boundaries for Regression:

• Row-estimate drift: if Actual Rows deviates more than 20% from Estimated Rows, run ANALYZE <table> or rebuild the index before trusting the plan.
• Operator substitution: watch for a Spatial Join downgrading to Hash Join or Nested Loop between runs — that means the optimizer lost spatial-predicate visibility, usually because a CTE stopped being materialized or a column gained a function wrapper.
• Window-function integration: when ranking nearest neighbours or computing spatial moving averages, combine ST_DWithin with partitioned windowing and keep a consistent ORDER BY inside OVER() to avoid a full sort. The partitioning strategies that preserve spatial locality are detailed in the window functions for geospatial reference.

For pipelines that hand join results back to GeoPandas, capture the plan on the DuckDB side before the DuckDB-to-GeoPandas sync boundary — a regression hidden behind the conversion is far harder to attribute.

Frequently Asked Questions #

Why does my spatial join return more rows than either input? #

Because a spatial join is many-to-many by nature: one point can fall inside several overlapping polygons, and one polygon can contain many points. If the polygons are supposed to be a partition and are not — overlapping administrative boundaries, duplicated features from a bad import — the fan-out is a data defect rather than a query defect. Check with a count of matches per input row before assuming the join is wrong.

Is ST_Intersects or ST_Contains faster? #

They cost about the same, because both run the same underlying topology evaluation. Choose on meaning, not speed: ST_Intersects is true for any shared point including a touching edge, while ST_Contains requires the interior. Getting this wrong shows up as boundary features appearing in two zones or in none, which is a correctness bug that no amount of tuning fixes.

Do I need both the bounding-box operator and the exact predicate? #

Yes, and in that order. The && operator alone gives false positives — envelopes overlap for shapes that do not — so the exact predicate is what makes the answer correct. The exact predicate alone gives the right answer at full cost, because there is nothing an index can serve. The pair is what makes it both correct and fast.

Why is ST_DWithin returning nothing? #

Nearly always a unit mismatch. On geographic coordinates the radius is in degrees, so a threshold of 50 meaning metres actually asks for fifty degrees, and a threshold of 0.0005 meaning degrees asks for roughly fifty metres. Confirm the coordinate range before adjusting the number: if the values fall inside plus or minus 180, the layer is still in degrees and the fix is a reprojection, not a different radius.

Should I index both sides of the join? #

Usually only the probe side benefits, because the build side is scanned in full to populate the hash table regardless. Indexing both costs build time and storage for no return on that particular query. Index the table that gets probed by a bounding-box window — typically the larger, more static layer — and leave the other alone.

How do I join points to the single nearest polygon rather than all matches? #

Rank within the match set rather than trying to express it as a join condition. Join with the two-stage predicate as usual, compute ST_Distance on the survivors, and use QUALIFY row_number() OVER (PARTITION BY point_id ORDER BY distance) = 1. Attempting to push “nearest” into the join predicate produces either a correlated subquery per row or a cross product, both of which are far slower.

Related #

See also:

• Tuning ST_DWithin radius and join order — deriving the radius from density, and choosing the build side.
• Optimizing point-in-polygon queries in DuckDB — the deep dive on the most common join shape.
• Vectorized aggregations — chaining proximity filters into aggregates without re-materializing geometry.
• Window functions for geospatial — nearest-neighbour ranking and spatial windowing.
• R-tree index internals — how the bounding-box index that powers && is built and rebuilt.
• CRS transformation — projecting to a metric CRS before distance predicates.

Up: Modern Spatial SQL Query Patterns

---

External Reference Standards:

• For CRS transformation accuracy and projection parameters, consult the PROJ Documentation and the EPSG Geodetic Parameter Dataset.
• For DuckDB spatial extension capabilities and configuration, review the official DuckDB Spatial Extension Docs.