DuckDB Spatial Analytics in Python: In-Process Spatial SQL

DuckDB brings database-grade spatial SQL into the Python process itself — no server, no connection pool, just an in-process engine that reads GeoParquet, Shapefiles, and GeoJSON through its spatial extension. For analytical workloads that are too big for a comfortable GeoDataFrame but don't justify a PostGIS deployment, it is often the fastest path. This stage of Spatial Analysis & Advanced Query Techniques covers the spatial extension end to end, and sits beside the two other out-of-core engines in this section: PostGIS Integration with Python for indexed, shared, transactional storage, and Scaling with Dask-GeoPandas for distributed, partitioned processing across many cores or machines.

DuckDB in-process spatial query DuckDB reads GeoParquet directly from disk or object storage, runs vectorized spatial SQL in-process, and returns results to GeoPandas without a server. GeoParquet local / S3 / R2 columnar DuckDB (in-process) spatial extension ST_* + R-tree vectorized, columnar predicate pushdown GeoDataFrame result rows no server
DuckDB runs columnar spatial SQL inside your Python process, reading GeoParquet directly — including remotely over HTTP range requests.

Architecture & Data Structures

DuckDB is a columnar, vectorized OLAP engine that runs embedded in your process, like SQLite but built for analytics. It executes queries a morsel at a time — small batches of column values pushed through the operator pipeline — so a spatial predicate runs across thousands of geometries per SIMD-friendly loop rather than row by row. The spatial extension layers geospatial capability on top of that engine:

Two properties make it fast on real data. Because it is columnar, it reads only the columns a query touches — a query selecting parcel_id and geom never pays to decode the other forty attributes in the file. And because it understands GeoParquet's per-row-group bounding-box statistics, it can skip whole chunks whose extent cannot satisfy a spatial filter, so a windowed query over a multi-gigabyte file touches only a fraction of it. For repeated in-memory lookups you can also build an explicit R-tree with CREATE INDEX ... USING RTREE (geom).

import duckdb

con = duckdb.connect()          # in-memory database; pass a path to persist
con.install_extension("spatial") # one-time download into the extension cache
con.load_extension("spatial")    # load into this connection

# Query a GeoParquet file directly — no load or import step
con.sql("""
    SELECT count(*) AS large_footprints
    FROM 'buildings.parquet'
    WHERE ST_Area(ST_Transform(geom, 'EPSG:4326', 'EPSG:32632')) > 500  -- m²
""").show()

Three concrete column types turn up in real result sets, and knowing which one you are holding saves an hour of confused debugging. GEOMETRY is the extension's own type and the only thing the ST_* family accepts directly. WKB_BLOB — or a plain BLOB — is a byte column that looks like geometry but is not: every function call against it fails with a binder error until you wrap it in ST_GeomFromWKB. BOX_2D is a four-double struct used for extents, produced by ST_Extent and its aggregate form ST_Extent_Agg; comparing two BOX_2D values costs a handful of float comparisons against the hundreds of coordinate tests a full ST_Intersects may run, which is exactly the asymmetry the query patterns further down exploit.

Which type a Parquet scan binds depends on the writer and the DuckDB version, so check rather than assume. Files written by GDAL or GeoPandas store geometry as WKB with geo metadata in the footer; recent DuckDB releases recognise that metadata and hand you a GEOMETRY, while older ones — and any file whose geo metadata was stripped by a Parquet rewrite that did not preserve key-value metadata — give you a BLOB. A one-line DESCRIBE settles it before you write the rest of the query:

-- Which type did the scan actually bind?
DESCRIBE SELECT * FROM 'parcels.parquet';
-- geom  GEOMETRY   -> ST_Intersects(geom, ...) works as written
-- geom  BLOB       -> wrap every reference: ST_GeomFromWKB(geom)

The R-tree deserves the same scepticism. CREATE INDEX ... USING RTREE (geom) only works on a materialized DuckDB table — you cannot index a Parquet glob in place — and the planner only uses it when the predicate compares the indexed column against a constant geometry, such as a literal envelope. A join predicate whose right side is another column will not touch the index, which is why a two-table spatial join falls back to a nested loop no matter how many indexes you build. That constraint shapes the whole scaling story below, and it is the single biggest structural difference from a server-side GiST index; the trade-off in full is laid out in DuckDB Spatial vs PostGIS for Analytics.

