Querying GeoParquet with DuckDB Spatial

This guide runs spatial SQL straight against a GeoParquet file — local or on object storage — and returns a GeoDataFrame, all without a database server or a full in-memory load. It is for analysts whose data has outgrown a comfortable GeoDataFrame but who don't want to stand up PostGIS. It sits under DuckDB Spatial Analytics in Spatial Analysis & Advanced Query Techniques.

Why This Approach / What Goes Wrong

GeoParquet stores geometry as well-known binary (WKB) in a columnar layout, with per-row-group bounding-box statistics written into the file footer. DuckDB exploits both facts at once. Because the format is columnar, a query that touches only building_id and geom never pays to decode the other attributes in the file. And because DuckDB reads the row-group bbox statistics, it can skip whole chunks whose extent cannot satisfy a spatial predicate — so a windowed query over a multi-gigabyte file reads only a fraction of the bytes on disk. This is predicate pushdown and row-group pruning working together, and it is why an in-process engine can out-run a naive full-file read into GeoPandas.

Column pruning and row-group pruning on a local GeoParquet file A local GeoParquet file is drawn as a grid: five attribute columns across, five row groups down, each row group carrying bounding-box statistics. A query selecting only building_id, area_m2 and geom lets DuckDB skip the height and class columns entirely (column pruning), while a spatial ST_Intersects window overlaps only row groups 2 and 3, so the other three row groups are skipped (row-group pruning). Only the cells at the intersection — the shaded blue cells — are read from disk, decoded by the spatial extension, encoded with ST_AsWKB and rebuilt into a GeoDataFrame with the CRS set to EPSG:4326. Only a fraction of the bytes on disk are ever touched. Two prunings at once: only matching columns × matching row groups are read Local GeoParquet file columnar layout · per-row-group bbox statistics building_id area_m2 geom height class RG 1 · skip RG 2 · hit RG 3 · hit RG 4 · skip RG 5 · skip read + decoded bbox-skipped column pruned only the shaded cells leave the disk — column pruning × row-group pruning matching rows only DuckDB spatial extension ST_Intersects bbox filter ST_AsWKB GeoDataFrame geometry from WKB CRS = EPSG:4326 a fraction of the bytes on disk
Column pruning and row-group pruning on a local GeoParquet file

The naive alternative — gpd.read_parquet() — loads every row and every column into memory before you filter, which defeats the whole point of a columnar format on a large file. The failure modes when you move to DuckDB are almost all mechanical:

Pass geometry across the boundary as WKB and set the CRS by hand, and the hand-off is clean and lossless.

There is a fourth failure mode that costs more than the other three combined, because it produces correct results and simply refuses to be fast: the file may have nothing to prune with. Row-group skipping works on ordinary Parquet column statistics — the min and max recorded for each column in each row group — and the minimum and maximum of a WKB byte string are meaningless as a spatial extent. Nothing about ST_Intersects can be pushed into a raw geometry column on its own. Pruning only happens when the file also carries per-row-group extent information: either a GeoParquet 1.1 bbox covering column (a struct of xmin, ymin, xmax, ymax whose numeric statistics are usable), or an equivalent set of plain numeric columns you added yourself. Files written by older tooling, or by a Parquet round-trip that dropped the extra columns, have neither — and a spatial filter over them reads every byte in the file while looking exactly like a query that should be selective. Step 7 shows how to tell the two cases apart in ten seconds, and step 8 how to fix a file that lacks the machinery.

Prerequisites

pip install "duckdb>=0.10" "geopandas>=0.14" "shapely>=2.0"

Version differences change what this guide's queries do, so it is worth knowing which side of them you are on. DuckDB 0.10 will read GeoParquet, but it binds the geometry column as a plain BLOB and every ST_* reference needs an explicit ST_GeomFromWKB(geom) wrapper; from 1.1 onward the reader recognises the geo footer metadata and hands you a real GEOMETRY, so the wrapper becomes unnecessary and the examples below work verbatim. ST_Hilbert, used in step 8, arrived with the 1.1-era extension. On the Python side, GeoPandas 1.0 moved GeoParquet I/O onto pyogrio and made GeoSeries.from_wkb significantly faster on large arrays, which matters once the result set is in the millions of rows. If you can choose, duckdb>=1.1 with geopandas>=1.0 is the pairing that needs the fewest workarounds:

