DuckDB Spatial Analytical SQL for GIS

DuckDB vs GeoPandas for In-Memory Analysis

GeoPandas does not get slower as the data grows — it stops, and the transition between those two states is what the comparison is actually about. This walkthrough, part of the engine comparisons reference, works through a six-step analysis in both engines, sets out the three arrangements for using them together, and names the four habits that do not survive the move.

Root-Cause Analysis: where the ceiling comes from

The memory wall is not one limit but four things multiplying, which is why it arrives sooner than the file size suggests.

  • Shapely objects are large. Every geometry becomes a Python object with its own header, and a frame of a million polygons holds a million of them. The in-memory footprint is routinely several times the file on disk.
  • Operations copy. Most frame operations return a new frame rather than mutating in place, so a chain of six steps can hold two or three full copies at its peak, not one.
  • There is no spill. When the working set exceeds available memory the process is killed. There is no degraded mode, which is what makes the failure a cliff rather than a slope.
  • Everything is single-threaded. The operations that would benefit most from parallelism — joins, overlays, dissolves — are the ones that get none, so the time to the ceiling is longer than it needs to be as well.
  • The index is rebuilt per operation. A spatial join builds an in-memory index over the right-hand frame each time it runs, so a repeated join pays that cost repeatedly.

The distinguishing question is whether the analysis ends small. If a large input reduces to a small result, the hybrid arrangement removes the ceiling entirely without a rewrite.

Six steps of one analysis, in each engine Read, filter, join, dissolve, compute and write — comparable at small size, categorically different once the frame exceeds memory. STEP GEOPANDAS DUCKDB read the file holds the whole frame streams filter to a region after materialising everything pushed into the scan spatial join single-threaded, in-memory index parallel, R-tree pruned dissolve every intermediate resident groupable, bounded intermediates compute an area column comparable comparable write the result comparable comparable

Two rows are equal, two differ by a factor, and two decide whether the job finishes.

Deterministic Configuration 

import duckdb, geopandas as gpd
from shapely import from_wkb

con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial")

# The engine now carries the memory the frame used to. Give it a spill target
# so the operation that used to kill the process degrades instead.
con.execute("SET memory_limit = '8GB'")
con.execute("SET threads = 8")
con.execute("SET temp_directory = '/var/tmp/duckdb_gpd'")

Optimized Execution Pattern 

# ANTI-PATTERN: everything in the frame. Each step holds a full copy, the join
# is single-threaded, and the process dies at whatever step first exceeds RAM.
parcels = gpd.read_parquet("parcels.parquet")            # whole layer resident
zones   = gpd.read_parquet("zones.parquet")
joined  = parcels.sjoin(zones, predicate="within")       # single-threaded
result  = joined.dissolve(by="zone_id", aggfunc="sum")   # every intermediate resident
# PATTERN: the reduction happens in the engine, and the frame is built from
# what is left — which is small enough that none of the above matters.
result = con.execute("""
    SELECT z.zone_id,
           ST_Union_Agg(p.geom) AS geom,
           sum(p.value)         AS total_value
    FROM read_parquet('parcels.parquet') p
    JOIN read_parquet('zones.parquet')   z ON z.geom && p.geom
    WHERE ST_Within(p.geom, z.geom)
    GROUP BY z.zone_id
""").fetch_arrow_table()

gdf = gpd.GeoDataFrame(
    result.to_pandas().drop(columns=["geom"]),
    geometry=from_wkb(result.column("geom").to_numpy(zero_copy_only=False)),
    crs="EPSG:27700",          # DuckDB does not carry this; supply it explicitly
)
Three arrangements, three purposes All-GeoPandas for exploration below the ceiling, all-DuckDB for scheduled pipelines, and the hybrid for analyses that start large and end small. ARRANGEMENT RIGHT WHEN WHY all GeoPandas exploratory, comfortably in memory iteration speed is unmatched all DuckDB a scheduled pipeline to a file no interactive step to make pleasant DuckDB reduces, GeoPandas receives the analysis starts large, ends small the common case, and the usual end state

The third row is where most ported workflows land, and it is not a compromise.

Choosing where the boundary goes

The boundary belongs at the last point where the row count drops sharply, and finding it is usually a matter of counting rather than judgement. Run the pipeline’s filters and joins as counts alone, without materialising anything, and look for the step after which the result is small enough to be comfortable in memory. That step is the handover.

Two mistakes recur. Putting the boundary too early means the frame still receives millions of rows and the ceiling is unchanged — a query that filters nothing has moved the work without reducing it. Putting it too late means an operation GeoPandas does well, such as producing a plot or calling a library that expects Shapely, has been reimplemented awkwardly in SQL for no gain. The right position is almost always immediately after the last aggregation.