Format readers divide along the same line. ST_Read routes through GDAL's OGR layer and materializes the whole layer before the rest of the plan sees a row, so a 4 GB Shapefile costs 4 GB of scan whatever the WHERE clause says — it is an ingestion path, not a query path. Native Parquet and CSV scans are the ones that prune. If a format only arrives through ST_Read, the productive move is to convert once with COPY (SELECT * FROM ST_Read('parcels.shp')) TO 'parcels.parquet' and query the columnar copy from then on, for the reasons set out in GeoParquet vs Shapefile for Storage.

Environment Configuration & Dependency Resolution

pip install "duckdb>=1.1" "geopandas>=1.0" "shapely>=2.0"
# The spatial extension is downloaded on first install_extension("spatial")

DuckDB ships as a single self-contained wheel with no system dependencies — a major reason to choose it on locked-down or Windows machines where compiling the GEOS/GDAL/PROJ stack is painful. That isolation is the opposite trade-off from GeoPandas, whose reproducibility hinges on aligning those C libraries; here the entire engine, including its bundled GDAL and PROJ, arrives in one artifact.

The spatial extension itself is fetched at runtime the first time you call install_extension("spatial"), then cached under ~/.duckdb/extensions. Pin duckdb>=1.1 so the extension's ABI matches the engine — a mismatched pair raises an Invalid Input Error on load. Three environment notes matter in production:

What DuckDB occupies on a single host The Python process holds the whole engine: a buffer manager capped by memory_limit, in which morsels of column values are processed with headroom left over, plus the loaded spatial extension and a hash join that spills once it exceeds the cap. Three locations on disk serve it: the extension cache under the DuckDB home directory, the spill space named by temp_directory, and the data itself, read in place from local GeoParquet or object storage. No server process is involved. Where DuckDB puts things on one host Your Python process Buffer manager · memory_limit = 8GB morsels of column values headroom spatial extension loaded per connection hash join spills past the cap load once spill scan extension cache ~/.duckdb/extensions spill space temp_directory='/var/tmp/duck' data, read in place parcels.parquet · s3://bucket/ The engine, its cache and its spill space all live beside your code — there is no daemon to reach.
Every DuckDB moving part sits on the host running the script: RAM bounded by memory_limit, a cached extension binary, a spill directory for oversized joins, and the data files themselves left where they are.
import duckdb

con = duckdb.connect()
for ext in ("spatial", "httpfs"):
    con.install_extension(ext)
    con.load_extension(ext)

con.execute("SET memory_limit='8GB'")            # cap RAM; spill beyond it
con.execute("SET temp_directory='/var/tmp/duck'") # where large joins spill

Two defaults bite in containers. memory_limit defaults to roughly 80% of what the engine believes the machine has, and on most container runtimes that reading comes from the host, not the cgroup — so a pod limited to 4 GB happily plans for 25 GB and gets killed by the OOM reaper with no DuckDB error to read. threads has the same blind spot: it defaults to the host core count, so a one-CPU pod spawns thirty-two worker threads that spend their time context-switching. Set both explicitly from whatever your orchestrator actually granted:

import os

con.execute(f"SET threads={os.cpu_count()}")     # replace with the cgroup quota in k8s
con.execute("SET preserve_insertion_order=false") # frees the engine to stream large scans

preserve_insertion_order=false is the cheapest memory win available on a big scan-and-write job: with ordering guaranteed, DuckDB must buffer results to reassemble the original row order, and turning it off lets the pipeline stream straight through. Only set it where row order genuinely does not matter — which for an aggregate or an ORDER BY-terminated query it never does.

Connections are the other place people get surprised. A duckdb.connect() object is not safe to share across threads; the supported pattern is one connection per process and con.cursor() per thread, each cursor carrying its own transaction state but sharing the buffer pool. Extensions load per connection, so a fresh cursor inherits the parent's loaded extensions, but a brand-new connect() in a worker process does not — pool workers each need their own load_extension("spatial"). And the module-level duckdb.sql(...) helper quietly operates on a hidden global connection, which is convenient in a notebook and a trap in a library, because settings you applied to your own connection do not apply to it.

