Optimizing Point-in-Polygon Queries in DuckDB

Point-in-polygon (PIP) assignment is the single most common spatial join shape, and it degrades the moment a query lets DuckDB Spatial evaluate exact containment on every candidate pair instead of pruning with a bounding-box index first. This walkthrough sits under Spatial Joins & Proximity Filters and isolates the exact bottleneck — ST_Within / ST_Contains over millions of points against thousands of polygons — then gives deterministic configuration, an indexed two-stage rewrite, plan-validation queries, and fallback routing for invalid geometry and out-of-memory (OOM) conditions. The goal is sub-second resolution that stays sub-second as the input grows, not a query that merely returns the right answer eventually.

Root-Cause Analysis of Execution Degradation #

DuckDB’s vectorized engine stores geometry internally as WKB and decodes it per batch; the layout trade-offs behind that representation are detailed in the ST_Geometry vs WKB reference. Without a bounding-box index, every point is tested against every polygon, forcing $O(N \times M)$ topology evaluations. PIP queries in production fail along four distinct axes, and each has a different fix:

• Missing bounding-box pre-filter. An exact ST_Contains with no cheap envelope test in front of it walks every polygon vertex for every point. The fix is an R-tree index plus the && overlap operator so the optimizer can prune candidates before topology runs.
• Predicate placement that defeats index routing. Putting the spatial predicate only in WHERE — or wrapping the indexed column in a function — prevents the planner from recognising an index-servable bounding-box stage. The MBR test must appear as && in the ON clause against the bare indexed GEOMETRY column.
• Unbounded memory during the join. DuckDB materialises decoded geometry before applying predicates, so a large unpruned join spikes the heap and spills. Whether the working set fits at all is governed by the in-memory vs disk storage boundaries for your dataset.
• CRS / unit mismatch. Mixing degree-space points with a projected polygon layer makes envelopes miss and any tolerance value meaningless. Align both inputs first; the rules are in the CRS mapping & transformations reference.

Deterministic Configuration & Memory Boundaries #

Enforce strict memory boundaries and load the extension before rewriting any query logic. DuckDB Spatial does not auto-index geometry columns, so the index step is mandatory, not optional. Confirm the prerequisites first:

• spatial extension installed and loaded in the session
• Points and polygons share one CRS (degree-space EPSG:4326 throughout, or both projected to a metric CRS)
• memory_limit set below total RAM with headroom for both inputs plus the index
• A fast local temp_directory available for graceful spill

INSTALL spatial;
LOAD spatial;

-- Cap resident memory below total RAM: the join holds both inputs, the
-- index, and the result region at once, so leave headroom for spill-free runs.
SET memory_limit = '4GB';

-- Match physical cores; R-tree build is memory-bandwidth bound, so
-- oversubscribing threads raises peak memory without speeding the build.
SET threads = 8;

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

-- Quiet the progress bar so EXPLAIN ANALYZE timings are not skewed by render overhead.
SET enable_progress_bar = false;

For datasets beyond 10M rows, partition by a spatial grid or materialise the larger relation into a temp table so the index has a stable column to attach to and I/O aligns with columnar cache boundaries.

Spatial Pre-Filtering & Index Materialization #

Create an R-tree index on the boundary table to enable index-assisted bounding-box filtering. A separate envelope column is not required: the R-tree on the native GEOMETRY column serves the minimum-bounding-rectangle (MBR) test directly.

-- Create an R-tree index for fast bounding-box lookups on the polygon layer.
CREATE INDEX idx_boundaries_geom ON boundaries USING RTREE (geom);

-- Verify the index exists and is the type you expect.
SELECT index_name, table_name, index_type
FROM duckdb_indexes()
WHERE table_name = 'boundaries';

Index materialisation caps memory residency, prevents OOM during large joins, and lets the planner push the cheap MBR intersection ahead of exact geometry evaluation. The internal node structure and rebuild cost of this index are covered in the R-tree index internals reference.

Optimized Execution Pattern #

Take a baseline PIP query assigning 5M GPS points to 10K administrative boundaries. The naive form names the predicate once and lets the engine choose a nested loop:

-- Naive: nested loop join, full geometry evaluation, no index utilisation.
SELECT p.id, p.lat, p.lon, b.name
FROM points p
JOIN boundaries b ON ST_Contains(b.geom, ST_Point(p.lon, p.lat));

The optimised form splits the predicate into two stages: a cheap MBR overlap the R-tree can serve, then exact containment only on the survivors.

