How to Move Data Between DuckDB and GeoPandas

On this page

Problem statement

DuckDB is the better engine for the heavy part — the join, the aggregate, the scan over millions of rows. GeoPandas is the better environment for everything after that: plotting, .explore(), contextily, the shapely API, and the rest of the Python geospatial ecosystem.

So most real workflows use both, and the boundary between them needs to be cheap and lossless. Two things routinely go wrong:

The geometry arrives as bytes. A DuckDB GEOMETRY column comes into pandas as an opaque object, not as a GeoSeries.
The CRS is lost. DuckDB tracks CRS on the column type and drops it as soon as anything is computed, so a GeoDataFrame built from a query result has crs=None unless you set it.

Both are one line to fix, and both are silent when you do not.

Quick answer

Convert to WKB in SQL, load it in GeoPandas, and set the CRS explicitly:

import duckdb
import geopandas as gpd
from shapely import from_wkb

con = duckdb.connect()
con.execute("install spatial; load spatial;")

frame = con.execute("""
    select name, admin, st_aswkb(geom) as geometry
    from st_read('provinces.shp')
    where admin = 'United States of America'
""").df()

gdf = gpd.GeoDataFrame(
    frame.drop(columns="geometry"),
    geometry=from_wkb(frame["geometry"].map(bytes)),   # DuckDB hands BLOBs over as bytearray
    crs="EPSG:4326",              # DuckDB will not tell you this — you must
)
print(len(gdf), gdf.crs, gdf.geom_type.value_counts().to_dict())

51 EPSG:4326 {'Polygon': 30, 'MultiPolygon': 21}

Going the other way is the mirror image:

con.execute("create table areas as select * from gdf_with_wkb")

with the geometry converted to WKB in pandas first.

Each crossing costs a serialisation, so there should be exactly one.

Step-by-step solution

1. Use WKB as the interchange format

Well-Known Binary is compact, exact and understood by both sides. ST_AsWKB in DuckDB produces it; shapely.from_wkb reads it, vectorised over a whole column. One wrinkle, measured with DuckDB 1.5.5: .df() delivers a BLOB column as Python bytearray objects, and from_wkb rejects them with TypeError: Expected bytes or string, got bytearray. .map(bytes) converts them; reading through Arrow avoids the problem, because Arrow's binary column converts to bytes.

The alternatives are worse: WKT is text and larger, GeoJSON is text and larger still, and passing DuckDB's internal geometry type into pandas gives you something that is not a shapely object.

2. Set the CRS yourself, every time

DuckDB attaches the source CRS to the column type when reading a file:

select typeof(geom) from st_read('provinces.shp') limit 1;
-- GEOMETRY('EPSG:4326')

but that annotation does not travel through ST_AsWKB, and ST_Transform strips it from its own result. So the CRS is knowledge you hold, not something the pipeline carries.

Write it down as a constant next to the query. A GeoDataFrame with crs=None will silently refuse to reproject and silently mis-join.

3. Prefer Arrow for large results

.df() builds a pandas DataFrame; .arrow() returns an Arrow table with no intermediate copy. For results of millions of rows the difference is measurable, and GeoPandas reads from Arrow directly:

sql = "select name, lat, lon from read_parquet('places.parquet') where country = 'GB'"
table = con.execute(sql).arrow()
if hasattr(table, "read_all"):                 # newer DuckDB returns a record batch reader
    table = table.read_all()
frame = table.to_pandas()
print(f"{table.num_rows:,} rows via Arrow")

109,179 rows via Arrow

Compare the transfer costs for a million two-column rows: JSON is 36.25 MB, 6.82 MB gzipped, Arrow with zstd 11.20 MB, raw float32 8.00 MB (Example 4) — on the wire, compression matters more than format. Inside a process there is nothing to compress and the format decides: fetching all 13.5 million GeoNames rows took 4.65 s with .df() and 1.21 s with .arrow().

4. Go the other way with a registered DataFrame

DuckDB can query a pandas DataFrame in the local scope by name, which makes pushing data in trivial:

places = gdf.to_wkb()                       # geometry column becomes bytes
result = con.execute("""
    select admin, count(*) as provinces from places group by 1
""").df()
print(result.to_string(index=False))

                   admin  provinces
United States of America         51

