This guide is part of Orchestrating Spatial ETL Pipelines, the reference for turning extract, clean, and load scripts into scheduled, restartable spatial data products.

PostGIS vs DuckDB Spatial for Analytical Loads

The Load stage is the one decision in a spatial pipeline that is expensive to reverse. Change your extraction library and you rewrite one task; change your load target after a year of production and you migrate every downstream query, dashboard, and materialized view that depends on it. For analytical spatial workloads — the ones that scan millions of features to compute coverage statistics, spatial joins, or tiled aggregates rather than serving single-feature lookups to a web map — the two engines that dominate the shortlist are PostGIS and the DuckDB spatial extension. They look interchangeable in a tutorial and behave nothing alike in production.

The confusion is understandable because both speak SQL, both call into GEOS for geometry predicates, and both can read a GeoParquet file. But one is a row-oriented transactional server built to serve thousands of concurrent point reads and writes, and the other is an embedded columnar engine built to scan a column of ten million geometries as fast as the disk can feed it. Choosing on syntax similarity rather than storage model is how teams end up running nightly analytical aggregations against a transactional heap that vacuum can never keep ahead of, or trying to serve concurrent web-map writes from a single-writer embedded file that locks on the second connection.

Two Storage Models, Two Load Targets

The architecture below shows why the same SELECT ... ST_Intersects(...) query has a completely different execution profile on each engine. PostGIS reads rows through a buffer cache and narrows the scan with a GiST index; DuckDB reads columns directly from GeoParquet, skipping row groups whose bounding-box statistics fall outside the query window before it decompresses a single geometry.

The rest of this guide works the difference into a decision you can defend, with a compatibility matrix, an indicative benchmark, and the load code for each path. The companion recipe on loading GeoDataFrames into PostGIS with GeoPandas covers the transactional write path in full once you have chosen PostGIS.

Prerequisites & Environment

You need Python 3.10+, both engines installed, and a GeoParquet sample large enough that the storage model actually matters — a few hundred rows will make both engines look identical and teach you nothing.

• PostGIS side — a PostgreSQL 15+ server with the PostGIS 3.4+ extension, reached through psycopg (3.x) and SQLAlchemy>=2.0. GeoPandas>=1.0 provides to_postgis.
• DuckDB side — duckdb>=1.1 with the bundled spatial extension (installed at runtime with INSTALL spatial; LOAD spatial;). No server, no daemon, no connection pool.
• Shared — GeoPandas>=1.0, shapely>=2.0 (vectorized array API), pyarrow>=15 for GeoParquet I/O.

pip install "geopandas>=1.0" "shapely>=2.0" "pyarrow>=15" \
    "duckdb>=1.1" "SQLAlchemy>=2.0" "psycopg[binary]>=3.1" "geoalchemy2>=0.15"

Verify both engines resolve their spatial stack before you benchmark anything — a missing GEOS or an unloaded spatial extension produces confusing “function does not exist” errors that look like data problems.

import duckdb

con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial;")
print(con.execute("SELECT ST_AsText(ST_Point(0, 0))").fetchone())  # -> ('POINT (0 0)',)
# PostGIS check: SELECT postgis_full_version(); over your SQLAlchemy engine

Version / Compatibility Matrix

This is the comparison matrix. Read each row as “what does this engine make cheap, and what does it make expensive?” — the answers are near mirror images, which is exactly why the two engines rarely compete for the same workload once you look past the shared SQL surface.

The single most important row is the first one. Everything else — indexing, concurrency, scan speed — follows from whether the engine stores data by row or by column, and whether it owns a running server or lives inside your worker process.

Step-by-Step Implementation

The fair way to compare load targets is to run the same analytical query — a spatial aggregation over a moderately large feature set — through each engine’s idiomatic path, then measure. The steps below build that harness.

Step 1 — Materialize a representative GeoParquet source

Both engines should start from the same on-disk GeoParquet so the comparison measures the engines, not the extraction. Write it once from a GeoDataFrame.

from pathlib import Path

import geopandas as gpd