pip install "duckdb>=1.1" "geopandas>=1.0" "shapely>=2.0" "pyarrow>=14"

Step-by-Step Implementation

1. Connect and load the spatial extension.

DuckDB opens an in-memory database by default; pass a path to connect() if you want the connection to persist cached extensions and temp tables.

import duckdb

con = duckdb.connect()          # in-memory; pass "analytics.duckdb" to persist
con.install_extension("spatial")
con.load_extension("spatial")   # required in every fresh session

2. Inspect the file's schema and geometry column name.

Different writers name the geometry column geometry, geom, or wkb_geometry. Check before you write the query rather than guessing.

schema = con.sql("DESCRIBE SELECT * FROM 'buildings.parquet'").df()
print(schema[["column_name", "column_type"]].to_string(index=False))
# column_name column_type
#  building_id      BIGINT
#     area_m2      DOUBLE
#        geom    GEOMETRY

3. Run a windowed spatial query, returning geometry as WKB.

ST_MakeEnvelope builds the query window from lon/lat corners, and ST_Intersects is the predicate DuckDB can push down against the row-group bounding boxes. Encoding the result geometry with ST_AsWKB keeps it in the one format GeoPandas can decode without ambiguity.

query = """
SELECT building_id, area_m2, ST_AsWKB(geom) AS wkb
FROM 'buildings.parquet'
WHERE ST_Intersects(
    geom,
    ST_MakeEnvelope(7.60, 45.00, 7.80, 45.10)   -- lon/lat window, EPSG:4326 (Turin)
)
AND area_m2 > 120
"""
result_df = con.sql(query).df()

4. Rebuild a GeoDataFrame, setting the CRS explicitly.

Decode the WKB column into a GeoSeries, drop the raw bytes, and stamp the CRS the file was authored in. The CRS is not carried in the WKB, so this line is mandatory, not optional.

Where the CRS is lost and where it is put back Four stages track one geometry column across the boundary. In the GeoParquet file the CRS lives in the geo file metadata, so it is present. Inside DuckDB the value is a GEOMETRY of planar coordinates with no CRS object attached, so the CRS is dropped. ST_AsWKB produces bytes in a plain columnar table, still with no datum. Only the final step restores it, because GeoSeries.from_wkb is given an explicit crs argument. The CRS never travels with the geometry itself. The CRS does not ride along with the geometry GeoParquet file metadata geo → EPSG:4326 DuckDB GEOMETRY planar coordinates no CRS object WKB bytes plain columnar table ST_AsWKB(geom) GeoDataFrame from_wkb + crs= crs="EPSG:4326" Is the coordinate reference system still attached? yes held in the footer no coordinates only no bytes carry no datum restored by hand, in step 4 Skip the crs= argument and the frame is silently unprojected for every downstream operation.
The CRS lives in the file footer and in the GeoDataFrame, never in between — the crs= argument in step 4 is the only thing that reconnects the two ends.
import geopandas as gpd

buildings = gpd.GeoDataFrame(
    result_df.drop(columns="wkb"),
    geometry=gpd.GeoSeries.from_wkb(result_df["wkb"]),
    crs="EPSG:4326",          # the CRS the GeoParquet was written in
)

5. Compute metric quantities in a projected CRS, not degrees.

If you need real areas or distances, reproject inside SQL with ST_Transform before measuring. For the Turin window, UTM zone 32N (EPSG:32632) is the correct metric CRS — never use Web Mercator for area or length. Note that ST_Transform in DuckDB takes source and target CRS as arguments and, like PyProj, expects lon/lat axis order for EPSG:4326 input here.

metric_sql = """
SELECT
    building_id,
    ST_Area(ST_Transform(geom, 'EPSG:4326', 'EPSG:32632')) AS area_m2_true
FROM 'buildings.parquet'
WHERE ST_Intersects(geom, ST_MakeEnvelope(7.60, 45.00, 7.80, 45.10))
"""
areas = con.sql(metric_sql).df()
print(areas["area_m2_true"].describe())   # values in square metres, not degrees

