GeoPandas to DuckDB Spatial Migration

Porting a GeoPandas workflow to DuckDB is less a line-by-line rewrite than a shift in execution model: GeoDataFrame operations are eager, single-threaded, and cap out the moment the frame no longer fits in RAM, whereas DuckDB streams the same spatial logic through vectorized SQL kernels across every core and spills to disk instead of dying. This guide, part of Migrating to DuckDB Spatial, maps the everyday GeoPandas operators — sjoin, dissolve, buffer, to_crs, clip, overlay — to their DuckDB equivalents, shows how to move geometry both directions over Arrow without a serialization tax, and treats the hybrid split (DuckDB for the heavy joins and aggregation, GeoPandas for last-mile visualization) as the destination rather than a wholesale replacement.

The Execution Model Shift #

Before touching syntax, internalize why the two runtimes behave differently, because most failed migrations come from carrying GeoPandas assumptions into SQL. A GeoDataFrame is a pandas frame with a geometry column of Shapely objects. Every operation is materialized eagerly: gdf.buffer(500) builds a full new GeoSeries in memory before the next line runs, gdf.sjoin(...) constructs an in-memory spatial index and evaluates it on a single Python thread, and the whole frame — decoded geometry, index, and result — coexists in the heap at peak. The ceiling is hard: when the working set exceeds RAM, the process is OOM-killed, not slowed.

DuckDB inverts every one of those properties. Geometry lives as WKB-backed GEOMETRY in a columnar store; operators execute as vectorized kernels over 2,048-value chunks; the optimizer pushes predicates down so an indexed join never materializes the full cross product; and work is streamed through a pipeline that spills over-budget state to temp_directory rather than failing. The practical consequence is that a dissolve or join that GeoPandas can only run by first downsampling becomes a single streaming SQL statement in DuckDB. The trade-off runs the other way for small, iterative, interactive geometry editing — where Shapely’s object model and matplotlib rendering are simply more convenient — which is exactly why the endpoint is a hybrid, not a migration to zero GeoPandas.

Runtime Configuration & Memory Guardrails #

The bridge is only as reliable as the session you open it in. GeoPandas users are accustomed to an implicit “all of RAM” budget; DuckDB wants that budget stated so it can decide when to spill rather than when to fault. Load the extension and set the guardrails before registering any frame.

import duckdb

con = duckdb.connect(":memory:")
con.execute("INSTALL spatial; LOAD spatial;")  # per-session; ST_ functions need it loaded

# Match physical cores. GeoPandas ran on one thread; the win here is parallelism, but
# hyperthread siblings share SIMD units, so oversubscribing slows kernels, not speeds them.
con.execute("SET threads = 8;")

# State the ceiling the GeoDataFrame never had. Too low → spill mid-join; too high on a
# shared host → the OS OOM-kills the whole Python process, taking the notebook with it.
con.execute("SET memory_limit = '12GB';")

# Drop row-order guarantees so joins and GROUP BY run in parallel. Re-impose order with an
# explicit ORDER BY only on the final result you hand back to GeoPandas.
con.execute("SET preserve_insertion_order = false;")

# Spill target for over-budget joins/unions. Point at local NVMe, never a network mount —
# a spilling ST_Union on a slow disk turns into an I/O cliff that looks like a hang.
con.execute("SET temp_directory = '/var/lib/duckdb/spill';")

Before you move a production pipeline across, confirm these preconditions hold — each maps to a failure mode that is silent until it is catastrophic:

• memory_limit exceeds the uncompressed footprint of the largest single relation the query holds in flight, not the size of the input file on disk.
• The geometry column is decoded to native GEOMETRY at ingestion with ST_GeomFromWKB, not parsed from WKT text per query — the difference is analyzed in ST_Geometry vs WKB.
• The CRS of every input is recorded out-of-band, because the SRID does not survive the Arrow crossing (see edge cases below).
• A local spill directory is configured so an over-budget dissolve degrades gracefully instead of aborting.

Trade-off: a generous memory_limit keeps joins spill-free but widens the blast radius on a shared host, where one runaway query can starve co-tenants. Prefer a per-connection ceiling plus SET max_temp_directory_size over trusting one global number; the durability and spill mechanics of the underlying store are covered in in-memory vs disk storage.

Moving Data Across the Arrow Bridge #

The migration lives or dies on the handoff. Never round-trip geometry through GeoJSON text or a .to_file() intermediate; move it as WKB over Apache Arrow, which both runtimes speak natively. The bridge is bidirectional and the shape of the pipeline is a hybrid: heavy set logic runs in DuckDB, results cross back to GeoPandas only for the last-mile rendering that GeoPandas is genuinely better at.

GeoPandas into DuckDB #

