How DuckDB Spatial Handles Coordinate Systems

DuckDB Spatial stores no inline SRID, so a GEOMETRY column that looks correct can still be silently misaligned with every layer you join it against. This page is the engine-level companion to CRS mapping and transformations within DuckDB Spatial Architecture & Fundamentals: it explains exactly how the engine resolves coordinate frames during ingestion, transformation, and indexing, and how to stop coordinate drift before it contaminates downstream joins.

Unlike PostGIS, which embeds a 4-byte SRID in every geometry, DuckDB’s GEOMETRY type stores raw WKB coordinates with no inline projection tag. The coordinate reference system (CRS) for a column is therefore an application-level concern: it must live in a companion srid column, a metadata side-table, or a documented schema contract. When no CRS is tracked, every spatial operator assumes whatever coordinate frame the bytes happen to be in — which is why a GEOMETRY column that round-trips cleanly through ST_AsText can still produce zero join matches against a layer that is one datum, one axis order, or one unit away.

Root-Cause Analysis: Where Coordinate Handling Silently Breaks #

Coordinate-system bugs in DuckDB almost never raise an error. They surface as empty result sets, geometries shifted by metres or hemispheres, or bounding boxes that pruning logic discards. The failure modes are distinct and each has a different fix:

The last two are the most expensive in production. The R-tree index internals construct the tree over raw bounding-box coordinates in whatever CRS the column already holds — the index has no concept of projection distortion, geodesic curvature, or unit normalisation. A predicate such as ST_DWithin or ST_Intersects therefore operates in untransformed coordinate space, so a cross-CRS join produces false negatives or catastrophic bounding-box mismatches. The rule that follows is non-negotiable: normalise every layer to one working CRS before you build an index and before you run a spatial join.

Ingestion is where the metadata is lost. GeoJSON follows a fixed convention — WGS 84 (EPSG:4326) with [longitude, latitude] order — so GeoJSON ingestion via st_read() or st_geomfromgeojson() parses coordinates directly, but no CRS tag survives into the column because DuckDB stores none. GeoParquet ingestion is richer: the extension reads the geo block in the Parquet footer (a PROJJSON object or an EPSG code) and uses it to inform ST_Transform, but it still does not embed the SRID into the geometry values. If the footer is malformed the parser can fall back silently, so confirm the embedded CRS before you trust it:

-- Inspect the GeoParquet footer to confirm the embedded CRS BEFORE ingestion.
-- A missing or malformed 'geo' key means you must supply the source CRS yourself.
SELECT key, value
FROM parquet_kv_metadata('s3://bucket/parcels.parquet')
WHERE key = 'geo';

Deterministic Configuration #

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

INSTALL spatial;
LOAD spatial;

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

-- Spill directory on fast local NVMe: large transforms overflow here rather
-- than triggering an OOM kill. Avoid network mounts (latency dominates).
SET temp_directory = '/mnt/nvme/duckdb_spill';

-- Match physical cores: the vectorized transform kernel parallelises until
-- memory bandwidth saturates; extra threads past that only raise peak memory.
SET threads = 4;

PROJ datum-shift grids (for example us_noaa_conus.tif) are cached in ~/.duckdb/proj on first use. In air-gapped environments, pre-stage the proj-data package and point PROJ_DATA at it, otherwise a transform that needs a grid will silently degrade to a less accurate sphere-based shift — the classic source of the consistent few-metre offset above.

Optimized Execution Pattern #

The naive pattern assumes a CRS and joins directly. The corrected pattern names both endpoints explicitly, normalises to one working frame, and only then builds the index or join.

-- BEFORE: assumes both layers share a CRS. Silently wrong if they don't —
-- the R-tree is built over coordinates in an unknown, possibly mixed frame.
CREATE INDEX idx_parcels ON parcels USING RTREE (geom);

SELECT p.id, s.id
FROM parcels p
JOIN sites s ON ST_DWithin(p.geom, s.geom, 100);   -- 100 *what*? degrees? metres?

-- AFTER: normalise every layer to one metric working CRS, validate, then index.
-- ST_Transform names source AND target explicitly because no SRID is stored.
CREATE OR REPLACE TABLE parcels_3857 AS
SELECT id, ST_Transform(geom, 'EPSG:4326', 'EPSG:3857') AS geom
FROM parcels;

CREATE OR REPLACE TABLE sites_3857 AS
SELECT id, ST_Transform(geom, 'EPSG:4326', 'EPSG:3857') AS geom
FROM sites;

CREATE INDEX idx_parcels_3857 ON parcels_3857 USING RTREE (geom);

-- 100 is now unambiguously metres because EPSG:3857 is a metric frame.
SELECT p.id, s.id
FROM parcels_3857 p
JOIN sites_3857 s ON ST_DWithin(p.geom, s.geom, 100);

