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.