-- Optimized: bbox filter in ON clause triggers an R-tree scan;
-- exact containment runs only on the surviving candidate set.
SELECT p.id, p.lon, p.lat, b.name
FROM points p
JOIN boundaries b
  ON b.geom && ST_Point(p.lon, p.lat)            -- Stage 1: MBR overlap (R-tree index)
WHERE ST_Contains(b.geom, ST_Point(p.lon, p.lat)); -- Stage 2: exact topology on survivors

The behavioural change is entirely in predicate placement: && in the ON clause is index-servable, so the planner replaces the nested loop with an index scan that hands a far smaller candidate set to ST_Contains. In practice the MBR stage prunes 85–99% of pairs before any topology check runs. Keep the indexed column (b.geom) bare on the left of &&; wrapping it in a function would hide it from the index.

Diagnostic Queries & Plan Validation #

Validate the rewrite with EXPLAIN ANALYZE rather than wall-clock guesswork:

EXPLAIN ANALYZE
SELECT p.id, b.name
FROM points p
JOIN boundaries b
  ON b.geom && ST_Point(p.lon, p.lat)
WHERE ST_Contains(b.geom, ST_Point(p.lon, p.lat));

Read the plan top-down and check three things:

• Join node type. You want an index/hash-style scan driven by the && predicate, not NESTED_LOOP_JOIN over the full polygon set. A nested loop is the signature of a defeated index.
• Estimated vs actual rows. Large drift between the planner’s estimate and the actual count at the join node means stale statistics or a non-selective MBR stage; re-check CRS alignment and index presence.
• Operator timing. On a pruned candidate set under ~500K rows, the exact-topology stage and any downstream vectorized aggregations should land within a couple hundred milliseconds. If the topology stage dominates total time, the MBR stage is not pruning.

One-liner monitors for live runs:

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

-- Watch peak allocation and spill while the join runs.
SELECT tag, memory_usage_bytes, temporary_storage_bytes
FROM duckdb_memory() ORDER BY memory_usage_bytes DESC;

If timing exceeds baseline, confirm && placement in the ON clause and verify the index appears in EXPLAIN output before changing anything else.

Geometry Validation & Fallback Routing #

Invalid polygons cause silent wrong answers or topology errors, so gate the boundary layer before indexing it:

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

Repair, then rebuild the index so the R-tree reflects the corrected envelopes:

-- Repair invalid polygons, then re-index so MBRs match the fixed geometry.
UPDATE boundaries SET geom = ST_MakeValid(geom) WHERE NOT ST_IsValid(geom);

DROP INDEX IF EXISTS idx_boundaries_geom;
CREATE INDEX idx_boundaries_geom ON boundaries USING RTREE (geom);

For GPS noise near a boundary edge, a strict containment test drops points that should match. Fall back to a small tolerance with ST_DWithin, keeping the && pre-filter in front so the index still drives the scan:

-- Tolerance-based fallback for edge noise; ~11m at the equator in EPSG:4326 degrees.
SELECT p.id, b.name
FROM points p
JOIN boundaries b
  ON b.geom && ST_Point(p.lon, p.lat)
WHERE ST_DWithin(b.geom, ST_Point(p.lon, p.lat), 0.0001);

When memory is still breached despite pruning, stop evaluating PIP per raw point. Aggregate points to a coarse grid cell first, then assign one representative point per cell — turning millions of exact checks into thousands:

-- Grid pre-aggregation: count points per cell, assign cell centroids to polygons.
WITH grid_counts AS (
    SELECT
        floor(lon / 0.01) AS cell_x,
        floor(lat / 0.01) AS cell_y,
        count(*)          AS point_count
    FROM points
    GROUP BY cell_x, cell_y
)
SELECT g.cell_x, g.cell_y, g.point_count, b.name
FROM grid_counts g
JOIN boundaries b
    ON b.geom && ST_Point(g.cell_x * 0.01 + 0.005, g.cell_y * 0.01 + 0.005)
WHERE ST_Contains(b.geom, ST_Point(g.cell_x * 0.01 + 0.005, g.cell_y * 0.01 + 0.005));

This reduces exact PIP evaluations by an order of magnitude while preserving analytical accuracy — the same count-per-region question, answered against cells instead of raw fixes. For row-bounded chunking patterns and N×N variants of this trade-off, see the sibling walkthrough on calculating distance matrices with SQL.

Related #

See also

• Calculating Distance Matrices with SQL — bounded two-stage filtering for the $O(N^2)$ pairwise case.
• R-tree index internals — how the bounding-box index this page relies on is built and maintained.

Up: Spatial Joins & Proximity Filters · Modern Spatial SQL Query Patterns