The key behavioural change is twofold: both ST_Transform endpoints are explicit (the engine never guesses a source), and the distance predicate now runs in a metric frame so the 100 threshold means metres rather than degrees. Choosing the working CRS is a trade-off — a global equal-area frame distorts local distance, while a local UTM zone is accurate but only over its band; pick the frame that matches the spatial extent of the join. The same ST_DWithin semantics drive the proximity-filter and point-in-polygon patterns, so normalising the CRS once pays off across every downstream query.

Diagnostic Queries & Plan Validation #

DuckDB exposes no ST_SRID function — there is no inline SRID to read — so you diagnose CRS state from coordinate ranges and column types rather than metadata. Geographic data in EPSG:4326 must fall within ±180 / ±90; anything outside that window is projected (metric) data masquerading as lon/lat.

-- Bucket rows by whether they look geographic. A table that splits across
-- both buckets is carrying mixed CRSs and will join incorrectly.
SELECT
    ST_XMin(geom) BETWEEN -180 AND 180
      AND ST_YMin(geom) BETWEEN -90 AND 90 AS looks_geographic,
    COUNT(*)                      AS row_count,
    MIN(ST_XMin(geom))            AS min_x,
    MAX(ST_XMax(geom))            AS max_x
FROM target_table
GROUP BY looks_geographic;

Confirm the index is actually used and the transform sits in the projection phase of the plan rather than re-running per comparison:

-- Expect an RTREE_INDEX_SCAN node and a single PROJECTION applying ST_Transform.
-- An anti-pattern is ST_Transform appearing inside the join FILTER, which
-- reprojects on every probe instead of once during materialisation.
EXPLAIN ANALYZE
SELECT p.id, s.id
FROM parcels_3857 p
JOIN sites_3857 s ON ST_DWithin(p.geom, s.geom, 100);

Key thresholds when reading the output: the index scan should prune to a small candidate set (if RTREE_INDEX_SCAN reports nearly the full table, the predicate is not selective or the frames still differ), and the ST_Transform cost should appear once, not scale with the join cardinality. A useful one-liner for continuous monitoring is the bounding-box extent of a transformed table — a sudden jump in the extent between runs flags a layer that was ingested in the wrong frame:

-- Monitoring one-liner: extent should stay stable run-to-run for a given layer.
SELECT ST_XMin(b), ST_YMin(b), ST_XMax(b), ST_YMax(b)
FROM (SELECT ST_Extent(geom) AS b FROM parcels_3857);

Geometry Validation & Fallback Routing #

A transform inherits the invalidity of its input: self-intersecting or non-closed rings reproject into geometry that fails topology predicates further down the plan. Guard at ingestion with ST_IsValid, repair with ST_MakeValid, and reject coordinates that fall outside their declared frame before they reach the index.

-- Validate and repair before transforming; drop rows outside valid geographic
-- bounds so a mis-projected layer cannot poison the index or the join.
CREATE OR REPLACE TABLE parcels_clean AS
SELECT id,
       ST_Transform(ST_MakeValid(geom), 'EPSG:4326', 'EPSG:3857') AS geom
FROM parcels_raw
WHERE ST_IsValid(geom)
  AND ST_XMin(geom) BETWEEN -180 AND 180
  AND ST_YMin(geom) BETWEEN -90  AND 90;

For datasets beyond roughly 10M rows, transformation memory scales with vertex count and can breach memory_limit. When the spill directory alone is not enough, fall back to chunked execution: reproject in row ranges and append, so peak memory tracks the chunk rather than the whole table. The trade-offs behind this disk-versus-memory boundary are covered in in-memory vs disk storage.

-- Chunked fallback for OOM-prone transforms: peak memory tracks one batch,
-- not the whole table. Widen the BETWEEN range as headroom allows.
COPY (
    SELECT id, ST_Transform(ST_MakeValid(geom), 'EPSG:32633', 'EPSG:4326') AS geom
    FROM large_table
    WHERE id BETWEEN 1 AND 5000000
) TO '/output/batch_1.parquet' (FORMAT PARQUET);

Where the source CRS is genuinely unknown, route the layer to a quarantine view rather than guessing: expose only the validated, single-frame geometry to consumers (DuckDB has no GRANT/row-level security, so a curated view plus a read-only ATTACH is the boundary). Once the frame is identified, apply ST_Transform with explicit endpoints, or a manual ST_Affine correction when the transform parameters are known but no EPSG code maps to them. Carrying a normalised, validated frame all the way through is also what keeps a handoff to pandas correct — see DuckDB-to-GeoPandas sync for how CRS metadata must be re-attached on the GeoPandas side, since the WKB that crosses the boundary still carries no SRID.

Related #

See also:

• ST_Geometry vs WKB storage — why neither representation carries an inline SRID, and how the type choice affects index construction.
• Optimizing point-in-polygon queries in DuckDB — the ST_DWithin / ST_Intersects patterns that depend on a normalised CRS.

Up: CRS Mapping & Transformations in DuckDB Spatial · DuckDB Spatial Architecture & Fundamentals

External Reference Standards: GeoJSON coordinate and axis-order rules follow RFC 7946 (WGS 84 / EPSG:4326, [longitude, latitude]).