DuckDB Spatial Analytical SQL for GIS

Converting Shapely Geometries to DuckDB Safely

Pushing geometry from Python into DuckDB is the direction with fewer articles written about it and more ways to go quietly wrong, because the conversion happens per object unless you stop it and the frame is lost unless you carry it. This walkthrough, part of the Shapely integration reference, covers the bulk conversion, the four things that do not travel, and the validation that belongs on the Python side rather than after the load.

Root-Cause Analysis: what goes wrong pushing geometry in

  • The conversion happens per object. Passing a column of Shapely objects means each one is serialised individually as the frame is read. Converting the whole array once with to_wkb is a single compiled pass and is the difference between a bulk copy and a Python loop.
  • The reference frame does not travel. Shapely objects carry no CRS at all, and a GeoDataFrame carries one that DuckDB has nowhere to store. Whatever frame the geometry is in has to be recorded separately or it is gone at the moment of transfer.
  • Invalid geometry arrives intact. to_wkb serialises whatever it is given, including self-intersecting rings. The failure then appears inside a later overlay, several stages from where it entered.
  • Empty and None are conflated. A column mixing None and empty geometries produces a mixture of SQL NULL and valid-but-empty values, which behave differently in every predicate and are easy to treat as one state.
  • The declared type widens with size. A WKB column that exceeds the 32-bit offset limit becomes a wider Arrow type, so the same pipeline sends a different type on a larger day.

The distinguishing question is whether the geometry crosses as an array of bytes or as a sequence of objects. Everything else on this page follows from that one decision.

Four ways in, by cost A Shapely column is ~120x, row-by-row insert ~300x, a bulk to_wkb column ~9x, and a declared Arrow table is the baseline. HOW IT CROSSES RELATIVE COST WHY a GeoDataFrame with Shapely geometry ~120× one serialisation per object row-by-row INSERT ~300× a statement per row as well a WKB column via to_wkb ~9× one compiled pass, then a bulk scan an Arrow table, declared schema 1× already in the engine’s layout

The bottom two rows differ by a factor. The top two differ from them by two orders of magnitude.

Deterministic Configuration 

import duckdb, pyarrow as pa
from shapely import to_wkb, from_wkb, is_valid, make_valid

con = duckdb.connect("gis.duckdb")
con.execute("INSTALL spatial; LOAD spatial")
con.execute("SET memory_limit = '6GB'")

# Declared once, used at every boundary — inference here produces a type
# that depends on how much data happened to move.
IN_SCHEMA = pa.schema([
    pa.field("parcel_id", pa.int64(),        nullable=False),
    pa.field("wkb",       pa.large_binary(), nullable=True),
])
SOURCE_CRS = "EPSG:27700"        # recorded here because nothing carries it

Optimized Execution Pattern

The pattern is one vectorised serialisation, one registration against a declared schema, and one SQL statement that converts and validates as it loads.

# ANTI-PATTERN: the geometry column crosses as objects, so every row is
# serialised individually — and the CRS is lost without anyone noticing.
con.register("staged", gdf)                      # Shapely objects, one at a time
con.execute("CREATE TABLE parcels AS SELECT * FROM staged")
# PATTERN: serialise once, declare the schema, and convert inside SQL.
wkb = to_wkb(gdf.geometry.values)                # one compiled pass over the array

tbl = pa.table({"parcel_id": gdf["parcel_id"].to_numpy(), "wkb": wkb}, schema=IN_SCHEMA)
con.register("staged", tbl)

con.execute("""
    CREATE OR REPLACE TABLE parcels AS
    SELECT parcel_id,
           ST_GeomFromWKB(wkb) AS geom,
           CASE WHEN wkb IS NULL THEN 'null' ELSE 'ok' END AS geom_state
    FROM staged
""")
# and record the frame, because the table does not carry it
con.execute("INSERT INTO layer_crs VALUES ('parcels', ?)", [SOURCE_CRS])

The layer_crs insert is two lines and it is the only thing standing between a correct pipeline and one that computes areas in the wrong units six months from now. Recording the frame at the moment of transfer is the only point at which the information is definitely available.

Four things that do not survive the transfer The CRS, the frame index, any validity guarantee, and the distinction between None and empty. DOES NOT SURVIVE WHY WHAT TO DO the coordinate frame Shapely never had it; DuckDB cannot store it record it as metadata the frame index SQL has no positional index an explicit key column any validity guarantee serialisation checks nothing validate before or gate after None versus empty easy to flatten in conversion preserve it deliberately

Only the second is inconvenient. The other three are silent.

Validating on the Python side

Validity is worth checking before the transfer rather than after it, for a reason that is about diagnosis rather than performance: on the Python side you still have the object, the row, and the context that produced it, so a failure names something actionable. After the load you have a row identifier and a reason string.

