Vectorizing Shapely Predicates over DuckDB

A Python loop that tests one Shapely geometry at a time — for geom in gdf.geometry: if region.contains(geom) — is the slowest possible way to evaluate a spatial predicate over a large layer, and this guide, part of the Shapely integration reference, shows how to replace it with DuckDB’s vectorized ST_ kernels and R-tree, keeping Shapely only for the custom operations DuckDB genuinely lacks.

Root-Cause Analysis of the Python-Loop Bottleneck #

The per-row Shapely predicate loop is not merely a little slow; it is architecturally the wrong shape for the work, and it fails along four independent axes at once. Understanding which axis dominates tells you whether the fix is pushing the predicate into SQL or batching the Shapely call — the two distinct remedies below.

• The GIL serialises every comparison. CPython holds a single interpreter lock, so a contains loop over a million geometries runs on exactly one core no matter how many the machine has. DuckDB’s ST_ kernels release into native code and fan the same work across every thread, turning a single-core crawl into a parallel scan.
• Per-object overhead dwarfs the geometry math. Each iteration boxes a GEOS geometry into a CPython object, dispatches a bound method, and — unless you prepared it explicitly — rebuilds GEOS’s internal predicate index on the test geometry from scratch. For small polygons the interpreter overhead exceeds the actual point-in-polygon test many times over.
• There is no spatial index. A bare Python loop compares every candidate against every target, forcing $O(N \times M)$ exact topology evaluations. DuckDB routes the same predicate through an R-tree index and the && bounding-box operator, pruning the overwhelming majority of pairs before any exact test runs — the same two-stage strategy detailed in optimizing point-in-polygon queries.
• Materialisation pins the whole layer in the Python heap. To loop, you first pull every geometry into memory as a Python list of objects — each with the full CPython object header — so a layer that occupies a compact WKB buffer in DuckDB balloons several-fold as boxed objects before the first comparison even runs.

shapely.vectorized.contains and the Shapely 2.x ufuncs remove the interpreter overhead of axis two, and they are the right tool once data is already in NumPy — but they still run in one process with no persistent spatial index, so on a join between two large layers they remain $O(N \times M)$. The moment a predicate spans two sizeable relations, the index is what matters, and that lives in DuckDB.

Deterministic Configuration #

Pin the connection, load the extension, and build the R-tree on the layer that plays the target role in the predicate (the polygons a point must fall inside, the regions a geometry must intersect). DuckDB Spatial does not auto-index geometry, so this step is mandatory for the pushdown to beat the loop.

import duckdb
import shapely
import numpy as np
import pyarrow as pa

con = duckdb.connect(config={
    "threads": 8,                     # fan ST_ kernels across cores; the loop could use only one
    "memory_limit": "12GB",           # hard ceiling below physical RAM; leave headroom for WKB buffers
    "preserve_insertion_order": False,# unlock parallel index build and scan
})
con.execute("INSTALL spatial; LOAD spatial;")

# Index the target layer so && prunes candidates before exact topology runs.
con.execute("CREATE INDEX idx_regions_geom ON regions USING RTREE (geom);")

Confirm the preconditions before rewriting any loop:

• Both layers share one CRS, so bounding boxes actually overlap and any tolerance is meaningful.
• An R-tree exists on the target layer and appears in duckdb_indexes().
• The predicate you need has a native ST_ equivalent — check before assuming a Shapely round-trip is required.
• Geometry crosses the boundary as WKB (ST_AsWKB / ST_GeomFromWKB), never stringified WKT.

Optimized Execution Pattern #

Case 1 — predicate loop to a DuckDB join #

The predicate exists in DuckDB, so the entire loop collapses into one indexed SQL join. Nothing iterates in Python.

# BEFORE — single-core, no index, O(N × M): tests every geometry against the region.
import geopandas as gpd

regions = gpd.read_file("regions.gpkg")
points = gpd.read_file("points.gpkg")
region = regions.union_all()            # one big polygon

hits = []
for geom in points.geometry:            # GIL-bound loop, one core
    if region.contains(geom):           # rebuilds prepared predicate each call
        hits.append(geom)

-- AFTER — vectorized, multi-threaded, index-pruned. && drives the R-tree,
-- ST_Contains runs only on the surviving candidate pairs.
SELECT p.point_id, r.region_id
FROM points p
JOIN regions r
  ON r.geom && p.geom              -- Stage 1: MBR overlap served by the R-tree
 AND ST_Contains(r.geom, p.geom);  -- Stage 2: exact topology on survivors only

The behavioural change is total: the for loop’s $O(N \times M)$ single-core topology sweep becomes a two-stage join where the && predicate lets the planner prune 85–99% of pairs through the index before ST_Contains ever runs, across every thread. Keep the indexed column (r.geom) bare on the left of &&; wrapping it in a function hides it from the index. The join-side tuning for this shape is developed in spatial joins and proximity filters.

Case 2 — object loop to a batched Arrow round-trip #

Sometimes the operation genuinely is not in DuckDB — an oriented minimum bounding rectangle, a shared-paths decomposition, a bespoke simplification. Here the fix is not to abandon Shapely but to stop iterating objects: pull the geometry as a WKB column, parse the whole column into a NumPy array with shapely.from_wkb, and run the Shapely 2.x ufunc across the array. One call, native inner loop, no Python-level iteration.

# BEFORE — object loop: boxes every geometry, dispatches per element under the GIL.
wkb_list = [row[0] for row in con.execute(
    "SELECT ST_AsWKB(geom) FROM parcels").fetchall()]

geoms = [shapely.from_wkb(b) for b in wkb_list]     # N Python calls
mrr   = [g.minimum_rotated_rectangle for g in geoms] # N Python calls, one core