For geometry, convert first: gdf.to_wkb() returns a DataFrame whose geometry column holds WKB, which DuckDB reads with ST_GeomFromWKB.

5. Keep the split at the right place

The division that works:

DuckDB — reading files, filtering, joining, aggregating, anything over a million rows.
GeoPandas — the result, once it is small enough to hold comfortably: plotting, interactive maps, shapely operations on a few thousand features, export to formats GDAL writes well.

Measured, the heavy end justifies itself: 13,472,195 points joined to 4,596 polygons took 50.9 s and 291 MB in DuckDB against 121.8 s and 23,813 MB in GeoPandas (benchmark). At the light end the choice hardly matters: on 7,342 points the two took 0.26 s and 0.44 s.

6. Do not round-trip more than once

Each crossing costs a serialisation. A pipeline that alternates — query, convert, operate, convert back, query — spends most of its time in WKB.

Decide where the boundary is, do all the heavy work on one side, cross once.

shapely.from_wkb over a whole column is a single C-level loop.

Code examples

The data used on this page

The points are GeoNames' full allCountries dump — every feature, about 13.5 million, rebuilt daily at download.geonames.org, so counts drift a little between runs — kept both as the 1.79 GB text file and as places.parquet (also linked as geonames.parquet) with name, latitude, longitude, feature class, country and population. The polygons are Natural Earth's 1:10m admin-1 states and provinces, 4,596 of them, as provinces.shp.

import os
import zipfile

import duckdb
import requests


def fetch(url, path):
    if not os.path.exists(path):
        with requests.get(url, stream=True, timeout=600) as reply:
            reply.raise_for_status()
            with open(path, "wb") as out:
                for chunk in reply.iter_content(1 << 20):
                    out.write(chunk)


fetch("https://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip")   # about 420 MB
if not os.path.exists("allCountries.txt"):
    zipfile.ZipFile("allCountries.zip").extract("allCountries.txt")
if not os.path.exists("places.parquet"):
    duckdb.sql("""
        copy (select column01 as name, column04 as lat, column05 as lon, column06 as feature_class,
                     column08 as country, column14 as population
              from read_csv('allCountries.txt', delim = '\\t', header = false, quote = '',
                            all_varchar = false, sample_size = -1))
        to 'places.parquet' (format parquet)""")
if not os.path.exists("geonames.parquet"):
    os.link("places.parquet", "geonames.parquet")
fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip", "provinces.zip")
if not os.path.exists("provinces.shp"):
    with zipfile.ZipFile("provinces.zip") as z:
        for member in z.namelist():
            if member.startswith("ne_10m_admin_1_states_provinces."):
                open("provinces" + os.path.splitext(member)[1], "wb").write(z.read(member))
n = duckdb.sql("select count(*) from read_parquet('places.parquet')").fetchone()[0]
print(f"{n:,} GeoNames points; places.parquet {os.path.getsize('places.parquet') / 1e6:.0f} MB, "
      f"allCountries.txt {os.path.getsize('allCountries.txt') / 1e9:.2f} GB")

13,472,195 GeoNames points; places.parquet 306 MB, allCountries.txt 1.79 GB

Example 1 — the two conversions, as functions

import duckdb
import geopandas as gpd
import pandas as pd
from shapely import from_wkb


def duckdb_to_geodataframe(con, sql, crs, geometry_column="geom"):
    """Run a query and return a GeoDataFrame with the CRS you specify."""
    frame = con.execute(sql).df()
    if geometry_column not in frame.columns:
        raise KeyError(
            f"no column {geometry_column!r} in the result; "
            f"select st_aswkb({geometry_column}) as {geometry_column}")
    geometry = from_wkb(frame[geometry_column].map(bytes))   # bytearray → bytes
    frame = frame.drop(columns=[geometry_column])
    return gpd.GeoDataFrame(frame, geometry=geometry, crs=crs)


def geodataframe_to_duckdb(con, gdf, table_name, geometry_column="geometry"):
    """Register a GeoDataFrame as a DuckDB table with a real geometry column."""
    staged = gdf.to_wkb()                                  # geometry → bytes
    con.register("_staging", staged)
    con.execute(f"""
        create or replace table {table_name} as
        select * exclude {geometry_column},
               st_geomfromwkb({geometry_column}) as geom
        from _staging
    """)
    con.unregister("_staging")
    rows = con.execute(f"select count(*) from {table_name}").fetchone()[0]
    print(f"{table_name}: {rows:,} rows  (CRS {gdf.crs} — record it, "
          f"DuckDB will not keep it)")


