Fixing CRS Drift in GeoDataFrame Conversion

A GeoDataFrame reconstructed from a DuckDB result almost always arrives with crs=None, because DuckDB Spatial’s GEOMETRY type stores coordinates only and neither the WKB wire format nor the Arrow buffer that ferries it preserves a coordinate reference — the drift this guide, part of the DuckDB to GeoPandas sync workflow, teaches you to detect and repair before it silently poisons every to_crs, plot, and distance computed downstream.

Root-Cause Analysis of CRS Drift Across the Boundary #

CRS drift is not one bug; it is a chain of five metadata drops, each of which discards the coordinate reference at a different membrane between the source table and the finished GeoDataFrame. Because the coordinates themselves survive unchanged, every drop is invisible until something tries to interpret those coordinates — and by then the geometry is either naive (crs=None) or, worse, silently tagged with a CRS it does not actually use.

• The GEOMETRY type has no SRID slot. DuckDB Spatial’s internal geometry representation is coordinate-only. Unlike a PostGIS geometry(Point, 4326), a DuckDB value carries no embedded SRID, so the moment data lands in a DuckDB column the coordinate reference lives only in your head or your ingestion code, never in the column. The representation trade-offs behind this are covered in the ST_Geometry vs WKB reference.
• ST_AsWKB emits plain WKB, not EWKB. Standard WKB encodes geometry type and coordinates but has no field for an SRID. ST_AsWKB(geom) therefore produces a byte buffer that is, by construction, coordinate-reference-free. There is no SRID to lose because none was ever serialized.
• The Arrow BLOB column carries no GeoArrow CRS metadata. When the WKB column crosses into Python as an Arrow binary field, it travels as an opaque BLOB, not as a GeoArrow extension type with a crs key in its field metadata. The Arrow schema knows the column is bytes; it does not know those bytes are geometry, let alone which datum they sit on.
• shapely.from_wkb has no CRS concept at all. Shapely geometries are pure coordinate containers — they never held a CRS and never will. Parsing the buffer yields correct shapes with zero coordinate-reference information attached.
• GeoDataFrame defaults crs to None. With nothing upstream supplying a CRS, GeoPandas constructs the frame naive. The dangerous variant is the wrong-assumption drop: code that unconditionally passes crs="EPSG:4326" to a frame whose coordinates are actually in a metre-based projection. That geometry is no longer naive — it is confidently mislabelled, and no later call will error.

A sixth failure mode compounds the rest: axis-order drift. EPSG:4326 is formally defined as latitude-longitude, but virtually all data and every DuckDB coordinate constructor treat it as longitude-latitude. If a reprojection honours the authority axis order while your coordinates use the traditional one, the result is transposed by up to the width of the globe with no exception raised.

Deterministic Configuration #

The only durable fix is to make the coordinate reference an explicit, tracked value that travels alongside the query rather than inside the geometry. Pin the connection, load the extension, and define the source CRS as a named constant next to the SQL that produced the geometry.

import duckdb
import geopandas as gpd
import pandas as pd
import shapely
import pyarrow as pa

con = duckdb.connect(config={
    "threads": 4,                     # match physical cores; ST_Transform is CPU-bound
    "memory_limit": "8GB",            # hard ceiling — reprojection doubles the coord buffer briefly
    "preserve_insertion_order": False,# unlock parallel scans; re-order in SQL if needed
})
con.execute("INSTALL spatial; LOAD spatial;")

# The CRS is metadata the geometry cannot carry — track it out-of-band, beside the query.
SOURCE_CRS = "EPSG:32633"   # the CRS the stored coordinates are ACTUALLY in
TARGET_CRS = "EPSG:4326"    # what you want the GeoDataFrame tagged / reprojected to

Before trusting any converted frame, confirm these preconditions:

• The true CRS of the source geometry is known and recorded as a constant, not inferred from the column.
• The geometry column is emitted with ST_AsWKB, so the boundary carries binary, not stringified WKT.
• Every GeoDataFrame construction passes an explicit crs=; no frame is left to default to None.
• If a metric CRS is required for distance or area, the reprojection happens in SQL with ST_Transform before the handoff, or in GeoPandas with to_crs after the CRS label is set — never by silently relabelling.

Optimized Execution Pattern #

The broken pattern looks correct: it parses geometry vectorized and builds a frame. It simply never supplies a CRS, so the first reprojection throws.