The check itself is vectorised in modern Shapely, so it costs a single compiled pass over the array rather than a loop. Repairing is likewise vectorised, and the three-way route — pass, repair, quarantine — is the same shape as the SQL-side gate, just applied where the context still exists.

# Vectorised, one pass, with the three outcomes separated before transfer.
geoms = gdf.geometry.values
ok    = is_valid(geoms)

repaired = make_valid(geoms[~ok])
kept     = is_valid(repaired) & (
    [g.geom_type for g in repaired] == [g.geom_type for g in geoms[~ok]]
)

print(f"clean={ok.sum()} repaired={kept.sum()} quarantined={(~ok).sum() - kept.sum()}")

Diagnostic Queries & Plan Validation

After the load, two counts confirm that what arrived is what was sent.

-- Row count and geometry-state census against what Python reported sending.
SELECT count(*)                                       AS rows_total,
       count(*) FILTER (WHERE geom IS NULL)           AS rows_null,
       count(*) FILTER (WHERE ST_IsEmpty(geom))       AS rows_empty,
       count(*) FILTER (WHERE NOT ST_IsValid(geom))   AS rows_invalid
FROM parcels;

A non-zero invalid count after a Python-side validation means something was repaired into a shape that still fails GEOS in the engine, which is rare and worth investigating rather than repairing again. A non-zero empty count usually means None and empty were conflated on the way in.

Three post-transfer symptoms A slow load means per-object conversion; areas wrong by a factor mean a lost frame; vanishing rows mean None and empty were conflated. SYMPTOM POINTS AT FIX the load is far slower than the volume geometry crossing as objects one vectorised to_wkb areas wrong by a consistent factor the frame was not recorded record and reproject rows present but absent from areal queries None and empty conflated preserve the distinction

Three symptoms, three of the four things that do not survive the transfer.

Streaming a frame too large to convert at once

to_wkb over a whole column allocates the whole serialised result, so a frame close to the memory ceiling can fail during the conversion rather than during the analysis. The fix is to convert in batches and register a record-batch reader rather than a table, which keeps the peak at one batch.

That also composes with the validation above: each batch can be validated, repaired and counted before it is handed over, so a large transfer reports its three-way census incrementally rather than all at the end. It is more code than the single-shot version and it is the version that survives a frame growing by an order of magnitude.

# Batch the conversion so the peak is one batch rather than the whole column.
def batches(gdf, size=100_000):
    for start in range(0, len(gdf), size):
        chunk = gdf.iloc[start:start + size]
        yield pa.record_batch(
            {"parcel_id": chunk["parcel_id"].to_numpy(),
             "wkb": to_wkb(chunk.geometry.values)}, schema=IN_SCHEMA)

reader = pa.RecordBatchReader.from_batches(IN_SCHEMA, batches(gdf))
con.register("staged", reader)
con.execute("CREATE OR REPLACE TABLE parcels AS SELECT parcel_id, ST_GeomFromWKB(wkb) AS geom FROM staged")

Geometry Validation & Fallback Routing

Where a geometry cannot be serialised at all, fail with the row rather than with the column.

# Serialise defensively so one unserialisable object names itself instead of
# failing the whole array with a message about the array.
out = []
for idx, g in zip(gdf.index, gdf.geometry.values):
    try:
        out.append(to_wkb(g))
    except Exception as exc:
        raise ValueError(f"row {idx}: {type(g).__name__} could not be serialised: {exc}") from exc

Frequently Asked Questions

Why is pushing a GeoDataFrame into DuckDB slow?

Because the geometry column crosses as Shapely objects, so each one is serialised individually as the frame is scanned. Converting the whole array once with to_wkb and passing a plain binary column is a single compiled pass, and the difference is more than an order of magnitude.

What happens to the CRS?

Nothing carries it. Shapely objects have never had one, and DuckDB has nowhere to store it, so the frame exists only in the GeoDataFrame you are reading from. Record it explicitly at the moment of transfer — that is the only point at which the information is definitely available.

Should I validate before or after the transfer?

Before, where possible. On the Python side you still have the object and the context that produced it, so a failure is actionable; after the load you have an identifier and a reason string. The check is vectorised, so it costs one compiled pass rather than a loop.

How do I keep None and empty distinct?

Decide what each means and encode it. None becomes SQL NULL; an empty geometry becomes a valid value that matches nothing. Flattening them is easy during conversion and produces rows that are present in the table and absent from every areal query, with nothing to explain why.

Is WKB or WKT the right format for this direction?

WKB, without exception outside debugging. WKT is several times larger, has to be parsed character by character on arrival, and loses precision at the last digits unless the formatting is handled carefully. WKB is the representation both sides already speak.

What if the frame is too large to convert in one call?

Batch it. to_wkb over a whole column allocates the entire serialised result, so a large frame can fail during conversion rather than during analysis. Converting in batches and registering a record-batch reader keeps the peak at one batch and lets the validation report incrementally.

Up: Shapely Integration with DuckDB