Version pinning matters more here than in most Python packages because the extension binary is version-matched to the engine. The spatial extension distributed for DuckDB 1.1 will not load into 1.2, and the failure surfaces as Invalid Input Error rather than anything mentioning versions. In a Docker image, pin the exact DuckDB patch version and bake the extension into the image at build time rather than downloading it on first request:

# Bake the extension into the image so cold starts never hit the network
python -c "import duckdb; c=duckdb.connect(); c.install_extension('spatial'); c.install_extension('httpfs')"

Vectorized Operations & Core Workflow

The core pattern is a three-step boundary crossing: read GeoParquet, filter and reduce with spatial SQL so only the rows you need materialize, then hand the result to a GeoPandas DataFrame for plotting or downstream work. The one detail that makes the hand-off clean is to serialize geometry as WKB at the SQL boundary with ST_AsWKB and rebuild it with GeoSeries.from_wkb, setting the CRS explicitly — DuckDB does not carry a Python-side CRS object across. The full recipe lives in Querying GeoParquet with DuckDB Spatial.

import duckdb
import geopandas as gpd

con = duckdb.connect()
con.load_extension("spatial")

# Spatial filter in SQL; return geometry as WKB for a clean GeoPandas hand-off.
# ST_MakeEnvelope takes (minx, miny, maxx, maxy) in the geometry's own CRS.
df = con.sql("""
    SELECT parcel_id, land_use, ST_AsWKB(geom) AS wkb
    FROM 'parcels.parquet'
    WHERE ST_Intersects(geom, ST_MakeEnvelope(7.6, 45.0, 7.8, 45.1))
""").df()

parcels = gpd.GeoDataFrame(
    df.drop(columns="wkb"),
    geometry=gpd.GeoSeries.from_wkb(df["wkb"]),
    crs="EPSG:4326",  # must match the CRS the parquet geometry is stored in
)

The envelope filter here is a cheap coarse pass — DuckDB uses the row-group bounding boxes to prune before it ever evaluates ST_Intersects on individual geometries. For an exact clip rather than a bounding-box hit, follow up with ST_Intersection in the same query, or express the tighter predicate directly (ST_Contains, ST_DWithin) so the engine still benefits from bbox pruning while returning geometrically correct rows.

The boundary crosses in the other direction too, and it is the half most workflows forget. Anything in your Python session that exposes an Arrow table or a pandas frame — including a GeoDataFrame whose geometry you have serialized to WKB — is queryable by variable name inside a DuckDB SQL string, with no copy and no temporary file. That makes the common "small curated layer joined against a huge file" pattern a one-liner: keep the study-area boundaries in memory where you edited them, and let DuckDB stream the large side off disk.

import duckdb
import geopandas as gpd

con = duckdb.connect()
con.load_extension("spatial")

# A small, hand-curated layer that lives in the notebook
study_areas = gpd.read_file("study_areas.gpkg").to_crs("EPSG:4326")
study_wkb = study_areas.assign(geom_wkb=study_areas.geometry.to_wkb()).drop(columns="geometry")

# Referenced by variable name; DuckDB reads it in place from the Python session
counts = con.sql("""
    SELECT a.area_name, count(*) AS sensor_count
    FROM study_wkb a
    JOIN 'sensors.parquet' s
      ON ST_Contains(ST_GeomFromWKB(a.geom_wkb), s.geom)
    GROUP BY a.area_name
""").df()

For results that are themselves too large to sit in one frame, do not call .df() at all. fetch_record_batch hands back an Arrow record-batch reader that pulls from the running query in chunks, so peak memory is one batch rather than the whole result — the same streaming discipline that makes pyogrio and Arrow fast on the GeoPandas side.

reader = con.sql(
    "SELECT parcel_id, ST_AsWKB(geom) AS wkb FROM 'parcels.parquet'"
).fetch_record_batch(100_000)

total_vertices = 0
for batch in reader:                       # one 100k-row batch resident at a time
    chunk = gpd.GeoSeries.from_wkb(batch.column("wkb").to_pylist(), crs="EPSG:4326")
    total_vertices += chunk.count_coordinates().sum()