# BEFORE — geometry is correct, but the frame is naive and every reprojection fails.
tbl = con.execute(
    "SELECT road_id, ST_AsWKB(geom) AS wkb FROM roads"
).fetch_arrow_table()

geoms = shapely.from_wkb(tbl.column("wkb").to_numpy(zero_copy_only=False))
attrs = tbl.drop_columns(["wkb"]).to_pandas(types_mapper=pd.ArrowDtype)
gdf = gpd.GeoDataFrame(attrs, geometry=geoms)   # no crs= → crs is None

print(gdf.crs)          # None
gdf.to_crs(epsg=3857)   # ValueError: Cannot transform naive geometries.
                        #   Please set a crs on the object first.

There are two correct repairs, and which one you use depends on whether you need the coordinates reprojected or merely labelled.

Repair A — label the frame with the CRS it already uses. When the stored coordinates are already in the CRS you want to work in, the fix is to attach the tag on construction. set_crs assigns or overrides the label without touching coordinates; to_crs (used next) is what actually transforms them.

# AFTER (A) — tag the geometry with the CRS it is genuinely in, then reproject safely.
gdf = gpd.GeoDataFrame(attrs, geometry=geoms, crs=SOURCE_CRS)  # explicit, no drift
# Equivalent post-hoc, and idempotent even if a wrong crs slipped in earlier:
gdf = gdf.set_crs(SOURCE_CRS, allow_override=True)

gdf_web = gdf.to_crs(epsg=3857)   # now transforms coordinates correctly

Repair B — reproject inside DuckDB before the handoff. When you need a metric CRS for the query itself (buffers, ST_Distance, grid binning), do the transform in SQL. DuckDB’s vectorized ST_Transform is faster than a per-geometry GeoPandas to_crs and keeps the whole reduction in the engine. Set always_xy so the traditional longitude-latitude order is honoured and the axis-order drift above cannot occur.

# AFTER (B) — reproject in SQL; the frame is then tagged with the TARGET CRS.
tbl = con.execute("""
    SELECT
        road_id,
        ST_AsWKB(
            ST_Transform(geom, ?, ?, always_xy := true)  -- source -> target, lon/lat order
        ) AS wkb
    FROM roads
""", [SOURCE_CRS, TARGET_CRS]).fetch_arrow_table()

geoms = shapely.from_wkb(tbl.column("wkb").to_numpy(zero_copy_only=False))
attrs = tbl.drop_columns(["wkb"]).to_pandas(types_mapper=pd.ArrowDtype)
gdf = gpd.GeoDataFrame(attrs, geometry=geoms, crs=TARGET_CRS)  # tag = the CRS SQL produced

The behavioural change in both repairs is the same one idea stated twice: the CRS is supplied from outside the geometry, because the geometry cannot supply it. The distinction that trips people is set_crs versus to_crs — set_crs changes the label and leaves coordinates alone, to_crs moves the coordinates to match a new label. Relabelling projected data as EPSG:4326 with set_crs does not fix drift; it manufactures it. The unit and datum rules that decide which transform you need are laid out in CRS mapping and transformations and, at coordinate-system level, in how DuckDB Spatial handles coordinate systems.

For the mechanics of parsing the WKB column without an extra copy — and where the CRS tag slots into a GeoArrow-aware handoff — see the zero-copy Arrow to GeoPandas handoff reference, and for the full conversion function this pattern extends, converting DuckDB queries to GeoDataFrame efficiently.

Diagnostic Queries & Plan Validation #

Drift is cheap to detect if you assert on it rather than discover it in a broken map. The first guard is a hard check that the frame is never naive.

def assert_crs(gdf, expected_epsg: int) -> None:
    """Fail loudly on drift: naive frame, or a CRS that is not the one intended."""
    if gdf.crs is None:
        raise ValueError("CRS drift: GeoDataFrame is naive — coordinates cannot be interpreted.")
    got = gdf.crs.to_epsg()
    if got != expected_epsg:
        raise ValueError(f"CRS drift: frame is EPSG:{got}, expected EPSG:{expected_epsg}.")

The subtler case is the mislabelled frame, which passes the check above because it has a CRS — just the wrong one. Coordinate ranges betray it: projected metre coordinates run to hundreds of thousands, while geographic degrees stay within [-180, 180] and [-90, 90]. A cheap range probe in SQL catches a frame that claims to be EPSG:4326 but holds projected values.