def write_source_geoparquet(gdf: gpd.GeoDataFrame, path: Path) -> Path:
    """Persist a GeoDataFrame as GeoParquet with row-group sizing tuned for scans."""
    gdf = gdf.to_crs(4326)  # normalize CRS so both engines agree on the SRID
    path.parent.mkdir(parents=True, exist_ok=True)
    # Smaller row groups improve DuckDB's bbox-based skipping on selective queries
    gdf.to_parquet(path, geometry_encoding="WKB", row_group_size=100_000)
    return path

Choosing GeoParquet as the interchange format is a pipeline decision in its own right; the tradeoffs against GeoPackage and FlatGeobuf are worked through in choosing a vector output format for spatial pipelines.

Step 2 — Load into PostGIS and build the GiST index

PostGIS cannot scan the GeoParquet in place; it must ingest rows into the heap first, and the analytical scan is only fast once a GiST index exists. Building the index after the bulk insert is materially faster than maintaining it row by row.

import geopandas as gpd
from sqlalchemy import Engine, text


def load_into_postgis(gdf: gpd.GeoDataFrame, engine: Engine, table: str) -> None:
    """Bulk-load a GeoDataFrame into PostGIS, then build the GiST index once."""
    gdf.to_postgis(table, engine, if_exists="replace", index=False, chunksize=50_000)
    with engine.begin() as conn:
        conn.execute(
            text(f'CREATE INDEX IF NOT EXISTS "{table}_geom_gix" '
                 f'ON "{table}" USING GIST (geometry)')
        )
        conn.execute(text(f'ANALYZE "{table}"'))  # refresh planner statistics

Step 3 — Register the same source in DuckDB with no import

DuckDB’s idiomatic path is the opposite: point a view at the GeoParquet and let the engine scan it. There is no load step and no index to build for a first analytical pass.

import duckdb


def register_geoparquet(con: duckdb.DuckDBPyConnection, path: str, view: str) -> None:
    """Expose a GeoParquet file to DuckDB as a scannable view — no import."""
    con.execute("INSTALL spatial; LOAD spatial;")
    con.execute(
        f"CREATE OR REPLACE VIEW {view} AS "
        f"SELECT * FROM read_parquet(?)",
        [path],
    )

Step 4 — Run the identical analytical query on both

The workload is a spatial aggregation: count features within a bounding window and summarize an attribute. Expressed almost identically in both dialects, it exercises each engine’s scan-and-filter path.

BBOX = (-74.05, 40.60, -73.85, 40.90)  # a query window in WGS84


def scan_postgis(engine, table: str, bbox) -> tuple:
    from sqlalchemy import text
    sql = text(
        f'SELECT count(*), avg(ST_Area(geometry)) FROM "{table}" '
        f"WHERE geometry && ST_MakeEnvelope(:x0, :y0, :x1, :y1, 4326)"
    )
    with engine.connect() as conn:
        return conn.execute(sql, dict(zip(("x0", "y0", "x1", "y1"), bbox))).one()


def scan_duckdb(con, view: str, bbox) -> tuple:
    sql = (
        f"SELECT count(*), avg(ST_Area(geometry)) FROM {view} "
        f"WHERE ST_Intersects(geometry, ST_MakeEnvelope(?, ?, ?, ?))"
    )
    return con.execute(sql, list(bbox)).fetchone()

The && bounding-box operator in PostGIS is what lets the planner use the GiST index; DuckDB reaches the same selectivity by skipping row groups whose stored bbox misses the window.

Advanced Patterns & Edge Cases

Incremental loads change the calculus

A one-shot benchmark rewards DuckDB’s zero-import scan, but production pipelines load deltas on a schedule, and there the transactional engine often wins on total cost. PostGIS can upsert only the changed features in one transaction; DuckDB’s columnar files are immutable, so an incremental update means rewriting whole partitions. If your pipeline appends a daily window rather than reloading the world, read scheduling and incremental spatial loads before you commit to either engine — the write pattern, not the scan, dominates the bill.

Concurrent writers and the single-writer ceiling