print("vertices scanned:", total_vertices)

Geometry / Data Processing Details

DuckDB's spatial functions mirror PostGIS naming, so SQL ports across the two engines with little change — the same ST_Contains, ST_DWithin, ST_Union, and ST_Area calls behave identically because both sit on GEOS. Where DuckDB earns its place is large aggregations that would be memory-heavy in pure GeoPandas: counting points in polygons, dissolving by attribute, computing per-zone statistics, and spatial joins between two files that together exceed RAM. For the predicate and overlay semantics themselves, see the sibling guide on Geometric Intersections & Overlays; DuckDB is the execution engine, not a different definition of the operations.

A point-in-polygon join is the canonical case. Both inputs stream from disk, the join predicate is a spatial containment test, and only the aggregated counts return:

import duckdb

con = duckdb.connect()
con.load_extension("spatial")

# Points-in-polygons aggregation across two GeoParquet files, streamed from disk
con.sql("""
    SELECT z.zone_id, count(*) AS sensor_count
    FROM 'zones.parquet' z
    JOIN 'sensors.parquet' s ON ST_Contains(z.geom, s.geom)
    GROUP BY z.zone_id
    ORDER BY sensor_count DESC
""").show()

A spatial join is a nested loop unless the planner can prune, so on large inputs give it something to prune with. Either build an R-tree on the smaller side, or add a bounding-box pre-filter that lets the planner discard non-overlapping pairs before the expensive ST_Contains:

# Dissolve parcels by land-use class, then keep only sizeable resulting blocks.
# ST_Union_Agg is the aggregate ("dissolve") form; ST_Area runs on projected CRS.
con.sql("""
    WITH dissolved AS (
        SELECT land_use,
               ST_Union_Agg(ST_Transform(geom, 'EPSG:4326', 'EPSG:32632')) AS geom
        FROM 'parcels.parquet'
        GROUP BY land_use
    )
    SELECT land_use, ST_Area(geom) AS area_m2
    FROM dissolved
    WHERE ST_Area(geom) > 10000
    ORDER BY area_m2 DESC
""").show()

Note the aggregate ST_Union_Agg rather than the binary ST_Union — dissolving many geometries into one per group is a reduction, and DuckDB's columnar aggregation is exactly where it beats an equivalent GeoDataFrame.dissolve() on data that no longer fits in memory.

Where the join actually stops scaling. Because the R-tree cannot serve a column-against-column predicate, JOIN ... ON ST_Contains(z.geom, s.geom) is a blockwise nested loop: every zone geometry is tested against every sensor geometry, and the cost is the product of the two row counts multiplied by the per-pair GEOS cost. Two hundred zones against ten million sensors is two billion predicate evaluations — minutes, but survivable. Two hundred thousand zones against the same ten million points is a trillion, and no amount of RAM makes that finish. The inflection is not a data size, it is a product: once left_rows × right_rows passes roughly 10⁹ you need to give the planner an equality or range condition it can hash on.

The reliable way to manufacture one is to materialize the bounding-box components as plain doubles and let the optimizer turn the four inequalities into a range join, keeping the exact predicate as a cheap residual filter over the survivors:

import duckdb

con = duckdb.connect()
con.load_extension("spatial")

# Explode each geometry's envelope into ordinary numeric columns once,
# so the four range conditions can drive the join instead of a nested loop.
con.sql("""
    CREATE TABLE zone_box AS
    SELECT zone_id, geom,
           ST_XMin(geom) AS xmin, ST_XMax(geom) AS xmax,
           ST_YMin(geom) AS ymin, ST_YMax(geom) AS ymax
    FROM 'zones.parquet';

    CREATE TABLE sensor_pt AS
    SELECT sensor_id, geom, ST_X(geom) AS x, ST_Y(geom) AS y
    FROM 'sensors.parquet';
""")

con.sql("""
    SELECT z.zone_id, count(*) AS sensor_count
    FROM zone_box z
    JOIN sensor_pt s
      ON s.x BETWEEN z.xmin AND z.xmax        -- range join: planner-friendly
     AND s.y BETWEEN z.ymin AND z.ymax
    WHERE ST_Contains(z.geom, s.geom)          -- exact test on survivors only
    GROUP BY z.zone_id
""").show()

