CRS Mapping & Transformations in DuckDB Spatial

Coordinate Reference System (CRS) handling is the most common silent-failure surface in analytical GIS pipelines: a transform that succeeds syntactically can still shift geometries by metres — or by entire hemispheres — when datum parameters, axis order, or units are mismatched. This guide covers the specific workflow of reprojecting mixed-CRS data to a single working CRS inside DuckDB Spatial Architecture & Fundamentals: how ST_Transform routes through PROJ, how to keep the operation in the projection phase of the plan, and how to detect drift before it contaminates downstream joins. It assumes you already understand the engine’s columnar layout and want deterministic, reproducible reprojection at scale.

DuckDB decouples coordinate storage from projection metadata — geometries are stored as raw WKB with no inline SRID. The CRS for a column must therefore be tracked out-of-band (a companion srid column, a side table, or documented schema contract), and every transform must name both endpoints explicitly. That single architectural fact drives every pattern below; the lower-level mechanics of how the engine resolves a projection are detailed in how DuckDB Spatial handles coordinate systems.

Runtime Configuration & Memory Guardrails #

Reprojection is a per-row vertex operation that materialises new geometry buffers, so the working set can briefly exceed the size of the source column. Configure hard boundaries before any bulk ST_Transform, and pre-warm the PROJ context cache so the first batch does not pay initialisation latency mid-scan.

INSTALL spatial;
LOAD spatial;

-- Cap resident memory BELOW total RAM: reprojection allocates fresh geometry
-- buffers, so leave headroom for the source column to stay cached.
SET memory_limit = '8GB';

-- More threads parallelise the vectorized transform kernel until memory
-- bandwidth saturates; past that point extra threads only raise peak memory.
SET threads = 8;

-- CRS output ordering is irrelevant for a reprojection staging step, so drop
-- order preservation to remove a buffering penalty on large materializations.
SET preserve_insertion_order = false;

-- Fast local NVMe for spill: a slow temp dir turns OOM-avoidance into I/O thrash
-- when a wide geometry expansion overflows memory_limit.
SET temp_directory = '/var/lib/duckdb/spill';

-- Pre-warm the PROJ context for each CRS pair you will use, so the first bulk
-- batch does not stall reinitializing the projection pipeline mid-scan.
SELECT ST_Transform(ST_Point(0, 0), 'EPSG:25832', 'EPSG:4326');

PROJ context initialisation carries a per-unique-pair latency cost. DuckDB Spatial caches these contexts in a thread-local pool, so a stable set of CRS pairs amortises cleanly. High-cardinality CRS distributions — say more than ~50 distinct input projections in one scan — will thrash that cache. Mitigate by normalising input CRSs upstream, or by partitioning the scan so each batch carries a single source CRS and the cache stays warm.

Primary Execution Patterns #

The canonical reprojection pipeline is a two-stage flow: validate and tag the source CRS, then transform and materialise into one target CRS before any spatial join runs. Materialising once is critical — DuckDB stores no inline SRID, so a transform left inside a join’s predicate re-runs PROJ for every candidate pair instead of once per row.

Normalize every layer to a single CRS before joining — DuckDB stores no inline SRID, so the transform must be explicit and must happen once, upstream of any join.

Stage one tags each row with a validated source CRS. ST_Transform requires VARCHAR CRS strings, so the EPSG code is concatenated at call time:

-- Stage 1: materialize source with explicit CRS tracking. The srid column
-- records the known source CRS per row; COALESCE supplies a documented default.
CREATE OR REPLACE TABLE staging_raw AS
SELECT
    id,
    COALESCE(srid, 4326) AS validated_srid,
    geom
FROM read_parquet('s3://bucket/raw/*.parquet');

-- Stage 2: deterministic transform to a single working CRS (EPSG:4326).
-- Guarding on ST_IsValid keeps malformed input out of the PROJ pipeline.
CREATE OR REPLACE TABLE staging_transformed AS
SELECT
    id,
    ST_Transform(geom, 'EPSG:' || validated_srid::VARCHAR, 'EPSG:4326') AS geom_4326,
    validated_srid AS source_srid
FROM staging_raw
WHERE ST_IsValid(geom);

The source format you read from materially changes this stage’s cost. Reprojecting straight off GeoParquet ingestion lets the reader carry embedded CRS metadata and prune columns before the transform ever sees a vertex, whereas GeoJSON ingestion forces a row-by-row parse and a larger transient footprint. Whether the staging table stays resident or spills is governed by the same trade-offs covered in in-memory vs disk storage; memory-resident staging tables give the transform kernel cache locality, while disk-backed scans need the memory_limit headroom configured above.

For programmatic pipelines, wrap both stages in an explicit transaction so a transform failure cannot leave a half-written staging table behind:

import duckdb
import logging