-- Find the boundary by counting rather than by materialising. The step after
-- which this drops sharply is where the handover belongs.
SELECT 'raw'      AS stage, count(*) FROM read_parquet('parcels.parquet')
UNION ALL SELECT 'filtered', count(*) FROM read_parquet('parcels.parquet') WHERE land_use = 'residential'
UNION ALL SELECT 'joined',   count(*) FROM read_parquet('parcels.parquet') p
                              JOIN read_parquet('zones.parquet') z ON z.geom && p.geom
UNION ALL SELECT 'grouped',  count(DISTINCT zone_id) FROM read_parquet('zones.parquet');
Four habits that do not survive the move Chained frames become one planned query; the index becomes an explicit key; dissolve defaults become named aggregates; printing becomes selecting from a CTE. GEOPANDAS HABIT WHAT REPLACES IT WHY IT IS BETTER chaining frame operations one query, planned as a whole the engine reorders and fuses steps the frame index as identity an explicit key column SQL has no positional index dissolve’s default aggregation naming an aggregate per column the default was rarely what you meant printing an intermediate selecting from a named CTE the whole dataset, not the first five

Three of the four are improvements once the initial friction passes.

Diagnostic Queries & Plan Validation 

# The comparison that matters is peak resident memory, not runtime — the
# frame version fails on memory long before it becomes slow.
import resource
peak_kb = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
print(f"peak RSS: {peak_kb / 1024:.0f} MB")

Diagnostic — did the reduction actually reduce? Compare the row count crossing the boundary against the input. A handover that passes millions of rows has moved the work without removing the ceiling, which is the most common way a hybrid rewrite disappoints.

Geometry Validation & Fallback Routing

The rebuilt frame needs three things supplied that the query result does not carry, and all three fail silently: the active geometry column, the CRS, and a deliberate index. Asserting them at construction costs three lines and catches the class of defect that otherwise appears as an area computed in the wrong units.

assert gdf.crs is not None, "CRS was not supplied — every later measurement is wrong"
assert gdf.geometry.name == "geometry", "the active geometry column is not the one intended"
assert len(gdf) == result.num_rows, "rows were lost rebuilding the frame"

What the hybrid costs in practice

The arrangement is not free, and the costs are worth naming so the decision is made with them rather than despite them. There are now two languages in one analysis, so a reader has to follow the boundary between them. The reduced result has to be rebuilt as a frame, which means supplying the CRS and the geometry column explicitly every time. And a bug can now live on either side of the boundary, which makes the first debugging session slower than it would have been in one language.

Against that: the ceiling disappears, the reduction parallelises, and the interactive part stays interactive. On any dataset large enough for the question to arise, that trade is decisively worth making — and the costs above shrink with familiarity while the benefits do not.

# A small helper removes most of the recurring cost, because the three things
# that must be supplied are supplied in exactly one place.
def to_gdf(arrow_tbl, crs: str, geom_col: str = "geom"):
    """Rebuild a GeoDataFrame from an Arrow result, with the CRS supplied
    explicitly because nothing in the result carries it."""
    df = arrow_tbl.to_pandas().drop(columns=[geom_col])
    geoms = from_wkb(arrow_tbl.column(geom_col).to_numpy(zero_copy_only=False))
    return gpd.GeoDataFrame(df, geometry=geoms, crs=crs)

Frequently Asked Questions

At what size does GeoPandas stop being viable?

There is no single number, because the footprint depends on geometry complexity rather than row count — a million points is comfortable where a million dense polygons is not. The practical signal is when a routine operation starts taking minutes or the kernel starts being killed; both mean the ceiling is close, and neither gives much warning.

Do I have to rewrite everything?

No, and most successful moves do not. The hybrid — DuckDB reducing, GeoPandas receiving the reduced result — removes the ceiling while keeping the interactive ergonomics, and it usually touches only the first few lines of an analysis. A full rewrite is worth it for a scheduled pipeline and rarely worth it for exploratory work.

Why did my CRS disappear?

It never crossed. DuckDB geometry carries no reference frame, so nothing was lost in transit — it has to be supplied when the frame is rebuilt, from whatever your pipeline records. A frame built with crs=None computes areas and reprojections in whatever units the coordinates happen to be, silently.

What happened to my frame index?

SQL has no positional index, so it does not survive a round trip. Anything that relied on the index as identity — a join, a lookup, an alignment between two frames — needs an explicit key column instead. This is the habit change that causes the most surprise and the least difficulty once it is made.

Is the SQL version always faster?

For anything that reduces the row count, substantially. For an operation on an already-small frame, no, and the query is longer to write — which is a perfectly good reason to leave it in pandas. The gains come from parallelism, columnar scanning and index pruning, none of which apply to a thousand rows.

Can I keep plotting with GeoPandas?

Yes, and that is one of the main reasons the hybrid arrangement is the usual end state. Plotting, explore(), and the whole surrounding ecosystem expect a GeoDataFrame, and the reduced result is exactly the size at which those tools are pleasant to use.

Up: DuckDB Spatial vs PostGIS, GeoPandas and Sedona