On elongated or L-shaped zones the envelope is a loose approximation and the residual filter throws away a large share of the candidates, but even a 10:1 false-positive rate leaves the exact test running on a hundredth of the pairs a nested loop would have visited. Watch for the opposite failure: if one layer's features have wildly different sizes — a national boundary sitting in the same table as city blocks — the giant envelopes match nearly everything and the range join degenerates back toward a cross product. Split the layer by feature size and run two joins in that case.

Where the aggregate stops scaling. Hash joins and sorts spill to temp_directory when they exceed memory_limit, but geometry aggregates do not: ST_Union_Agg builds the union of an entire group in memory through GEOS, and there is no spill path for a half-built polygon. Dissolving a million parcels into one landmass will exhaust RAM regardless of the memory limit you set, because the limit governs the buffer manager and not the C++ allocations inside the aggregate. The workaround is hierarchical: union within a coarse spatial key first, materialize, then union the — far fewer, far simpler — intermediate polygons. Dropping the vertex count with ST_Simplify before the union is the other lever, and often the larger one, since union cost tracks vertices rather than rows.

# Two-stage dissolve: union inside coarse grid cells, then union the cell results.
con.sql("""
    WITH cell AS (
        SELECT land_use,
               floor(ST_X(ST_Centroid(geom)) / 0.25) AS cx,
               floor(ST_Y(ST_Centroid(geom)) / 0.25) AS cy,
               ST_Union_Agg(ST_Simplify(geom, 0.00002)) AS geom
        FROM 'parcels.parquet'
        GROUP BY land_use, cx, cy
    )
    SELECT land_use, ST_Union_Agg(geom) AS geom
    FROM cell
    GROUP BY land_use
""").show()

The same staged reduction underpins Dissolving and Aggregating Features by Attribute on the GeoPandas side; DuckDB simply lets each stage read from disk instead of RAM. When even the staged version will not fit on one machine, that is the honest signal to move the job to Scaling with Dask-GeoPandas rather than to buy a larger instance.

CRS Alignment & Projection Pipeline

DuckDB's spatial extension treats coordinates as planar numbers. It does not attach a live coordinate reference system to a GEOMETRY column the way GeoPandas does, and it will not silently reproject between systems in a join — two layers in different CRSs simply produce a join that returns nothing, because their coordinates never overlap numerically. The rule is therefore explicit: confirm every input shares one CRS before you evaluate a spatial predicate, and call ST_Transform (backed by the bundled PROJ) when they don't. Establish and document that canonical CRS upstream with Coordinate Systems with PyProj, then keep it consistent through the whole query.

import duckdb

con = duckdb.connect()
con.load_extension("spatial")

# Reproject 4326 → UTM 32N inside SQL before a metric area calculation.
# ST_Transform(geom, source, target) takes explicit EPSG authority strings.
con.sql("""
    SELECT id, ST_Area(ST_Transform(geom, 'EPSG:4326', 'EPSG:32632')) AS area_m2
    FROM 'fields.parquet'
""").show()

Two gotchas dominate here. First, computing ST_Area or ST_Length on raw EPSG:4326 geometry yields square (or linear) degrees, not metres — a meaningless number that grows or shrinks with latitude. Reproject to an appropriate projected CRS such as the local UTM zone first, exactly as you would in GeoPandas, and never use Web Mercator (EPSG:3857) for measurement. Second, DuckDB's ST_Transform expects axis order in the traditional longitude/latitude sense for EPSG:4326 input, so coordinates that were stored latitude-first will transform to the wrong place; normalize axis order when the data is ingested rather than mid-query. Pick the target zone from the data's centroid so the distortion stays small across the extent.

Axis order is worth pinning down explicitly rather than trusting the default. ST_Transform accepts an optional fourth argument that forces the traditional longitude-then-latitude interpretation, and passing it removes any dependence on how the bundled PROJ resolved the CRS definition's authority-declared axis order:

-- Force lon/lat interpretation regardless of the CRS definition's declared axis order
SELECT ST_Transform(geom, 'EPSG:4326', 'EPSG:32632', always_xy := true) AS geom_utm
FROM 'fields.parquet';

