Replacing GeoPandas sjoin with DuckDB

Swapping geopandas.sjoin() for a DuckDB spatial join removes the single hardest ceiling in a GeoPandas pipeline — a spatial index built entirely in one process’s heap and walked on one thread — and this walkthrough, part of the GeoPandas to DuckDB Spatial migration guide, replaces it with a two-stage RTREE-and-&& join that prunes candidates in parallel across every core and streams instead of materializing the whole result.

Root-Cause Analysis: Why sjoin Hits a Wall #

geopandas.sjoin(left, right, predicate="within") is convenient because it hides three expensive steps behind one call, and each is a distinct scaling limit:

• The index is built in memory, every time. sjoin constructs an STRtree (via the Shapely 2.0 / GEOS bindings) over the right-hand geometries on each call. Nothing persists between calls, so an interactive loop rebuilds the same tree repeatedly, and the tree plus both decoded geometry sets must all coexist in the heap. When that sum exceeds RAM, the process is OOM-killed outright rather than degrading.
• Evaluation is single-threaded. The query-and-refine loop over the STRtree runs on one Python thread. A 24-core box gives sjoin no advantage; wall-clock time scales with the candidate count on a single core.
• The refinement predicate runs per candidate in Python. After the bounding-box query, the exact within/intersects test is a GEOS call per surviving pair, dispatched through Python object overhead. The pruning is good but the survivors are still checked one Python call at a time.

DuckDB attacks all three. The R-tree is a persistent on-disk index you build once; the && MBR pre-filter and the exact predicate both execute as vectorized kernels the optimizer parallelizes across threads; and the join streams, spilling to temp_directory under budget pressure rather than requiring the whole result to fit in memory. The layout that makes the geometry decode cheap enough to vectorize is analyzed in ST_Geometry vs WKB.

The complexity story is the same on paper and different in practice. Both engines aim to avoid the naive $O(N \times M)$ all-pairs comparison by using a bounding-box index to reach roughly $O(N \log M)$ candidate lookups. GeoPandas achieves that bound but pays it on one core with Python-object dispatch on every refinement; DuckDB achieves the same asymptotic bound and then divides the constant factor by the core count while keeping the refinement in a compiled kernel. On a join whose candidate set is large, that constant-factor difference is the whole migration: identical big-O, an order of magnitude in wall-clock, and — more importantly — a hard memory ceiling replaced by a soft spill boundary.

Deterministic Configuration #

Load the extension and pin the session so the rewrite is reproducible. Only the knobs the join actually depends on are set here.

import duckdb

con = duckdb.connect("migration.duckdb")   # on-disk DB so the R-tree persists between runs
con.execute("INSTALL spatial; LOAD spatial;")

# Match physical cores. This is the whole point of the migration: sjoin used one thread,
# the DuckDB join uses all of them. Oversubscribing past physical cores just contends.
con.execute("SET threads = 8;")

# The join holds both inputs, the R-tree, and the result region at once. Cap below total
# RAM so an unpruned spike spills to disk instead of triggering an OS OOM-kill.
con.execute("SET memory_limit = '8GB';")

# Unlock parallel index construction by dropping order preservation; the join output was
# never ordered anyway, so re-impose ORDER BY only if a caller needs it.
con.execute("SET preserve_insertion_order = false;")

Confirm the preconditions before rewriting query logic:

• spatial extension loaded in the session
• Both layers share one CRS — DuckDB does not read the GeoDataFrame.crs, so align before ingesting, per CRS mapping and transformations
• The right-hand (indexed) geometry is a native GEOMETRY column, decoded once with ST_GeomFromWKB
• An R-tree index exists on that column before the join runs

Optimized Execution Pattern #

Start from the GeoPandas baseline: assign each of several million points to the administrative polygon that contains it.

# Before — GeoPandas: in-memory STRtree, single-threaded refinement.
import geopandas as gpd

points  = gpd.read_file("gps_points.gpkg")     # ~5M points
regions = gpd.read_file("admin_regions.gpkg")  # ~10K polygons
joined = gpd.sjoin(points, regions, how="inner", predicate="within")

The DuckDB replacement is a build-index-once, then two-stage-join pattern. Ingest both frames over WKB, index the polygon layer, and split the predicate into a cheap MBR overlap the R-tree serves and an exact ST_Within that runs only on survivors:

con.register("pts_df", points.to_wkb())
con.register("rgn_df", regions.to_wkb())
con.execute("""
    CREATE TABLE points  AS SELECT * REPLACE (ST_GeomFromWKB(geometry) AS geometry) FROM pts_df;
    CREATE TABLE regions AS SELECT * REPLACE (ST_GeomFromWKB(geometry) AS geometry) FROM rgn_df;

    -- Build the R-tree ONCE on the polygon layer; it persists in the on-disk database.
    CREATE INDEX idx_regions_geom ON regions USING RTREE (geometry);
""")