DuckDB permits exactly one writer per database file. That is a non-issue for a batch transform that opens the file, writes, and closes, but it is fatal for a serving layer where a web map and a nightly load both write. PostGIS’s MVCC handles that mix natively. Do not discover this at 2 a.m. when the second connection blocks — decide up front whether the load target is also a serving target.

Geometry function coverage gaps

PostGIS remains the broadest geometry surface: topology, raster, ST_Subdivide for splitting monster polygons, and the <-> KNN operator with index support. DuckDB spatial covers the everyday predicates and measures but has thinner topology and raster support. If your load target must also serve nearest-neighbor queries or topological validation, the coverage row of the matrix, not the scan speed, decides it.

The hybrid pattern: DuckDB scans, PostGIS serves

The two engines are not mutually exclusive, and the strongest production architectures use both for what each does best. A common shape is to let DuckDB do the heavy analytical passes directly over a partitioned GeoParquet lake — building tiled aggregates, running one-off coverage studies, computing spatial joins across the full history — and to load only the small, curated result into PostGIS for the serving layer that a web map or an API queries with concurrent writes. DuckDB can even read from and write to PostGIS through its postgres extension, so a transform task can ATTACH a live PostGIS database, scan a GeoParquet file, and upsert a summarized layer in one SQL script without moving data through Python. Treated this way, “PostGIS vs DuckDB” stops being an either/or and becomes a division of labour: the columnar engine owns the scan-heavy batch analytics, and the transactional engine owns the concurrent, low-latency serving surface.

Performance Optimization

The table below is an indicative benchmark, not a guarantee — a single-node run over a synthetic 10-million-feature polygon dataset in GeoParquet, executing the Step 4 bounding-box aggregation. Absolute numbers depend entirely on hardware, row-group sizing, geometry complexity, and cache state; treat the shape of the result as the signal, not the milliseconds.

Two patterns emerge and both follow from the storage model. On a cold, full-scan analytical query the columnar engine pulls ahead because it reads only the columns the query touches and decompresses geometry lazily. On the incremental delta the transactional engine pulls ahead because it mutates only the affected rows instead of rewriting a partition. Optimize accordingly: give DuckDB smaller row groups and a bbox-friendly query, and give PostGIS a staged upsert plus a post-load ANALYZE so the planner keeps choosing the GiST index.

import time
from contextlib import contextmanager


@contextmanager
def timed(label: str):
    """Minimal wall-clock timer for A/B load-target benchmarking."""
    start = time.perf_counter()
    yield
    print(f"{label}: {time.perf_counter() - start:.3f}s")


with timed("duckdb_full_scan"):
    scan_duckdb(con, "features", (-180, -90, 180, 90))
with timed("postgis_full_scan"):
    scan_postgis(engine, "features", (-180, -90, 180, 90))

Integration into ETL Pipelines

The load target is not an isolated choice; it dictates the shape of the Load task in your orchestrator. If you pick PostGIS, the Load stage is a transactional upsert keyed on a stable identifier — the exact pattern the loading GeoDataFrames into PostGIS with GeoPandas recipe implements with a staging table and a GiST index. If you pick DuckDB, the Load stage is often “write a partitioned GeoParquet and register it,” which pushes the durability question onto object storage rather than a server.

Either way, keep the interchange format stable between stages: emit GeoParquet from the Transform step, and let the Load step decide whether to import it (PostGIS) or scan it in place (DuckDB). That keeps the two engines swappable behind one contract, so a pipeline that outgrows an embedded scan can graduate to a transactional server without rewriting anything upstream of Load. Route features that fail geometry or schema gates to a dead-letter store before the Load task, so neither engine ever ingests an invalid geometry that would abort a transaction or poison a scan.

Troubleshooting Reference

Related

• Loading GeoDataFrames into PostGIS with GeoPandas — the transactional write path with SRID enforcement, chunked inserts, staging-table upserts, and GiST indexing
• Scheduling and Incremental Spatial Loads — watermark tracking and delta loads, where the write pattern decides between the two engines
• Choosing a Vector Output Format for Spatial Pipelines — why GeoParquet is the interchange format that keeps load targets swappable
• Orchestrating Spatial ETL Pipelines — the top-level reference covering how the Load stage fits the wider DAG