The symptom of getting this wrong is unmistakable once you know it: transformed coordinates land in the Gulf of Guinea, or a metre-scale grid returns values off by hundreds of kilometres, because latitude was fed in where longitude was expected. The same rule and the same failure appear in the PyProj Transformer API, and the reasoning behind both is set out in EPSG vs PROJ String vs WKT CRS Formats.

There is one accuracy difference from PyProj that surprises people migrating a pipeline. The spatial extension ships its own PROJ database inside the extension binary, and that bundle contains the CRS definitions but not the optional datum-shift grid files that PROJ downloads separately. Transformations that need a grid — NAD27 to NAD83 in North America, OSGB36 to ETRS89 in Britain, most national realizations of a shifting datum — silently fall back to a lower-accuracy Helmert approximation instead of raising. The error is metre-scale, not kilometre-scale, so it will not look wrong on a map; it will just quietly disagree with the same transformation done in PyProj with grid files installed. If your work is cadastral, survey-grade, or has to reconcile against an authority dataset, do the datum change in PyProj and use DuckDB only for the same-datum projection step.

Transform cost is worth budgeting for as well. ST_Transform is evaluated per row, and unlike the predicates it cannot be pruned away — a query that transforms before it filters pays the projection cost on every row in the file. Filter first in the source CRS, then transform the survivors:

-- Filter in the stored CRS, transform only what survives
SELECT parcel_id,
       ST_Area(ST_Transform(geom, 'EPSG:4326', 'EPSG:32632', always_xy := true)) AS area_m2
FROM 'parcels.parquet'
WHERE ST_Intersects(geom, ST_MakeEnvelope(7.60, 45.00, 7.80, 45.10));

On a fifty-million-row file that reordering is routinely the difference between a query that returns in seconds and one that runs for several minutes, and it is the same principle that governs Coordinate Reference System Transformations at ingestion time: reproject the smallest set of coordinates that will answer the question.

Production Export & Integration

DuckDB is most valuable as the analytical hop between cloud-native storage and whatever consumes the result — a notebook, a PostGIS table, or a web map. The integration checklist:

One export detail catches almost everyone once. A plain COPY ... TO 'out.parquet' (FORMAT PARQUET) writes a perfectly good Parquet file, but it does not write the geo key-value metadata that makes it GeoParquet — the geometry column lands as an anonymous binary blob, and gpd.read_parquet refuses it or reads it as bytes. Route the write through the GDAL Parquet driver, which records the CRS in the footer, whenever anything other than DuckDB will read the file back; the syntax appears in DuckDB Spatial vs PostGIS for Analytics.

Credentials for object storage belong in a secret rather than in SET statements, which land in query logs and shell history. A persistent secret survives across sessions in a database file, and the CREDENTIAL_CHAIN provider picks up the same instance role or environment credentials your other AWS tooling already uses:

CREATE SECRET parcels_store (
    TYPE S3,
    PROVIDER CREDENTIAL_CHAIN,   -- instance role, env vars, or ~/.aws/credentials
    REGION 'eu-central-1'
);
SELECT count(*) FROM 's3://city-open-data/parcels/*.parquet';

For datasets you will query repeatedly by area, write them out partitioned rather than as one file. Hive-style partitioning on an administrative key puts each partition in its own directory, and DuckDB skips whole directories from the path predicate alone — before it reads a single footer, let alone a geometry.

COPY (
    SELECT *, ST_Transform(geom, 'EPSG:4326', 'EPSG:32632', always_xy := true) AS geom_utm
    FROM 'parcels.parquet'
) TO 'parcels_by_region'
  (FORMAT PARQUET, PARTITION_BY (region_code), OVERWRITE_OR_IGNORE);

-- Later: only the matching directory is ever opened
SELECT count(*) FROM 'parcels_by_region/*/*.parquet' WHERE region_code = 'ITC1';

Partition on a column with tens to low hundreds of distinct values. Partitioning on something high-cardinality — a parcel id, a full postcode — produces thousands of tiny files, and the per-file footer reads then cost more than the scan you were trying to avoid, especially over object storage where each one is a separate round trip. The same file-size arithmetic drives partition sizing in Partitioning Strategies for Dask-GeoPandas.

