How to Run ST_Intersects in DuckDB on GeoParquet

This guide is a complete, runnable recipe for loading DuckDB’s spatial extension, reading a GeoParquet file straight from disk or S3, and executing an ST_Intersects spatial join accelerated by an R-tree index, with a verification step that proves the index was used.

Context and prerequisites

ST_Intersects returns true when two geometries share at least one point, and it is the workhorse predicate for point-in-polygon and polygon-overlap joins. Run this recipe on DuckDB 1.0 or later (validated on 1.1.x) with the spatial extension; for s3:// inputs you also need httpfs. This page sits under DuckDB geospatial analytics on lakehouse tables, which covers the broader engine-selection picture — here we focus on getting one join correct and fast. The GeoParquet inputs are assumed to follow the standard WKB column encoding described in GeoParquet encoding standards.

Complete working solution

import duckdb

con = duckdb.connect()

# 1. Load spatial (geometry types + ST_* functions) and httpfs (for s3:// paths)
con.execute("INSTALL spatial; LOAD spatial;")
con.execute("INSTALL httpfs; LOAD httpfs;")

# If reading from S3, create a secret (skip for local files)
con.execute("""
CREATE OR REPLACE SECRET s3_src (
    TYPE S3, PROVIDER credential_chain, REGION 'us-east-1'
);
""")

# 2. Load points GeoParquet, decoding the WKB geometry column into GEOMETRY
con.execute("""
CREATE TABLE sensors AS
SELECT
    sensor_id,
    ST_GeomFromWKB(geometry) AS geom
FROM read_parquet('s3://lakehouse/sensors/*.parquet');
""")

# 3. Load polygons GeoParquet
con.execute("""
CREATE TABLE districts AS
SELECT
    district_id,
    name,
    ST_GeomFromWKB(geometry) AS geom
FROM read_parquet('s3://lakehouse/districts/*.parquet');
""")

# 4. Build an R-tree on the larger (points) table so the join can prune by bbox
con.execute("CREATE INDEX sensors_rtree ON sensors USING RTREE (geom);")

# 5. The ST_Intersects spatial join: assign each sensor to its district
result = con.execute("""
SELECT
    d.name           AS district,
    count(*)         AS sensor_count
FROM sensors s
JOIN districts d
  ON ST_Intersects(s.geom, d.geom)
GROUP BY d.name
ORDER BY sensor_count DESC;
""").fetchall()

for district, n in result:
    print(f"{district:20s} {n}")

con.close()

For a single-file local run, replace the s3:// globs with a filesystem path such as '/data/sensors.parquet' and drop the secret and httpfs lines — everything else is identical.

Step-by-step walkthrough

1. Load extensions. LOAD spatial registers the GEOMETRY type and every ST_* function, including ST_Intersects and ST_GeomFromWKB. LOAD httpfs is only needed for object storage; local files work with spatial alone. Extensions are cached after the first INSTALL, so repeated runs skip the download.

2. Decode WKB into GEOMETRY. GeoParquet stores geometry as Well-Known Binary in a Parquet BLOB column. ST_GeomFromWKB(geometry) reinterprets those bytes as DuckDB’s native geometry type. Skipping this step leaves you with an opaque blob and ST_Intersects will raise a binder error.

3. Materialize the polygon side. The districts table is small, so loading it into a base table (rather than a view) lets the optimizer see its cardinality and treat it as the build side of the join.

4. Create the R-tree. USING RTREE (geom) builds a bounding-box index. During the join, DuckDB first checks R-tree bounding boxes for overlap — a cheap integer comparison — and only calls the expensive GEOS ST_Intersects on candidate pairs that pass. Indexing the larger table gives the biggest win because most of its rows are eliminated at the bbox stage.

5. Run the join. ON ST_Intersects(s.geom, d.geom) is the join condition. Because the predicate is a plain two-argument ST_Intersects over the indexed column, the optimizer can extract a bounding-box comparison and drive it through the R-tree. Wrapping either geometry in a transform here would defeat the index.