-- After — DuckDB: index-assisted, parallel, streaming spatial join.
SELECT p.id, p.metric, r.region_id, r.name
FROM points p
JOIN regions r
  ON r.geometry && p.geometry            -- Stage 1: MBR overlap, served by the R-tree
WHERE ST_Within(p.geometry, r.geometry); -- Stage 2: exact topology, survivors only

The load-bearing change is entirely in predicate placement. The && operator against the bare indexed column in the ON clause is index-servable, so the planner replaces what would be a full scan with an R-tree scan that hands ST_Within a candidate set one to two orders of magnitude smaller. GeoPandas did the same bounding-box pruning internally, but serially and in one heap; here the pruning and the exact check both parallelize and stream.

Matching how and predicate variants #

sjoin exposes how and predicate; both map cleanly to SQL once you see the correspondence.

The how="left" case has one trap: to keep unmatched left rows, the exact predicate must move into the ON clause, because a WHERE ST_Within(...) filters the null-padded rows back out and silently turns your left join into an inner one.

-- how="left": keep every point, null region columns where nothing contains it.
SELECT p.id, r.region_id
FROM points p
LEFT JOIN regions r
  ON r.geometry && p.geometry
 AND ST_Within(p.geometry, r.geometry);   -- exact test in ON, NOT in WHERE

For predicate="dwithin", keep the && overlap in front so the R-tree still drives the scan, then apply ST_DWithin; the proximity variants and their index behaviour are covered in spatial joins and proximity filters.

Diagnostic Queries & Plan Validation #

Prove the index is doing the work rather than guessing from wall-clock time. Confirm the index is registered, then read the measured plan:

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

EXPLAIN ANALYZE
SELECT p.id, r.region_id
FROM points p JOIN regions r
  ON r.geometry && p.geometry
WHERE ST_Within(p.geometry, r.geometry);

Check three things in the output:

• Join type. An index/hash-style scan driven by && is correct. A NESTED_LOOP_JOIN over the full regions table means the index was not used — usually because it was never built, or the indexed column was wrapped in a function (ST_Buffer(r.geometry, …) && …) and hidden from the planner. Keep r.geometry bare on the left of &&.
• Estimated vs actual rows at the join. Large drift signals stale statistics or a non-selective MBR stage; re-check that both layers share a CRS and that the R-tree exists. A well-pruned join hands only a small fraction of pairs to ST_Within.
• Operator timing. If the exact-topology node dominates total time, the MBR stage is not pruning — the classic cause is degree-space envelopes against a projected layer. The same two-stage timing analysis, with thresholds, is worked through in optimizing point-in-polygon queries.

Watch memory and spill live while a large join runs:

SELECT tag, memory_usage_bytes, temporary_storage_bytes
FROM duckdb_memory() ORDER BY memory_usage_bytes DESC;

Geometry Validation & Result Equality #

Two things must hold before you retire the GeoPandas code: the polygon layer must be topologically valid, and the DuckDB result set must equal the sjoin result exactly.

Invalid polygons make both the MBR and the exact test unreliable, so quarantine and repair before indexing:

-- Isolate offenders, repair deterministically, then rebuild so MBRs match fixed geometry.
UPDATE regions SET geometry = ST_MakeValid(geometry) WHERE NOT ST_IsValid(geometry);

DROP INDEX IF EXISTS idx_regions_geom;
CREATE INDEX idx_regions_geom ON regions USING RTREE (geometry);

Then assert set-equality against the GeoPandas baseline — the migration is only trustworthy once the two produce the identical (point, region) membership:

import geopandas as gpd

gpd_pairs  = set(map(tuple, gpd.sjoin(points, regions, how="inner",
                                      predicate="within")[["id", "region_id"]].to_numpy()))
duck_pairs = set(map(tuple, con.execute("""
    SELECT p.id, r.region_id
    FROM points p JOIN regions r
      ON r.geometry && p.geometry
    WHERE ST_Within(p.geometry, r.geometry)
""").df().to_numpy()))

assert gpd_pairs == duck_pairs, (
    f"sjoin migration diverged: only-GeoPandas={len(gpd_pairs - duck_pairs)}, "
    f"only-DuckDB={len(duck_pairs - gpd_pairs)} rows"
)

A non-empty symmetric difference almost always traces to one of three causes: a CRS mismatch (degree-space points tested against a projected layer), points sitting exactly on a boundary where GEOS within is edge-sensitive, or empty geometries that the two engines handle differently — filter them with WHERE NOT ST_IsEmpty(geometry) on both sides before comparing. When memory is still breached on a very large join despite pruning, chunk the point side and union the results, or push the reduction upstream with the grid pre-aggregation shown in point-in-polygon optimization, and hand only the reduced result back to GeoPandas via DuckDB-to-GeoPandas sync.

Related #

See also

• Spatial joins and proximity filters — the full family of MBR-pruned join predicates the predicate variants map onto, including point-in-polygon optimization.
• DuckDB-to-GeoPandas sync — moving the joined result back into a GeoDataFrame for rendering.

Up: GeoPandas to DuckDB Spatial migration · Migrating to DuckDB Spatial