def execute_crs_pipeline(conn: duckdb.DuckDBPyConnection, source_path: str) -> None:
    try:
        conn.execute("BEGIN TRANSACTION;")
        conn.execute(f"""
            CREATE OR REPLACE TABLE staging_raw AS
            SELECT id, COALESCE(srid, 4326) AS validated_srid, geom
            FROM read_parquet('{source_path}');
        """)
        conn.execute("""
            CREATE OR REPLACE TABLE staging_transformed AS
            SELECT
                id,
                ST_Transform(geom, 'EPSG:' || validated_srid::VARCHAR, 'EPSG:4326') AS geom_4326,
                validated_srid AS source_srid
            FROM staging_raw
            WHERE ST_IsValid(geom);
        """)
        conn.execute("COMMIT;")
    except Exception as e:
        conn.execute("ROLLBACK;")
        logging.error(f"CRS pipeline failed: {e}")
        raise

Once staging_transformed holds a single CRS, downstream proximity work — for example spatial joins and proximity filters — operates on consistent coordinates and the planner can route bbox pre-filters through an R-tree without per-pair reprojection.

Execution Plan Validation #

The defining anti-pattern in CRS work is a transform that drifts out of the projection phase and into a join or filter, where it re-evaluates far more often than once per row. Confirm placement with EXPLAIN before promoting any pipeline:

EXPLAIN
SELECT id, ST_Transform(geom, 'EPSG:' || validated_srid::VARCHAR, 'EPSG:4326') AS geom_4326
FROM staging_raw
WHERE ST_IsValid(geom);

ST_Transform should resolve inside a PROJECTION node sitting directly above the table scan. If it instead appears inside a FILTER, or worse inside a NESTED_LOOP_JOIN, the planner is reprojecting candidate rows repeatedly — a CRS pipeline that should be $O(N)$ collapses toward $O(N \times M)$. The usual cause is a transform written into a join’s ON/WHERE predicate instead of a pre-materialised column.

Use EXPLAIN ANALYZE to attach real timings and row counts, and watch two thresholds:

EXPLAIN ANALYZE
SELECT id, ST_Transform(geom, 'EPSG:' || validated_srid::VARCHAR, 'EPSG:4326') AS geom_4326
FROM staging_raw
WHERE ST_IsValid(geom);

• Row-estimate drift. If the estimated cardinality on the PROJECTION node diverges from actual by more than ~2x, the optimizer is mis-sizing intermediate buffers and may over-allocate during the transform. Re-ANALYZE the source table to refresh statistics.
• Transform share of wall time. On a healthy plan the scan dominates and ST_Transform is a thin slice of total time. If the transform node alone exceeds ~40% of wall time, you are almost certainly thrashing the PROJ context cache (too many distinct CRS pairs in one scan) rather than bottlenecking on vertices.

Performance Trade-offs #

Reprojection cost is dominated by two factors: how the data is read, and how many distinct CRS pairs the cache must juggle. The ingestion format sets the floor on both parse latency and the transient memory the transform must fit inside:

Three quantified rules follow from production pipelines:

• Single-CRS batches beat mixed batches. Partitioning a scan so each batch carries one source CRS keeps the PROJ cache warm and typically removes 15–30% of transform overhead versus interleaving many projections in one pass.
• Geographic→projected transforms are cheap; datum-shifting transforms are not. A pure axis-and-scale transform (e.g. EPSG:4326 → EPSG:3857) is a closed-form calculation. A datum shift that needs a grid (NTv2/TOWGS84) can be several times slower per vertex — reserve those for data that genuinely needs the accuracy.
• Materialise once, reuse many. Reprojecting into a staging table and joining against it is strictly cheaper than reprojecting inside each query, because the transform runs once per row instead of once per candidate pair. The break-even is immediate for any geometry touched by more than one downstream query.

Apply the lightest variant that preserves the accuracy your analysis requires: skip the transform entirely when source and target already match, prefer projected CRSs for distance/area work so results come back in metres rather than degrees, and only invoke datum grids when sub-metre fidelity is contractual.

Edge Cases & Anti-Patterns #

Silent drift occurs when a transform returns geometry that is syntactically valid but geometrically wrong. Establish a validation boundary at each stage rather than trusting the absence of an exception.

Axis-order inversion (lon/lat vs lat/lon). The single most frequent CRS bug. EPSG:4326 is formally defined as latitude-first, but modern PROJ (and DuckDB Spatial’s default) treats it as longitude-first. Data exported from a lat-first toolchain lands with coordinates swapped, and the geometry appears valid while sitting in the wrong hemisphere. Catch it with a range gate on the output:

-- Reject rows whose transformed coordinates fall outside valid geographic range.
SELECT id,
       ST_XMin(geom_4326) AS min_lon,
       ST_YMin(geom_4326) AS min_lat
FROM staging_transformed
WHERE NOT (ST_XMin(geom_4326) BETWEEN -180 AND 180
       AND ST_YMin(geom_4326) BETWEEN  -90 AND  90);

Ambiguous or deprecated EPSG codes. Transform against an allow-list instead of accepting whatever the source claims; an unknown code is a data-quality signal, not a transform input:

SELECT id, geom
FROM staging_raw
WHERE validated_srid NOT IN (4326, 3857, 25832, 32601, 32632)  -- allowed codes
   OR ST_IsEmpty(geom)
   OR NOT ST_IsValid(geom);

