CRS in DuckDB: Why ST_Transform Moves Your Data to the Wrong Place
Problem statement
Three measured results, all from DuckDB 1.5.5 with the spatial extension, all silent:
-- 1. the transform that lands in the wrong hemisphere
select st_astext(st_transform(st_point(-0.1276, 51.5072),
'EPSG:4326', 'EPSG:3857'));
POINT (5733755.276187301 -14204.378766821068) -- London, apparently
-- 2. the join that returns nothing, with no error
select count(*) from a join b on st_intersects(a.geom, b.geom);
0
-- 3. the distance that is 17.5% too long
select st_distance_sphere(st_point(-0.1276, 51.5072), -- London
st_point(2.3522, 48.8566)); -- Paris
403552.4605518389 -- true: 343,530 mNone of these raise. All three come from the same root: DuckDB's spatial extension follows the EPSG database's axis order, and almost every other tool in the Python GIS stack does not.
Quick answer
Pass always_xy := true to every ST_Transform, and never pass a lon/lat point to the sphere or spheroid functions:
-- correct: longitude first, as every GIS file stores it
select st_astext(st_transform(st_point(-0.1276, 51.5072),
'EPSG:4326', 'EPSG:3857', always_xy := true));
POINT (-14204.367025221705 6711506.705400525) -- correct Web Mercator
select st_astext(st_transform(st_point(-0.1276, 51.5072),
'EPSG:4326', 'EPSG:27700', always_xy := true));
POINT (530043.1949814068 180358.20862000482) -- correct British National GridFor distances, swap the coordinates or use a projected CRS:
-- ST_Distance_Sphere expects POINT(latitude longitude)
select st_distance_sphere(st_point(51.5072, -0.1276),
st_point(48.8566, 2.3522));
343529.8654949107 -- correct
Step-by-step solution
1. Know where the CRS lives
DuckDB's spatial extension carries the CRS on the column type, not on each value:
select typeof(geom) from st_read('provinces.shp') limit 1;
-- GEOMETRY('EPSG:4326')
select typeof(st_point(0, 0));
-- GEOMETRY <- constructed geometries have no CRS
select typeof(st_transform(geom, 'EPSG:4326', 'EPSG:3857')) from st_read('x.shp') limit 1;
-- GEOMETRY <- the transform drops itSo the CRS is present when a file supplies it and absent as soon as you compute anything. Treat the type annotation as a helpful label, not as a tracked property.
2. Understand the axis-order problem
EPSG:4326 is officially defined with latitude first. Almost every file format, every GeoJSON, every POINT(x y) and every points_from_xy(lon, lat) in the Python ecosystem uses longitude first.
PROJ can follow either convention, and DuckDB's ST_Transform defaults to the authority's — latitude first. That is why transforming a normal lon/lat point without always_xy produces a coordinate in the wrong place: it interpreted your longitude as a latitude.
Measured for London (−0.1276, 51.5072) to EPSG:3857:
always_xy default (false) POINT (5733755, -14204) wrong
always_xy := true POINT (-14204.4, 6711506.7) correct
The second matches what pyproj and GeoPandas produce, because they default to always_xy=True.
3. Watch for the silent zero-row join
Joining geometries in different coordinate systems is not an error in DuckDB — it is arithmetic on numbers that happen not to overlap:
with a as (select geom from st_read('places.shp') limit 10),
b as (select st_transform(geom, 'EPSG:4326', 'EPSG:3857') g
from st_read('places.shp') limit 10)
select count(*) from a join b on st_intersects(a.geom, b.g);
-- 0Ten points joined against the same ten points, and the answer is zero. PostGIS raises on mixed SRIDs; DuckDB does not, so the check has to be yours.
4. Never pass lon/lat to the sphere and spheroid functions
ST_Distance_Sphere, ST_Distance_Spheroid, ST_Area_Spheroid and ST_Perimeter_Spheroid all expect POINT(latitude longitude). Measured:
input result truth
POINT(lon lat) London → Paris 403,552 m 343,530 m (+17.5%)
POINT(lat lon) London → Paris 343,530 m 343,530 m correct
ST_Area_Spheroid on a lon/lat polygon NaN
ST_Area_Spheroid on a lat/lon polygon 269,154,549,884 m²
The NaN is the merciful case — it is at least visible. The 17.5% distance error is not, and it is exactly the kind of number that survives into a report.
5. Check ST_Area and ST_Length units
ST_Area and ST_Length are planar and use whatever units the coordinates are in. On a lon/lat polygon that means square degrees:
select st_area(geom) from st_read('provinces.shp') where name = 'Colorado';
-- 27.998104200139107 <- square degrees, not square metresA square degree is not a unit of area anywhere on Earth. Reproject before measuring, and use ST_Area_Spheroid — with lat/lon input — if you need geodesic area without projecting.
6. Adopt a convention and enforce it
The convention that survives contact with the rest of the stack:
-- 1. always pass always_xy := true
-- 2. keep everything in one projected CRS for analysis
-- 3. transform to 4326 only at the very edge, for output
create table analysis as
select st_transform(geom, 'EPSG:4326', 'EPSG:27700', always_xy := true) as geom, *
from st_read('input.shp');Once every geometry is in one projected CRS, ST_Area, ST_Length, ST_Distance and ST_Buffer all return metres and the axis question disappears.
Code examples
Example 1 — a transform wrapper that cannot be called wrongly
def transform_sql(column, source_crs, target_crs):
"""Always always_xy. There is no case in a lon/lat file pipeline where the
authority axis order is what you want."""
return (f"st_transform({column}, '{source_crs}', '{target_crs}', "
f"always_xy := true)")
def reproject_table(con, source, target_table, source_crs, target_crs,
geom_column="geom"):
con.execute(f"""
create or replace table {target_table} as
select * exclude {geom_column},
{transform_sql(geom_column, source_crs, target_crs)} as {geom_column}
from {source}
""")
sample = con.execute(f"select st_astext({geom_column}) from {target_table} "
f"limit 1").fetchone()[0]
print(f"{target_table}: {sample[:60]}")
print(f" CRS is no longer tracked on the type — record it: {target_crs}")Example 2 — asserting the CRS before a join
def assert_same_crs(con, *tables, geom_column="geom"):
"""DuckDB will not raise on a CRS mismatch. This will."""
types = {}
for table in tables:
t = con.execute(f"select typeof({geom_column}) from {table} limit 1").fetchone()
types[table] = t[0] if t else "empty"
distinct = set(types.values())
for table, t in types.items():
note = " <- no CRS on the type" if t == "GEOMETRY" else ""
print(f"{table:24} {t}{note}")
if len(distinct) > 1:
raise ValueError(
f"geometry column types differ: {types}. A join across these returns "
f"zero rows without an error — reproject one side first.")
if distinct == {"GEOMETRY"}:
print("! neither side declares a CRS; the type check cannot help you here")
return TrueThe second warning matters as much as the exception. Once a transform has stripped the CRS from both sides, the types match and prove nothing.
Example 3 — a sanity check with a known coordinate
KNOWN = {
# (lon, lat) -> expected easting, northing after transform
"EPSG:27700": ((-0.1276, 51.5072), (530043, 180358), 5), # central London
"EPSG:3857": ((-0.1276, 51.5072), (-14204, 6711507), 5),
}
def verify_transform(con, target_crs):
"""One known point catches an axis-order mistake immediately."""
(lon, lat), (want_x, want_y), tolerance = KNOWN[target_crs]
got = con.execute(f"""
select st_x(p), st_y(p) from (
select st_transform(st_point({lon}, {lat}), 'EPSG:4326', '{target_crs}',
always_xy := true) as p)
""").fetchone()
dx, dy = abs(got[0] - want_x), abs(got[1] - want_y)
ok = dx < tolerance and dy < tolerance
print(f"{target_crs}: got ({got[0]:,.0f}, {got[1]:,.0f}), "
f"expected ({want_x:,}, {want_y:,}) {'ok' if ok else 'AXIS ORDER WRONG'}")
if not ok:
print(" → add always_xy := true, or your coordinates are already lat/lon")
return okA single known point is the cheapest possible regression test for this whole class of bug, and it belongs in the test suite of any pipeline that transforms coordinates.
Explanation
Why the axis-order ambiguity exists at all
The EPSG registry defines EPSG:4326 with latitude as the first axis, because that is how geodesists write coordinates. GIS software, GeoJSON, most APIs and every (x, y) convention put longitude first, because that is how a Cartesian plane works.
Both are defensible and both are in wide use. PROJ therefore supports both and offers always_xy to force the GIS convention. pyproj and GeoPandas default to always_xy=True; DuckDB's ST_Transform does not, which is why the same transform gives different answers in two tools on one machine.
Why the wrong result is so far away rather than slightly off
Swapping London's coordinates makes latitude −0.1276 and longitude 51.5072 — a point in the Gulf of Guinea, five thousand kilometres east of the meridian. Web Mercator then maps that to an easting of 5.7 million metres and a northing near zero.
A slight error would be easier to miss and harder to cause. This one is visible on any map, which is the single mercy of the whole situation: the failure is loud once anything plots it.
Why the zero-row join is the more dangerous failure
The transform error is visible. A join returning zero rows looks like a legitimate answer: no points fell in any polygon. Nothing on the screen distinguishes "correctly empty" from "coordinate systems do not overlap".
PostGIS prevents this by storing an SRID with each geometry and raising on a mismatch. DuckDB's type annotation is not enforced across operations, so a count(*) of zero after a join deserves an explicit CRS check before it is believed.
Why the sphere functions take latitude first
ST_Distance_Sphere is documented as taking latitude in the first coordinate, following the same geodetic convention as the EPSG axis order. It is internally consistent and externally surprising, because ST_Point(x, y) — the function you used to build the point — is longitude first.
So one line constructs a lon/lat point and the next line reads it as lat/lon, inside one query, with no error. The measured consequence is a distance 17.5% too long; the fix is either to swap the arguments or, better, to work in a projected CRS where ST_Distance is plain metres.
Edge cases or notes
always_xy := trueon every transform. There is no case in a normal GIS pipeline where the authority order is what you want.ST_Transformdrops the CRS from the result type, so the annotation cannot be relied on downstream.- A shapefile with no
.prjyields a plainGEOMETRY— no CRS, no warning. ST_Areaon lon/lat returns square degrees. Colorado came out as 28.0.ST_Area_Spheroidon lon/lat returnsNaNrather than a wrong number, which is the one helpful failure here.ST_Bufferon lon/lat buffers in degrees, producing an ellipse that shrinks towards the poles.- Record the CRS in the table or column name if you cannot rely on the type.
- Test with one known coordinate. London to EPSG:27700 should be about (530043, 180358).
Internal links
- DuckDB spatial explained: a spatial database that is just a file — the wider engine
- Fixing coordinates in the wrong order in an API response — the same trap elsewhere
- Projected versus geographic CRS explained — why degrees are not metres
- My points plot in the ocean: latitude and longitude swapped — the classic symptom
- How to reproject spatial data in GeoPandas — the same operation, other defaults
- Fixing a CRS mismatch in GeoPandas — the error GeoPandas gives you and DuckDB does not
- How to run a spatial join in DuckDB — where the zero-row failure appears
- How to measure distance accurately — geodesic distance done properly
FAQ
Why does ST_Transform put my point in the wrong place?
Because it follows the EPSG axis order — latitude first — by default. Pass always_xy := true and it matches what GeoPandas and pyproj produce.
Does DuckDB track the CRS of a geometry?
Partly. ST_Read puts it on the column type as GEOMETRY('EPSG:4326'), but constructed geometries have none and ST_Transform drops it from its result.
Why does my spatial join return zero rows?
Very often a CRS mismatch. DuckDB does not raise on mixed coordinate systems — the geometries simply do not overlap. Check the column types before believing an empty result.
Why is ST_Distance_Sphere giving the wrong distance?
It expects POINT(latitude longitude). With a normal lon/lat point, London to Paris came back as 403,552 m instead of 343,530 m — 17.5% too far.
Why does ST_Area return a tiny number?
Because it is planar and your coordinates are degrees, so the answer is in square degrees. Colorado measured 28.0. Reproject before measuring area.
What convention should I adopt?
always_xy := true on every transform, one projected CRS for the whole analysis, and a conversion to EPSG:4326 only at the output edge.