geodataframe_to_duckdb(con, gdf, "us_areas")

us_areas: 51 rows  (CRS EPSG:4326 — record it, DuckDB will not keep it)

Example 2 — the heavy-then-light pattern

def heavy_then_light(points_parquet, polygons_path, crs="EPSG:4326", top=20):
    """DuckDB does the millions of rows; GeoPandas gets the summary."""
    con = duckdb.connect()
    con.execute("install spatial; load spatial;")
    con.execute("set enable_progress_bar = false")

    sql = f"""
        select a.name,
               count(*)               as points,
               st_aswkb(any_value(a.geom)) as geom
        from read_parquet('{points_parquet}') p
        join st_read('{polygons_path}') a
          on st_intersects(st_point(p.lon, p.lat), a.geom)
        group by a.name
        order by points desc
        limit {top}
    """
    gdf = duckdb_to_geodataframe(con, sql, crs=crs)
    print(f"{len(gdf)} areas crossed the boundary into GeoPandas")
    return gdf


import time

t = time.perf_counter()
busiest = heavy_then_light("places.parquet", "provinces.shp")
print(f"{time.perf_counter() - t:.1f} s; top three:")
print(busiest[["name", "points"]].head(3).to_string(index=False))

20 areas crossed the boundary into GeoPandas
52.5 s; top three:
         name  points
Uttar Pradesh  172436
   California  126373
        Texas  113915

any_value(a.geom) keeps one geometry per group without a second join. The whole point of the shape is that a 13-million-row input becomes a 20-row GeoDataFrame before anything is materialised in Python.

Example 3 — a round-trip test that catches silent loss

def assert_roundtrip(con, gdf, crs=None):
    """Geometry and CRS should survive a trip through DuckDB unchanged."""
    crs = crs or gdf.crs
    geodataframe_to_duckdb(con, gdf, "_roundtrip")
    back = duckdb_to_geodataframe(
        con, "select * exclude geom, st_aswkb(geom) as geom from _roundtrip", crs=crs)

    assert len(back) == len(gdf), f"row count changed: {len(gdf)} → {len(back)}"
    assert back.crs == gdf.crs, f"CRS changed: {gdf.crs} → {back.crs}"

    same = gdf.geometry.reset_index(drop=True).geom_equals(
        back.geometry.reset_index(drop=True))
    assert same.all(), f"{(~same).sum()} geometries differ after the round trip"
    print(f"round trip clean: {len(gdf):,} features, CRS {crs}")


assert_roundtrip(con, gdf)

_roundtrip: 51 rows  (CRS EPSG:4326 — record it, DuckDB will not keep it)
round trip clean: 51 features, CRS EPSG:4326

Run this once when you set the boundary up. WKB is lossless, so a failure means the CRS assumption is wrong or a geometry was invalid to begin with.

Example 4 — what a million rows cost to move

import io
import json

import numpy as np
import pyarrow as pa
import pyarrow.ipc as ipc

million = con.execute("select lon, lat from read_parquet('places.parquet') limit 1000000").df()
as_json = json.dumps(million.to_dict(orient="records")).encode()
buffer = io.BytesIO()
with ipc.new_file(buffer, pa.Table.from_pandas(million, preserve_index=False).schema,
                  options=ipc.IpcWriteOptions(compression="zstd")) as writer:
    writer.write_table(pa.Table.from_pandas(million, preserve_index=False))
raw32 = million.to_numpy(dtype=np.float32).tobytes()
import gzip

for label, size in (("JSON records", len(as_json)), ("JSON, gzipped", len(gzip.compress(as_json))),
                    ("Arrow IPC, zstd", buffer.tell()), ("raw float32", len(raw32))):
    print(f"{label:16} {size / 1e6:6.2f} MB")

for label, fetch_ in (("to pandas via .df()", lambda: con.execute("select * from read_parquet('places.parquet')").df()),
                      ("to Arrow via .arrow()", lambda: con.execute("select * from read_parquet('places.parquet')").arrow())):
    t = time.perf_counter()
    out = fetch_()
    if hasattr(out, "read_all"):
        out = out.read_all()
    print(f"{label:22} {time.perf_counter() - t:5.2f} s for {len(out):,} rows")