Round-trip precision loss. Transforming a high-precision projected CRS (e.g. a UTM zone) to a global CRS and back can introduce sub-metre rounding. Enforce a tolerance with a round-trip check and quarantine anything that moves too far:

-- Flag geometries that shift more than 0.5 units after a round trip.
SELECT id
FROM staging_transformed
WHERE ST_Distance(
    ST_Transform(geom_4326, 'EPSG:4326', 'EPSG:' || source_srid::VARCHAR),
    (SELECT geom FROM staging_raw sr WHERE sr.id = staging_transformed.id)
) > 0.5;

Invalid geometry into PROJ. Self-intersecting or empty geometries can make ST_Transform return NULL or raise PROJ: invalid projection. Guard at the boundary and route repairable inputs through ST_MakeValid before retrying, rather than dropping them blindly:

-- Repair-then-transform fallback for inputs PROJ rejected.
CREATE OR REPLACE TABLE staging_repaired AS
SELECT id, validated_srid, ST_MakeValid(geom) AS geom
FROM staging_raw
WHERE NOT ST_IsValid(geom) AND geom IS NOT NULL;

Assuming a GRANT-style access model. DuckDB has no roles or GRANT/REVOKE. Exposing reprojected data safely means publishing a curated view and sharing the underlying database file read-only, not revoking column privileges:

-- Expose a single-CRS view; distribute it via a READ_ONLY attachment and keep
-- the writable staging table in a separate, restricted database file.
CREATE SCHEMA IF NOT EXISTS analytics;
CREATE OR REPLACE VIEW analytics.crs_4326 AS
SELECT id, geom_4326
FROM staging_transformed
WHERE source_srid IN (4326, 3857, 32632);

ATTACH 'staging.duckdb' AS staging (READ_ONLY);

A further deployment anti-pattern is plaintext cloud credentials embedded in read_parquet URIs; use CREATE SECRET so the connection string never carries the key. Centralising all reprojection on one node simplifies CRS standardisation but makes that node a single point of failure under concurrent load — keep ingestion and ST_Transform materialisation on the primary and serve analytical reads from read-only replicas.

Query Regression Analysis #

CRS pipelines regress quietly: a library upgrade changes a default axis interpretation, or a new source introduces a CRS pair that thrashes the cache, and the plan silently moves the transform into a join. Capture the plan as JSON, snapshot a baseline, and diff on every run so a placement or cost regression fails the pipeline instead of shipping drift.

import duckdb
import json

def capture_transform_plan(conn: duckdb.DuckDBPyConnection, sql: str) -> dict:
    """Return the JSON execution plan for a reprojection query."""
    conn.execute("PRAGMA enable_profiling = 'json';")
    plan = conn.execute("EXPLAIN (FORMAT JSON) " + sql).fetchone()[0]
    return json.loads(plan)

def assert_transform_in_projection(plan: dict) -> None:
    """Fail if ST_Transform has drifted out of a PROJECTION node into a join."""
    def walk(node):
        name = node.get("name", "")
        body = json.dumps(node.get("extra_info", node))
        if "ST_Transform" in body and ("JOIN" in name.upper() or "FILTER" in name.upper()):
            raise AssertionError(f"ST_Transform leaked into {name}; expected PROJECTION")
        for child in node.get("children", []):
            walk(child)
    walk(plan)

def diff_against_baseline(plan: dict, baseline_path: str) -> None:
    """Compare node names + estimated cardinalities against a stored baseline."""
    with open(baseline_path) as fh:
        baseline = json.load(fh)
    def shape(node):
        return {
            "name": node.get("name"),
            "card": node.get("cardinality"),
            "children": [shape(c) for c in node.get("children", [])],
        }
    if shape(plan) != shape(baseline):
        raise AssertionError("Plan shape drifted from baseline — review before deploy")

Run assert_transform_in_projection in CI against the staging query, and refresh the stored baseline deliberately (never automatically) whenever an intended optimisation changes the plan. Pair this with the round-trip tolerance query above as a data-level regression gate, so both the plan and the numeric output are guarded. When the same reprojected output feeds Python tooling, validate the CRS survives the handoff in DuckDB-to-GeoPandas sync, where a missing crs assignment on the GeoDataFrame is a classic post-transform drift point.

Related #

See also

• How DuckDB Spatial handles coordinate systems — the projection-resolution mechanics underneath ST_Transform.
• GeoParquet parsing — reprojecting straight off embedded CRS metadata with column pruning.
• In-memory vs disk storage — whether your staging table stays resident or spills.
• Spatial joins and proximity filters — what runs after every layer shares one CRS.
• DuckDB-to-GeoPandas sync — preserving CRS across the Python handoff.

Up: DuckDB Spatial Architecture & Fundamentals

---

External Reference Standards #

For authoritative CRS definitions, transformation parameters, and format specifications, consult the PROJ documentation, the EPSG Geodetic Parameter Dataset, and the GeoParquet specification.