6. (Optional) Query remote GeoParquet over HTTP range requests.

The httpfs extension lets DuckDB read a remote file with range requests, pulling only the row groups the predicate needs. Keep the filter spatial so pruning still applies over the network.

Whole-file fetch compared with ranged reads of a remote GeoParquet file Two bars represent the same 2.1 gigabyte remote parcels file. Without range support the whole bar is transferred: every row group crosses the network before any filter runs. With httpfs loaded and a spatial predicate, the file is drawn as twenty-four row-group cells; only the three whose bounding boxes intersect the query window are fetched, plus the footer that holds the statistics. The result is roughly forty megabytes on the wire instead of 2.1 gigabytes. Range requests fetch row groups, not whole files Full fetch no range support 2.1 GB across the wire — every row group, then filter locally Ranged reads httpfs + ST_Intersects 3 row groups whose bbox hits the window footer: bbox statistics Same query, roughly 40 MB on the wire instead of 2.1 GB
Over HTTP the pruning becomes a bandwidth saving: DuckDB reads the footer statistics, then requests byte ranges for only the row groups the spatial predicate can match.
con.install_extension("httpfs")
con.load_extension("httpfs")

remote_count = con.sql(
    "SELECT count(*) FROM 'https://data.example.com/parcels.parquet' "
    "WHERE ST_Intersects(geom, ST_MakeEnvelope(7.60, 45.00, 7.80, 45.10))"
).fetchone()[0]
print("Matching parcels:", remote_count)

7. Prove that pruning actually happened.

Wall-clock time is a poor test — a fast disk hides a full scan of a 2 GB file. Read the row counts instead. Profiling prints how many rows each operator emitted, and a PARQUET_SCAN node that emitted the file's entire row count did no pruning at all, no matter how selective the WHERE clause looks.

con.execute("PRAGMA enable_profiling='query_tree'")
con.sql("""
    SELECT count(*) FROM 'buildings.parquet'
    WHERE ST_Intersects(geom, ST_MakeEnvelope(7.60, 45.00, 7.80, 45.10))
""").fetchall()
con.execute("PRAGMA disable_profiling")

# PARQUET_SCAN  ... Rows: 41,208     <- pruned: only matching row groups decoded
# PARQUET_SCAN  ... Rows: 8,914,330  <- NOT pruned: whole file scanned, then filtered

The metadata function tells you the same thing without running the query, and it is the faster diagnosis when you are deciding whether a file is worth querying remotely at all. If parquet_metadata reports one enormous row group, or parquet_schema shows no bbox column, there is nothing to skip:

rg = con.sql("""
    SELECT row_group_id, row_group_num_rows
    FROM parquet_metadata('buildings.parquet')
    GROUP BY 1, 2 ORDER BY 1 LIMIT 5
""").df()
print(rg)
# One row group of 8.9M rows means row-group pruning can never help.
# Aim for roughly 100k–1M rows per group on a spatially sorted file.

cols = con.sql("SELECT name FROM parquet_schema('buildings.parquet')").df()
print("bbox covering column present:", cols["name"].str.startswith("bbox").any())

8. Make a slow file prunable: sort spatially, then add a bbox column.

Pruning is only as good as the clustering of the data. If the rows arrived in database insertion order, every row group spans the whole study area, every bounding box overlaps every query window, and nothing can be skipped even with perfect statistics. Sorting by a space-filling curve fixes that: ST_Hilbert maps each geometry's centre to a single integer such that points close in space land close in the ordering, so consecutive rows — and therefore row groups — end up spatially compact. Write out the bbox components alongside it so the numeric statistics exist for the planner to use.

import duckdb

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

bounds = con.sql(
    "SELECT ST_Extent_Agg(geom) FROM 'buildings.parquet'"
).fetchone()[0]