The most portable path serializes the active geometry to WKB, registers the resulting pandas frame, and decodes it once into a native column. GeoDataFrame.to_wkb() returns an ordinary DataFrame in which each geometry column is replaced by WKB bytes, which DuckDB reads directly:

import geopandas as gpd

gdf = gpd.read_file("parcels.gpkg")          # active geometry column: "geometry"
wkb_df = gdf.to_wkb()                          # geometry column becomes WKB bytes
con.register("parcels_df", wkb_df)             # zero-copy view over the pandas frame

# Decode WKB → native GEOMETRY exactly once, at ingestion, so every later ST_ call is
# a kernel over decoded WKB rather than a per-row parse.
con.execute("""
    CREATE TABLE parcels AS
    SELECT * REPLACE (ST_GeomFromWKB(geometry) AS geometry)
    FROM parcels_df
""")

When you would rather avoid the pandas hop, GeoDataFrame.to_arrow(geometry_encoding="WKB") yields an Arrow stream whose geometry column is already WKB; wrap it with pyarrow.table(...) and register that instead. The deeper zero-copy mechanics — chunk alignment, when a copy is actually elided, and how to keep the buffer from being duplicated — are the subject of Arrow interop patterns and the zero-copy Arrow to GeoPandas handoff.

DuckDB back to GeoPandas #

Coming the other way, emit geometry as WKB with ST_AsWKB, fetch the Arrow table, and rebuild the GeoDataFrame from that binary column. The step that GeoPandas users forget is re-attaching the CRS by hand, because DuckDB’s GEOMETRY type carries no SRID for it to hand back:

arrow_tbl = con.execute("""
    SELECT region_id, total_area, ST_AsWKB(geometry) AS geometry
    FROM dissolved_regions
""").fetch_arrow_table()

df = arrow_tbl.to_pandas()
result = gpd.GeoDataFrame(
    df,
    geometry=gpd.GeoSeries.from_wkb(df["geometry"]),  # decode WKB → Shapely
    crs="EPSG:3857",                                   # RE-ATTACH: DuckDB did not carry it
)

The end-to-end conversion, including the streaming fetch_record_batch variant that never holds the whole result in memory, is worked through in the DuckDB-to-GeoPandas sync guide and its converting DuckDB queries to a GeoDataFrame efficiently walkthrough.

Mapping GeoPandas Operations to DuckDB SQL #

Most GeoPandas pipelines are a short chain of the same half-dozen verbs. Each has a direct, more scalable SQL form. The table is the fast reference; the notes below it flag where the semantics differ enough to break a naive translation.

A dissolve is the clearest illustration of the model shift, because GeoPandas must group and union entirely in memory while DuckDB streams the grouping:

-- gdf.dissolve(by="region", aggfunc="sum") →
SELECT
    region_id,
    ST_Union(geom)        AS geometry,      -- topological dissolve, GEOS-backed
    SUM(population)        AS total_pop
FROM parcels
GROUP BY region_id;

Reduce precision ahead of a union to strip IEEE-754 slivers and speed the merge — ST_Union(ST_ReducePrecision(geom, 0.001)) — the same discipline that keeps vectorized aggregations stable at scale.

clip and overlay are the two operators most often mistranslated, because both are geometry-producing rather than boolean, and both need the && pre-filter that GeoPandas hides. A clip against a single mask geometry is an intersection kept only where the inputs actually overlap:

-- gdf.clip(mask_polygon) → intersect and keep the portion inside the mask.
SELECT id, ST_Intersection(geom, (SELECT geom FROM mask)) AS geometry
FROM features, mask
WHERE geom && mask.geom                        -- prune before the costly intersection
  AND NOT ST_IsEmpty(ST_Intersection(geom, mask.geom));

An overlay(..., how="intersection") is the same idea against a whole second layer rather than one mask — a spatial join that emits the clipped intersection of every overlapping pair:

-- left.overlay(right, how="intersection") → join, then produce the shared geometry.
SELECT a.id AS a_id, b.id AS b_id,
       ST_Intersection(a.geom, b.geom) AS geometry
FROM layer_a a
JOIN layer_b b
  ON a.geom && b.geom                          -- Stage 1: MBR overlap (R-tree servable)