# AFTER — one bulk parse, one vectorized ufunc over a NumPy array, then bulk serialize back.
tbl = con.execute("SELECT parcel_id, ST_AsWKB(geom) AS wkb FROM parcels").fetch_arrow_table()

# from_wkb over a numpy array of bytes: a single vectorized parse, no per-row Python.
geoms = shapely.from_wkb(tbl.column("wkb").to_numpy(zero_copy_only=False))

# Shapely 2.x exposes GEOS as numpy ufuncs — the loop lives in C, not the interpreter.
mrr = shapely.minimum_rotated_rectangle(geoms)   # custom op DuckDB lacks, still vectorized

# Serialize the whole array back to WKB in one call and hand it to DuckDB.
out_wkb = shapely.to_wkb(mrr)                    # numpy array of bytes
con.register("mrr_arrow", pa.table({
    "parcel_id": tbl.column("parcel_id"),
    "wkb": pa.array(out_wkb),
}))
con.execute("""
    CREATE OR REPLACE TABLE parcel_mrr AS
    SELECT parcel_id, ST_GeomFromWKB(wkb) AS geom FROM mrr_arrow
""")

The change that matters is the dtype: shapely.from_wkb and shapely.to_wkb operate on NumPy arrays, so the GEOS work runs in a compiled ufunc loop while the Python interpreter is untouched. A comprehension over shapely.from_wkb(b) per element throws away exactly that advantage. Register the result as an Arrow table and re-ingest it with ST_GeomFromWKB so the custom output rejoins the columnar engine; the zero-copy mechanics of that registration are covered in Arrow interop patterns and the zero-copy Arrow to GeoPandas handoff.

Diagnostic Queries & Plan Validation #

Confirm the pushdown actually engaged the index rather than silently falling back to a nested loop — the failure mode that makes an SQL rewrite no faster than the loop it replaced. EXPLAIN ANALYZE is the arbiter.

EXPLAIN ANALYZE
SELECT p.point_id, r.region_id
FROM points p
JOIN regions r
  ON r.geom && p.geom
 AND ST_Contains(r.geom, p.geom);

Read the plan for three signals:

• Join node type. You want an index-driven scan on the && predicate, not a NESTED_LOOP_JOIN over the full region set — a nested loop is the signature of a defeated index, and it means the rewrite bought nothing.
• Candidate pruning. Compare actual rows entering ST_Contains against the raw pair count. If the exact-topology stage sees nearly $N \times M$ rows, the && stage is not pruning; re-check that the index exists and both layers share a CRS.
• No ROW_EXECUTION around the geometry. A row-at-a-time node wrapping an ST_ call means the kernel is not vectorizing — usually a scalar UDF or an implicit cast slipped into the projection.

Two one-liners confirm the runtime state:

-- Confirm the R-tree is registered before trusting any plan.
SELECT index_name, index_type FROM duckdb_indexes() WHERE table_name = 'regions';

-- 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;

For the Case 2 round-trip, the equivalent check is a timing sanity test: a shapely ufunc over a NumPy array should run in tens of milliseconds per hundred-thousand geometries. If it does not, you are almost certainly still handing it a Python list rather than the array from_wkb returns — the ufunc silently falls back to element-wise dispatch on non-array input.

Geometry Validation & Fallback Routing #

Both paths break on invalid geometry, but in different ways. In the SQL join, a self-intersecting polygon makes ST_Contains throw a GEOS error mid-scan; in the Shapely round-trip, an invalid ring parses fine through from_wkb and then produces a wrong or empty result from the ufunc. Validate at the boundary in both cases.

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

-- Repair in place, then rebuild the index so MBRs match the corrected geometry.
UPDATE regions SET geom = ST_MakeValid(geom) WHERE NOT ST_IsValid(geom);
DROP INDEX IF EXISTS idx_regions_geom;
CREATE INDEX idx_regions_geom ON regions USING RTREE (geom);

For the batched Shapely path, screen the array with the vectorized validity predicate and route only the offenders through repair, keeping the fast path fast:

valid = shapely.is_valid(geoms)          # vectorized boolean mask, no Python loop
if not valid.all():
    geoms = np.where(valid, geoms, shapely.make_valid(geoms))  # repair only the failures
result = shapely.minimum_rotated_rectangle(geoms)

When the target layer is too large to hold both inputs plus the R-tree inside memory_limit, do not iterate points in Python to save memory — that reintroduces the very loop this guide removes. Instead chunk the driving layer in SQL with a bounded LIMIT/OFFSET window or a partition key, running the same indexed predicate over each chunk so the R-tree stays resident and only the candidate set is bounded:

-- Chunked predicate join: the index stays hot, only the probe side is bounded.
SELECT p.point_id, r.region_id
FROM (SELECT * FROM points ORDER BY point_id LIMIT 1000000 OFFSET ?) p
JOIN regions r
  ON r.geom && p.geom
 AND ST_Contains(r.geom, p.geom);

This keeps memory flat while preserving the vectorized, index-pruned execution — the same chunking discipline the batch processing pipelines guide applies to writes. When the predicate output must land back in Python as a GeoDataFrame rather than stay in the engine, hand it off through the binary path in DuckDB to GeoPandas sync instead of rebuilding geometries row by row.

Related #

See also

• Optimizing point-in-polygon queries in DuckDB — the two-stage && plus exact-topology pattern the pushdown relies on.
• Arrow interop patterns — moving the WKB column across the boundary without an extra copy for the Shapely round-trip.
• Spatial joins and proximity filters — join-side tuning for predicates pushed into SQL.

Up: Shapely Integration · Python & DuckDB Integration Workflows