con.execute("""
    COPY (
        SELECT *,
               ST_XMin(geom) AS bbox_xmin, ST_YMin(geom) AS bbox_ymin,
               ST_XMax(geom) AS bbox_xmax, ST_YMax(geom) AS bbox_ymax
        FROM 'buildings.parquet'
        ORDER BY ST_Hilbert(geom, ?::BOX_2D)
    ) TO 'buildings_sorted.parquet'
      (FORMAT GDAL, DRIVER 'Parquet', SRS 'EPSG:4326', ROW_GROUP_SIZE 200000)
""", [bounds])

Query the sorted copy with the numeric columns in the predicate first and the exact test second. The four inequalities are ordinary column comparisons, so the planner prunes row groups on them, and ST_Intersects then runs only on what survived:

fast = con.sql("""
    SELECT building_id, ST_AsWKB(geom) AS wkb
    FROM 'buildings_sorted.parquet'
    WHERE bbox_xmin <= 7.80 AND bbox_xmax >= 7.60
      AND bbox_ymin <= 45.10 AND bbox_ymax >= 45.00
      AND ST_Intersects(geom, ST_MakeEnvelope(7.60, 45.00, 7.80, 45.10))
""").df()
print("rows:", len(fast))

On a nationwide building footprint file the difference is not marginal. An unsorted 9-million-row export answers a city-sized window by decoding all nine million geometries; the Hilbert-sorted copy touches four or five row groups and returns in a second or two. The one-time sort costs a few minutes and is worth it for any file you will query more than twice — the same reasoning that justifies pre-tiling in Generating PMTiles from GeoParquet.

Verification

Confirm the result decoded to real geometries with the expected CRS, and that the spatial filter actually narrowed the set rather than silently returning everything.

print("Rows returned:", len(buildings))            # Rows returned: 4127
print("CRS:", buildings.crs.to_epsg())             # CRS: 4326
print("Geom type:", buildings.geometry.geom_type.unique())  # ['Polygon']
assert buildings.geometry.notnull().all()

# All features fall inside the requested window
minx, miny, maxx, maxy = buildings.total_bounds
assert 7.60 <= minx and maxx <= 7.80, "Filter window not respected"

Edge Cases & Debugging

Frequently Asked Questions

Should I query GeoParquet with DuckDB or just use gpd.read_parquet with a filter? gpd.read_parquet accepts filters that prune on ordinary attribute columns, so if your selection is land_use == 'residential' and the file fits in RAM, GeoPandas is simpler and there is no reason to add an engine. Reach for DuckDB when the selection is spatial, when the file is larger than memory, when it lives on object storage, or when the result you actually want is an aggregate — in that last case DuckDB returns a few hundred rows where GeoPandas would have to materialize millions first.

Does DuckDB read the geometry faster than GeoPandas? Per geometry, no — both end up decoding WKB, and GeoPandas 1.0 with pyogrio is competitive at that. The advantage is entirely in how many geometries get decoded. A well-sorted file with row-group statistics lets DuckDB decode one percent of the file; GeoPandas decodes everything it reads. That is why step 8 matters more than any engine choice: on a badly laid-out file, DuckDB's advantage largely disappears.

Why does the CRS have to be set by hand every single time? Because WKB is a coordinate encoding, not a spatial reference. The CRS lives in the GeoParquet footer, DuckDB's GEOMETRY type has no field to carry it, and the result set that reaches Python is a plain columnar table. Read the file's declared CRS once with gpd.read_parquet(path, columns=[]).crs or from the geo metadata and store it in a constant, rather than hard-coding EPSG:4326 and hoping — the second habit is what produces silently unprojected frames.

Can I write the query result straight back to GeoParquet without touching Python? Yes, and you should when the result is large. COPY (SELECT ...) TO 'out.parquet' (FORMAT GDAL, DRIVER 'Parquet', SRS 'EPSG:4326') keeps the whole round trip inside the engine, never materializes a frame, and produces a file GeoPandas can read directly. Only pull the result into Python when Python is going to do something with it.

Is it safe to run this against a file another job is writing? No. Parquet has no in-place update path, so a writer that is rewriting a file mid-query gives DuckDB a truncated or inconsistent footer, usually surfacing as a decode error rather than wrong answers. Write to a new path and swap, or use a partitioned directory where each write lands a new partition — the same discipline the Cloud-Native Geospatial Formats workflow expects.

Related