WHERE ST_Intersects(a.geom, b.geom);           -- Stage 2: exact test before intersecting
``` The `to_crs` translation is the one most likely to silently corrupt results: `ST_Transform` needs both the source and target authority codes spelled out because, unlike a `GeoDataFrame`, a DuckDB column has no attached CRS to infer from. The reprojection rules and axis-order gotchas are detailed in [CRS mapping and transformations](/duckdb-spatial-architecture-fundamentals/crs-mapping-transformations/).

The `sjoin` translation is deep enough to warrant its own treatment — index construction, the `how="inner"`/`"left"` variants, and validating that the parallel join returns exactly the GeoPandas result set are covered in [replacing GeoPandas sjoin with DuckDB](/migrating-to-duckdb-spatial/geopandas-to-duckdb-migration/replacing-geopandas-sjoin-with-duckdb/). Where the pipeline genuinely needs per-geometry Python logic that has no `ST_` equivalent, keep Shapely in the loop with the vectorized-callback pattern in [Shapely integration](/python-duckdb-integration-workflows/shapely-integration/) rather than falling back to a row-at-a-time UDF.

## Execution Plan Validation

The whole reason to migrate is parallel, index-assisted execution; `EXPLAIN ANALYZE` is how you prove you got it rather than a serial nested loop that merely happens to run inside DuckDB. Validate the translated join before trusting it:

```sql
EXPLAIN ANALYZE
SELECT l.id, r.region_id
FROM points l
JOIN regions r
  ON r.geom && l.geom                       -- Stage 1: MBR overlap (R-tree servable)
WHERE ST_Within(l.geom, r.geom);            -- Stage 2: exact topology on survivors

Read the measured plan and assert three things:

• Join operator. You want an index/hash-style scan driven by &&, not a NESTED_LOOP_JOIN over the full right table. A nested loop is the exact behaviour you left GeoPandas to escape — it means the R-tree was never built or the indexed column was wrapped in a function and hidden from the planner.
• No ROW_EXECUTION fallback. A ROW_EXECUTION node or a scalar Python UDF inside the join phase bypasses SIMD kernels and reintroduces the single-threaded ceiling. Replace the UDF with the native ST_ equivalent.
• Cardinality drift. Compare actual vs estimated rows at the join. Drift beyond ~30% signals skewed geometry density or stale statistics and predicts plans that flip between builds. The MBR pre-filter should prune the candidate set to a small fraction before topology runs; the two-stage mechanics are worked out in point-in-polygon optimization.

Capture the machine-readable plan with EXPLAIN (FORMAT JSON) … so the assertion can run in CI rather than by eye — the harness at the end of this guide does exactly that.

Performance Trade-offs #

Migrating is a set of bounded trade-offs, not a free win in every direction. Quantify against your own data, but the typical ranges and the decision that governs each are:

The honest boundary: GeoPandas remains the right tool when the frame is small enough to sit in memory comfortably and the work is interactive — plotting, tweaking a single geometry, or feeding a scikit-learn model that already wants a DataFrame. The migration target is not “delete GeoPandas” but “stop asking GeoPandas to do the multi-gigabyte join.” Everything above roughly a million features, or any pipeline that OOM-kills today, belongs on the DuckDB side of the Arrow bridge; the crossover point where DuckDB overtakes an in-memory frame is measured directly in the performance crossover benchmarks.

Edge Cases & Anti-Patterns #

CRS metadata does not cross the Arrow boundary (silent wrong answers). A GeoDataFrame stores its CRS in gdf.crs; DuckDB’s GEOMETRY type stores dimensionality but no SRID. Neither to_wkb() nor to_arrow() propagates the CRS into DuckDB, and ST_AsWKB carries none back. If you skip the manual crs= on the return trip, GeoPandas will happily plot metre-space coordinates on a degree-space basemap. Record the CRS out-of-band and re-attach it on every crossing.

# Anti-pattern: CRS is silently lost, so the returned frame is unprojected.
gdf_back = gpd.GeoDataFrame(df, geometry=gpd.GeoSeries.from_wkb(df["geometry"]))

# Fix: re-attach the CRS you tracked yourself.
gdf_back = gpd.GeoDataFrame(
    df, geometry=gpd.GeoSeries.from_wkb(df["geometry"]), crs="EPSG:3857"
)

The recurring drift this causes and how to detect it are catalogued in fixing CRS drift in GeoDataFrame conversion.

Geometry column naming mismatch. GeoPandas defaults its geometry column to geometry, but a frame read from a shapefile or renamed by a prior step may use geom, the_geom, or wkb_geometry. to_wkb() serializes whatever column the frame considers geometry; if your CREATE TABLE hard-codes geometry, the decode silently reads the wrong column or errors. Inspect gdf.geometry.name before templating the SQL, and alias explicitly.

Active-geometry ambiguity with multiple geometry columns. A GeoDataFrame can hold several geometry columns but has exactly one active one, set by set_geometry. to_wkb() converts all geometry-typed columns to WKB, not just the active one, so the DuckDB table gets several WKB columns and you must decode each with its own ST_GeomFromWKB. Conversely, when building the frame on the way back, GeoPandas has no notion of “active” until you tell it — pass geometry= or call set_geometry("geometry"), or .plot() and .crs will raise.

Empty and null geometries. gpd.GeoSeries.from_wkb maps SQL NULL to a None geometry, but an empty geometry (GEOMETRYCOLLECTION EMPTY) survives as a valid-but-empty Shapely object that many downstream ops treat inconsistently. Filter with WHERE NOT ST_IsEmpty(geom) on the DuckDB side before the handoff if your GeoPandas code assumes non-empty features.