Common errors and fixes

Verification

Confirm the R-tree is actually driving the join, not a full scan, with EXPLAIN ANALYZE.

plan = con.execute("""
EXPLAIN ANALYZE
SELECT count(*)
FROM sensors s JOIN districts d
  ON ST_Intersects(s.geom, d.geom);
""").fetchall()
for row in plan:
    print(row[1])

Look for an RTREE_INDEX_SCAN operator on the sensors side and a low Rows Scanned count relative to the table size. As an independent correctness check, verify that every sensor total is preserved across the join grouping:

total_joined = con.execute("""
SELECT sum(sensor_count) FROM (
    SELECT count(*) AS sensor_count
    FROM sensors s JOIN districts d ON ST_Intersects(s.geom, d.geom)
    GROUP BY d.name
);
""").fetchone()[0]
print("sensor-district pairs:", total_joined)

The two-stage filter-then-refine pattern shown above is why the R-tree matters so much: the exact GEOS predicate only ever runs on the handful of candidate pairs whose bounding boxes already overlap. To push this further — skipping whole row groups before decode by filtering on a numeric bbox covering column — see predicate pushdown optimization, and to run the same style of join against an Iceberg table instead of loose GeoParquet, see querying Iceberg tables with the DuckDB spatial extension. The canonical function reference is the DuckDB spatial functions documentation.

Making the Join Faster Still

The R-tree gets the join working; three further steps get it fast, and each is independent of the others.

The bounding-box pre-filter is worth applying even with the R-tree in place, because the two operate at different stages: the numeric predicate eliminates row groups during the Parquet read, before any geometry exists in memory, while the R-tree eliminates candidate pairs after the rows have been materialised. Doing both means the index is only ever asked about rows that survived the file-level prune.

Verifying Correctness, Not Just Speed

The right-hand case is the one that produces quiet undercounts. Rows with null geometry never satisfy any spatial predicate, so an inner join drops them without comment and a total that should have been complete is short by however many there were. Counting them before the join, and deciding explicitly whether they belong in the output, turns a silent loss into a stated one.

Scaling the Same Recipe

The recipe above holds up well past the scale most people expect, and the two adjustments that extend it further are both about avoiding materialisation.

The first is to skip the intermediate tables when the source is already well laid out. CREATE TABLE ... AS SELECT reads and stores everything before the join begins; querying read_parquet directly in the join lets DuckDB push the bounding-box predicate into the Parquet reader and never materialise the rows it will discard. The intermediate table is worth creating only for the small side, where the optimiser benefits from knowing the cardinality.

The second is to stream the output rather than collecting it. fetchall builds the entire result in Python memory, which is fine for an aggregate and wasteful for a join that returns millions of rows. Writing the result straight to Parquet with COPY ... TO keeps peak memory flat and produces an output that the rest of the pipeline can read without a further conversion.

Beyond that, the boundary is the one described in the section overview: when the join’s working set exceeds the node, the recipe stops applying and the workload belongs on a cluster. Everything up to that point — and it is a lot of data — runs comfortably in a single process with no infrastructure at all, which is the reason this pattern is worth having in the toolkit even on platforms whose default answer is distributed compute.

For reading a governed table rather than loose files, and for the predicate split that keeps the table planner involved, see querying Iceberg tables with the DuckDB spatial extension. That guide picks up exactly where this one leaves off, with the same two-stage predicate structure applied across the table boundary. The layout guidance that makes the bounding-box columns available in the first place is in spatial partitioning and indexing strategies. Without those columns, the R-tree is doing all the work alone and the file-level prune never happens.

The two mechanisms compose, and a table that supports both will answer this query in a fraction of the time either alone would achieve. Check both in the plan before concluding the query is as fast as it can be. A plan that shows one but not the other has headroom left in it. Both mechanisms are cheap to add and neither is enabled by default.