-- If a frame is tagged 4326, X/Y must fall in degree ranges. Values in the
-- tens of thousands mean the geometry is projected and the tag is a lie.
SELECT
    min(ST_X(ST_Centroid(geom))) AS min_x,
    max(ST_X(ST_Centroid(geom))) AS max_x,
    max(abs(ST_Y(ST_Centroid(geom)))) AS max_abs_y   -- > 90 ⇒ not geographic
FROM roads;

When you push the reprojection into SQL (Repair B), validate that ST_Transform did not degrade the scan into per-row evaluation. A healthy plan shows the transform inside a vectorized PROJECTION above a plain TABLE_SCAN, not a row-at-a-time fallback.

EXPLAIN
SELECT road_id, ST_AsWKB(ST_Transform(geom, 'EPSG:32633', 'EPSG:4326', always_xy := true))
FROM roads;

Assert two things in the output: the ST_Transform call appears inside a PROJECTION (not a SCALAR_FUNCTION_EVAL in a ROW_EXECUTION node), and the row estimate at the scan matches the table cardinality — a large gap points to stale statistics that will make the transform’s cost unpredictable across data versions. The same plan-fingerprinting discipline used across the DuckDB to GeoPandas sync guide applies unchanged here.

Geometry Validation & Fallback Routing #

Reprojection is where invalid or out-of-range geometry finally errors, because ST_Transform evaluates coordinates against the projection’s area of use. A longitude outside a projected CRS’s valid domain transforms to infinity or NaN, which then silently propagates into every distance and bounding box. Guard the input before transforming and quarantine anything that cannot be projected cleanly.

-- Isolate geometry that will not survive the transform: invalid topology, or
-- coordinates outside the target projection's usable domain.
SELECT road_id, geom
FROM roads
WHERE NOT ST_IsValid(geom)
   OR ST_X(ST_Centroid(geom)) NOT BETWEEN -180 AND 180;

Repair topology with ST_MakeValid before the transform runs, so a self-intersection cannot throw mid-projection:

UPDATE roads SET geom = ST_MakeValid(geom) WHERE NOT ST_IsValid(geom);

After reprojection, check that no coordinate became non-finite. GeoPandas will happily hold inf in a geometry and only fail much later, so route non-finite rows out at the boundary:

import numpy as np

xs = shapely.get_coordinates(gdf.geometry)          # (n, 2) numpy array
finite = np.isfinite(xs).all(axis=1)
if not finite.all():
    bad = gdf.loc[~gdf.geometry.map(lambda g: np.isfinite(shapely.get_coordinates(g)).all())]
    gdf = gdf.drop(bad.index)   # quarantine, then reconcile separately

When the source table is large enough that a single reprojected result would breach memory_limit, transform and hand off in bounded chunks rather than materialising the whole frame — the transform doubles the coordinate buffer for the duration of the projection, so a result that fits as input can spill as output. Stream record batches, tag each chunk with the target CRS on construction, and concatenate:

result = con.execute("""
    SELECT road_id, ST_AsWKB(ST_Transform(geom, ?, ?, always_xy := true)) AS wkb
    FROM roads
""", [SOURCE_CRS, TARGET_CRS])

frames = []
for batch in result.fetch_record_batch(500_000):
    if batch.num_rows == 0:
        continue
    t = pa.Table.from_batches([batch])
    g = shapely.from_wkb(t.column("wkb").to_numpy(zero_copy_only=False))
    a = t.drop_columns(["wkb"]).to_pandas(types_mapper=pd.ArrowDtype)
    frames.append(gpd.GeoDataFrame(a, geometry=g, crs=TARGET_CRS))  # tag every chunk

gdf = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), crs=TARGET_CRS)

Tagging each chunk with crs=TARGET_CRS and re-asserting it on the concatenated frame closes the loop: no chunk can drift, and the final pd.concat cannot silently reset the CRS to None.

Related #

See also

• Converting DuckDB queries to GeoDataFrame efficiently — the full conversion function this CRS handling plugs into.
• Zero-copy Arrow to GeoPandas handoff — where a GeoArrow-aware boundary can preserve the CRS tag instead of dropping it.
• CRS mapping and transformations — the datum and unit rules that decide whether you relabel or reproject.

Up: DuckDB to GeoPandas Sync · Python & DuckDB Integration Workflows