Row-at-a-time UDF as a sjoin crutch. Registering a Python function that calls Shapely per row to reproduce a GeoPandas predicate reintroduces exactly the single-threaded ceiling you migrated to escape. Use the native ST_ predicate, or the batched-callback pattern in vectorizing Shapely predicates over DuckDB.

Query Regression Analysis #

A migration is not done when the SQL runs — it is done when the SQL provably returns the same features the GeoDataFrame did, on every build. Wrap the two implementations in a regression harness that compares geometry set-equality (after normalization, since operation order and precision can differ) and asserts the DuckDB plan never regresses to a serial nested loop. This slots into the orchestration covered in the Python and DuckDB integration workflows guide.

import duckdb
import json
import shapely
import geopandas as gpd

con = duckdb.connect(":memory:")
con.execute("INSTALL spatial; LOAD spatial;")
con.execute("SET threads = 8; SET memory_limit = '12GB';")

# --- Reference result: the legacy GeoPandas path -------------------------------
left = gpd.read_file("points.gpkg")
right = gpd.read_file("regions.gpkg")
gpd_out = left.sjoin(right, predicate="within")[["id", "region_id"]]

# --- Candidate result: the migrated DuckDB path --------------------------------
con.register("pts", left.to_wkb())
con.register("rgn", right.to_wkb())
con.execute("""
    CREATE TABLE points  AS SELECT * REPLACE (ST_GeomFromWKB(geometry) AS geometry) FROM pts;
    CREATE TABLE regions AS SELECT * REPLACE (ST_GeomFromWKB(geometry) AS geometry) FROM rgn;
""")
JOIN_SQL = """
SELECT l.id, r.region_id
FROM points l JOIN regions r
  ON r.geometry && l.geometry
WHERE ST_Within(l.geometry, r.geometry)
"""
duck_out = con.execute(JOIN_SQL).df()

# --- 1. Assert identical membership -------------------------------------------
gpd_pairs  = set(map(tuple, gpd_out.to_numpy()))
duck_pairs = set(map(tuple, duck_out.to_numpy()))
assert gpd_pairs == duck_pairs, (
    f"Row-set mismatch: only-GeoPandas={gpd_pairs - duck_pairs}, "
    f"only-DuckDB={duck_pairs - gpd_pairs}"
)

# --- 2. Assert the plan stayed parallel and index-assisted --------------------
plan = json.loads(con.execute(f"EXPLAIN (FORMAT JSON) {JOIN_SQL}").fetchone()[1])
BANNED = {"NESTED_LOOP_JOIN", "ROW_EXECUTION"}

def walk(node):
    yield node.get("name", "")
    for child in node.get("children", []):
        yield from walk(child)

nodes = set(walk(plan[0] if isinstance(plan, list) else plan))
assert not (nodes & BANNED), f"Plan regressed to serial execution: {nodes & BANNED}"

For a geometry-producing operation such as a dissolve or overlay, membership equality is not enough — compare the geometries themselves after normalizing vertex order and snapping precision, since GEOS in GeoPandas and GEOS in DuckDB can differ in the last floating-point digit:

# Normalize before comparing so equivalent geometries that differ only in vertex
# ordering or ULP-level precision are treated as equal.
def normalized_wkb(geom):
    return shapely.to_wkb(shapely.normalize(shapely.set_precision(geom, 1e-6)))

gpd_geoms  = {r: normalized_wkb(g) for r, g in zip(gpd_out.region_id, gpd_out.geometry)}
duck_geoms = {r: normalized_wkb(g) for r, g in zip(duck_out.region_id,
                                                    gpd.GeoSeries.from_wkb(duck_out.geometry))}
assert gpd_geoms == duck_geoms, "Dissolved geometry diverged after normalization"

Three signals are worth tracking across builds: the symmetric difference of the two row sets (any non-empty set is a hard regression), the appearance of NESTED_LOOP_JOIN or ROW_EXECUTION in the plan (a lost index or a stray UDF), and per-operator operator_timing deltas that localize a slowdown to a specific node. Run this harness in CI on a fixed fixture so a divergence surfaces in a pull request rather than in a production map.

Related #

See also

• Replacing GeoPandas sjoin with DuckDB — the full index-assisted rewrite of the single most common GeoPandas operation, with how and predicate variants.
• DuckDB-to-GeoPandas sync and converting DuckDB queries to a GeoDataFrame efficiently — the return-trip conversion in depth, including streaming batches.
• Arrow interop patterns — the zero-copy mechanics of the bridge both directions rely on.
• Shapely integration — keeping per-geometry Python logic in the pipeline without a row-at-a-time UDF.

Up: Migrating to DuckDB Spatial