DuckDB vs PostGIS vs Dask-GeoPandas selection matrix A four-row comparison. Deployment: DuckDB is embedded in-process, PostGIS is a shared server, Dask-GeoPandas runs on a multi-node cluster. Workload: read-mostly analytics, transactional serving, distributed batch respectively. Data location: files and object storage, indexed tables, partitioned across nodes. Best fit: big files on one machine with no server, concurrent indexed queries, and one job spread across many cores. Choosing the engine at a glance Selection axis DuckDB PostGIS Dask-GeoPandas Deployment Embedded in-process, no server Shared server always-on daemon Multi-node cluster of workers Workload Read-mostly analytical queries Transactional concurrent serving Distributed out-of-core batch Data location Files & object storage Indexed tables GiST in the database Partitioned across nodes Best fit Big files, one box, no server to run Concurrent indexed shared queries One job across many cores or nodes
Pick DuckDB for embedded single-machine analytics over files and object storage; PostGIS for shared, indexed, transactional serving; Dask-GeoPandas for a single job fanned out across nodes.

Choosing DuckDB vs the alternatives. Reach for DuckDB when the work is embedded, read-mostly analytics over files and object storage on a single machine — it needs no server and no ingestion step. Move to PostGIS when many clients need shared, transactional, GiST-indexed storage they can query concurrently, and to Dask-GeoPandas when a single job must fan out across partitions on many cores or nodes.

Windows / Platform Edge Cases & Debugging

Frequently Asked Questions

Should I load the data into a DuckDB database file or query the Parquet directly? Query the Parquet directly for anything you run a handful of times — there is no ingest step, the files stay the shared source of truth, and re-running after an upstream refresh needs no reload. Materialize into a .duckdb file only when you need something a file scan cannot give you: an R-tree index, a primary key, repeated joins over the same tables where re-parsing footers dominates, or a persistent secret. The cost of being wrong is small in one direction and large in the other, so start with the files.

Why is my spatial join slower in DuckDB than in GeoPandas? Almost certainly because GeoPandas built an R-tree and DuckDB did not. sjoin constructs an STRtree over the right-hand frame automatically; DuckDB has no equivalent for a column-to-column predicate and runs a nested loop. Below a few hundred thousand features on each side, in-memory GeoPandas usually wins outright. DuckDB pulls ahead when at least one side is too large to hold, or when the join is preceded by a selective filter that pruning can exploit — and it needs the bounding-box range-join rewrite to stay ahead as both sides grow.

Can several processes read the same DuckDB database at once? Not a database file — it takes a single-writer lock and readers of an unlocked file are limited to read-only mode. Parquet files have no such constraint: any number of processes can scan the same GeoParquet concurrently, which is why the file-first pattern also happens to be the concurrency-friendly one. If you need genuinely concurrent readers and writers against shared state, that is the boundary where PostGIS Integration with Python becomes the right answer rather than a heavier one.

Do I still need GeoPandas if DuckDB can do the spatial SQL? Yes, for anything that is not a set-oriented query. Plotting, per-feature Python logic, integration with scikit-learn or spatial clustering algorithms, and the whole CRS object model live on the GeoPandas side. The productive division is DuckDB for the filter-join-aggregate stage that shrinks the data, GeoPandas for everything you do with the result once it fits comfortably in memory.

How do I know a spatial filter actually pruned anything? Turn on profiling and read the row counts, not the wall clock — PRAGMA enable_profiling='query_tree' prints the rows each operator emitted, and a Parquet scan that emitted the full file's row count did no pruning at all. The usual cause is a file written without per-row-group bounding-box statistics, and the fix is to rewrite it; the diagnosis is walked through in Querying GeoParquet with DuckDB Spatial.

Is DuckDB's GEOMETRY type the same thing as PostGIS geometry? The function names match and both sit on GEOS, so predicate and overlay semantics agree, but the type systems differ in one consequential way: PostGIS geometry carries an SRID, DuckDB's does not. There is no spatial_ref_sys table and no way to ask a column what CRS it is in, so the CRS is something you track outside the engine — in file metadata, in a naming convention, or in the code that reads the result. Every CRS bug in a DuckDB pipeline traces back to that missing field.

Related