JSON records      36.25 MB
JSON, gzipped      6.82 MB
Arrow IPC, zstd   11.20 MB
raw float32        8.00 MB
to pandas via .df()     4.65 s for 13,472,195 rows
to Arrow via .arrow()   1.21 s for 13,472,195 rows

Explanation

Why WKB and not the native type

DuckDB's GEOMETRY is an internal representation optimised for its own operations. When it crosses into pandas it becomes an object pandas cannot interpret, and shapely cannot read it either.

WKB is the lingua franca: a documented binary encoding that both sides implement natively and vectorised. shapely.from_wkb on a whole column is a single C-level loop, which is why the conversion is fast enough not to think about.

Why the CRS has to be an explicit constant

DuckDB's CRS tracking is a type annotation, and annotations are lost by most operations — including the ST_AsWKB that gets the data out. GeoPandas needs a real CRS object to reproject, to join safely and to plot with a basemap.

The consequence is that the CRS must live in your code, next to the query, as a constant. That is not as bad as it sounds: a pipeline should know what CRS it works in, and writing it down is a documentation improvement as much as a technical necessity.

Why the direction of the split matters

Both engines can do everything. The reason to put the heavy work in DuckDB is memory: it streams, and GeoPandas materialises.

Measured on the global join, GeoPandas peaked at 23,813 MB against DuckDB's 291 MB — an 82× difference on identical output.

At small sizes the difference stops mattering: on 7,342 points against 4,596 polygons DuckDB took 0.26 s and GeoPandas 0.44 s, because query planning and structure building dominate both.

Why one crossing is the target

Serialising a million geometries to WKB and back is not free, and a pipeline that alternates between the two libraries does it repeatedly.

The design that avoids it: express the whole heavy stage as one query — joins, filters, aggregates and all — and cross once with a result small enough that the conversion cost is irrelevant. That is also the stage where SQL is at its most readable, because the whole transformation is in one place.

Use .arrow() rather than .df() above about a million rows.

Edge cases or notes

st_aswkb() in the select list, or the geometry arrives unusable.
crs= is not optional. A GeoDataFrame with crs=None will not reproject and will join wrongly.
gdf.to_wkb() returns a plain DataFrame with bytes in the geometry column, ready for con.register.
.arrow() avoids a copy and is worth it above about a million rows.
Empty and null geometries survive WKB but need handling on both sides.
Mixed geometry types are fine in DuckDB and awkward in some GeoPandas exports.
Register temporary frames and unregister them, or the connection holds references.
Do not put a WKB column in a Parquet file and call it GeoParquet — that needs the geo metadata key, which DuckDB's native writer does add.

Internal links

DuckDB spatial explained: a spatial database that is just a file — when to use which engine
How to run a spatial join in DuckDB — the heavy stage
CRS in DuckDB: why ST_Transform moves your data to the wrong place — why the CRS is not carried
How to benchmark DuckDB against GeoPandas honestly — deciding where the boundary is
WKT, WKB and GeoJSON explained — the interchange formats
What is a GeoDataFrame — the destination
Spatial SQL versus GeoPandas — the wider comparison
How to read shapefiles, GeoJSON and GeoParquet in DuckDB — getting data in

FAQ

How do I get a GeoDataFrame out of a DuckDB query?

Select st_aswkb(geom), read the result with .df(), and build the GeoDataFrame with shapely.from_wkb and an explicit crs=.

Why is my geometry column full of bytes?

Because WKB is bytes. Convert it with shapely.from_wkb, which is vectorised over the whole column.

Why does my GeoDataFrame have no CRS?

DuckDB tracks CRS on the column type and drops it in most operations, including ST_AsWKB. Set crs= explicitly from a constant in your code.

Should I use `.df()` or `.arrow()`?

.arrow() above about a million rows — it avoids an intermediate copy. Below that the difference is not worth the extra line.

How do I push a GeoDataFrame into DuckDB?

gdf.to_wkb(), register the resulting DataFrame, then st_geomfromwkb() in a create table as select.

Where should the boundary between the two be?

After the heavy stage. DuckDB streams a 13.5-million-point join in 291 MB where GeoPandas needs 23,813 MB (measured); cross once, with a result small enough that the conversion does not matter.