Coastal & Marine Spatial Analysis Pipelines

Debugging Axis-Order Inversion in pyproj

A coastal dataset that should sit off the Delaware coast suddenly plots in the Indian Ocean, or a UTM easting comes back as a nonsensical negative number, and nothing raised an exception. This is almost always an axis-order inversion: latitude and longitude have been silently swapped somewhere in a pyproj transform. It is one of the sharpest edges in the CRS alignment workflow for coastal GIS projects, because the transform “works” — it returns numbers, they are just the wrong numbers by thousands of kilometres. This page reproduces the swap with real coordinates, explains why the EPSG authority is technically correct and your intuition is not, and gives a production normaliser that pins axis order at the ingestion boundary so it can never drift again.

Same lon,lat Tuple, Two Axis Contracts in pyproj One coordinate tuple written as (lon=-75.09, lat=38.79) is sent through two pyproj Transformers. The top path uses the default (no always_xy): pyproj honours EPSG:4326 authority order, which is latitude first, so it reads the tuple as (lat=-75.09, lon=38.79). The transform succeeds but the easting/northing land far outside the survey bounding box, and the bounds sanity check raises. The bottom path uses always_xy=True: the tuple is read as lon,lat, the easting/northing fall inside the survey bbox, and the bounds check passes. A centre label reads "same numbers, different axis contract". input tuple (-75.09, 38.79) meant as lon, lat Transformer (default) honours authority order EPSG:4326 = lat, lon ✗ outside survey bbox bounds check raises Transformer(always_xy) forces traditional order reads tuple as lon, lat ✓ inside survey bbox bounds check passes same numbers, different axis contract

Root cause: the authority says lat first, your code assumes lon first

Every geodetic CRS defines an axis order. The EPSG registry declares EPSG:4326 (WGS 84 geographic) with latitude first, then longitude — the order a geodesist writes a position. Most GIS tooling, GeoJSON (RFC 7946), and the old PROJ.4 string world instead use longitude first, then latitude (the “x, y” convention). Both are internally consistent; they simply disagree about which number comes first.

Modern pyproj (PROJ 6+) made a deliberate, correctness-driven choice: a Transformer honours the authority axis order by default. So Transformer.from_crs("EPSG:4326", "EPSG:32618") expects its first argument to be latitude. Feed it a (lon, lat) tuple — which almost every coastal pipeline produces — and it transforms the wrong pair. The reproduction is exact:

from pyproj import Transformer

# Cape Henlopen, DE — written the way most pipelines store it: (lon, lat)
lon, lat = -75.0941, 38.7870

# DEFAULT transformer: honours EPSG:4326 authority order (lat, lon)
bad = Transformer.from_crs("EPSG:4326", "EPSG:32618")   # no always_xy
easting, northing = bad.transform(lon, lat)             # we pass (lon, lat) by habit
print(round(easting, 1), round(northing, 1))
# -> 40329954.6 4595537.9   <-- easting is ~40 million: nonsense, wrong hemisphere

An easting of forty million metres is physically impossible inside a UTM zone (valid eastings are roughly 100,000–900,000 m). The transform did exactly what it was told: it read -75.0941 as a latitude and 38.7870 as a longitude, and dutifully projected a point that does not exist near your survey. Nothing raised, because from PROJ’s perspective nothing was wrong.

This is the same class of silent corruption described in Fixing CRS Mismatch in QGIS and GeoPandas: the horizontal reference is technically defined but not the one you meant. Resolving which datum and axis convention a source truly uses is the discipline formalised in Geodetic Datums & EPSG Management.

Formally, the authority-ordered and traditional-ordered coordinate vectors are related by a swap matrix. If pauth=(φ,λ)\mathbf{p}_\text{auth} = (\varphi, \lambda) is (latitude, longitude) and pxy=(λ,φ)\mathbf{p}_\text{xy} = (\lambda, \varphi) is (longitude, latitude), then

pxy=S pauth,S=[0110],S=S−1.\mathbf{p}_\text{xy} = \mathbf{S}\,\mathbf{p}_\text{auth}, \qquad \mathbf{S} = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}, \qquad \mathbf{S} = \mathbf{S}^{-1}.

Because S\mathbf{S} is its own inverse, applying the wrong order once and never correcting it leaves the data cleanly transposed — which is exactly why a swapped dataset often looks almost plausible (a point near $45^\circ$N, $45^\circ$W maps onto another valid-looking location) instead of obviously broken.

Step-by-step fix with production code

1. Confirm the authority axis order for the CRS

Before changing anything, prove the diagnosis by reading what the CRS actually declares. pyproj exposes the axis metadata directly.

import logging

from pyproj import CRS

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(message)s",
    datefmt="%Y-%m-%dT%H:%M:%S",
)
logger = logging.getLogger(__name__)


def log_axis_order(epsg: str) -> tuple[str, ...]:
    """Log the declared axis order of a CRS and return the abbreviations."""
    crs = CRS.from_user_input(epsg)
    order = tuple(ax.abbrev for ax in crs.axis_info)
    directions = tuple(f"{ax.abbrev}:{ax.direction}" for ax in crs.axis_info)
    logger.info("CRS %s declares axis order %s (%s)", epsg, order, directions)
    return order


log_axis_order("EPSG:4326")   # -> ('Lat', 'Lon')  (latitude first)
log_axis_order("EPSG:32618")  # -> ('E', 'N')       (easting first)

Seeing ('Lat', 'Lon') printed is the confirmation: the default Transformer will treat your first argument as latitude.

2. Rebuild the Transformer with always_xy=True

The fix is a single keyword. always_xy=True tells pyproj to ignore authority order and accept/return coordinates in traditional GIS order — longitude/easting first, latitude/northing second — for both the input and output CRS.

from pyproj import Transformer

TARGET_CRS = "EPSG:32618"  # UTM Zone 18N (WGS84) — NE US coastal ops


def make_transformer(src: str, dst: str = TARGET_CRS) -> Transformer:
    """Build a lon,lat / easting,northing transformer with pinned axis order."""
    transformer = Transformer.from_crs(src, dst, always_xy=True)
    logger.info("Transformer %s -> %s built with always_xy=True", src, dst)
    return transformer


good = make_transformer("EPSG:4326")
easting, northing = good.transform(-75.0941, 38.7870)  # (lon, lat) in, (E, N) out
print(round(easting, 1), round(northing, 1))
# -> 492013.9 4293201.2   <-- a sane UTM 18N easting/northing near Cape Henlopen

The easting is now ~492 km into the zone and the northing ~4,293 km north of the equator — both physically reasonable for the mid-Atlantic coast.

3. Detect swapped data with a bounds sanity check

A keyword you have to remember is a keyword someone will forget. Encode the invariant as a check that raises when a transformed point escapes the known survey extent, and put it on the hot path so an inversion cannot pass silently.

from dataclasses import dataclass


@dataclass(frozen=True)
class SurveyBBox:
    """Expected projected extent of the survey, in the target CRS units (metres)."""
    min_e: float
    min_n: float
    max_e: float
    max_n: float

    def contains(self, easting: float, northing: float) -> bool:
        return (
            self.min_e <= easting <= self.max_e
            and self.min_n <= northing <= self.max_n
        )


# Delaware Bay approaches, UTM 18N, with generous margin.
DELAWARE_BBOX = SurveyBBox(min_e=400_000, min_n=4_200_000,
                           max_e=560_000, max_n=4_400_000)


def transform_checked(
    transformer: Transformer,
    lon: float,
    lat: float,
    bbox: SurveyBBox = DELAWARE_BBOX,
) -> tuple[float, float]:
    """Transform one point and raise if it escapes the survey extent."""
    easting, northing = transformer.transform(lon, lat)
    if not bbox.contains(easting, northing):
        raise ValueError(
            f"Transformed point ({easting:.1f}, {northing:.1f}) falls outside "
            f"survey bbox {bbox}. Likely lat/lon axis inversion at ingestion."
        )
    return easting, northing

4. Normalise axis order once, at the ingestion boundary

The durable fix is to stop letting individual call sites decide axis order. Force every incoming coordinate stream into one canonical (lon, lat) order the moment it enters the pipeline, log the coercion, and never revisit it downstream.

import numpy as np


def normalise_to_lonlat(
    coords: np.ndarray,
    source_epsg: str,
) -> np.ndarray:
    """Return an (N, 2) lon,lat array regardless of the source's authority order.

    `coords` is (N, 2) in the source CRS's *own* axis order. If the source
    declares latitude first, the columns are swapped so the output is always
    (lon, lat). Raises on malformed input rather than guessing.
    """
    arr = np.asarray(coords, dtype="float64")
    if arr.ndim != 2 or arr.shape[1] != 2:
        raise ValueError(f"coords must be (N, 2); got shape {arr.shape}")

    order = tuple(ax.abbrev for ax in CRS.from_user_input(source_epsg).axis_info)
    if order[0].lower().startswith("lat"):
        logger.info("Source %s is lat-first; swapping to lon,lat", source_epsg)
        arr = arr[:, ::-1]
    elif not order[0].lower().startswith(("lon", "e")):
        raise ValueError(f"Unrecognised axis order {order} for {source_epsg}")
    return np.ascontiguousarray(arr)

With a normaliser at the boundary, every internal transform can safely assume (lon, lat) and build its Transformer with always_xy=True. This is the ingestion contract the parent CRS alignment workflow expects every dataset to satisfy before any spatial operation runs.

Verification and acceptance test

A visual “it plots in the right ocean now” is not verification. Round-trip a known landmark and assert the projected coordinates land inside the survey bounding box, then wire it into CI so a future refactor that drops always_xy fails loudly.

import pytest


def test_landmark_lands_in_survey_bbox():
    # Landmark: Cape Henlopen lighthouse approach, lon/lat (WGS84).
    lon, lat = -75.0941, 38.7870
    transformer = make_transformer("EPSG:4326")

    easting, northing = transform_checked(transformer, lon, lat)

    assert DELAWARE_BBOX.contains(easting, northing)
    # Independent hard bounds: any valid UTM easting is 100k–900k m.
    assert 100_000 < easting < 900_000, f"easting {easting:.1f} not a valid UTM value"


def test_default_transformer_would_swap():
    # Guard/regression: the un-fixed transformer must fail the bounds check,
    # proving the check actually catches an inversion.
    lon, lat = -75.0941, 38.7870
    bad = Transformer.from_crs("EPSG:4326", TARGET_CRS)  # no always_xy
    with pytest.raises(ValueError, match="axis inversion"):
        transform_checked(bad, lon, lat)


def test_normaliser_swaps_lat_first_source():
    latlon = np.array([[38.7870, -75.0941]])       # (lat, lon) as EPSG:4326 declares
    out = normalise_to_lonlat(latlon, "EPSG:4326")
    assert out[0, 0] == pytest.approx(-75.0941)     # lon now first
    assert out[0, 1] == pytest.approx(38.7870)

As a CLI cross-check outside Python, PROJ’s own cs2cs/projinfo report the axis order the library will use, independent of your code:

projinfo EPSG:4326 | grep -iA1 "Axis"     # shows "Latitude" listed before "Longitude"

Edge cases and gotchas

  • WMS/WMTS 1.3.0 flip. OGC bumped WMS from 1.1.1 to 1.3.0 specifically to honour authority axis order, so a BBOX for EPSG:4326 in WMS 1.3.0 is lat,lon while the same request in 1.1.1 is lon,lat. When stitching tiled coastal basemaps, pin the version and swap the bbox accordingly, or tiles land transposed even though every individual transform is correct.
  • GeoPandas vs raw shapely. GeoDataFrame.to_crs() handles axis order correctly through its attached CRS, but the moment you drop to bare shapely geometries and feed their .x/.y into a hand-built Transformer, you are back to choosing the order yourself — always pass always_xy=True there. GeoJSON output is always lon,lat by RFC 7946 regardless of the CRS’s authority order.
  • PostGIS ST_FlipCoordinates. If upstream data was ingested swapped, ST_FlipCoordinates(geom) transposes X and Y in place — useful for a one-off repair, dangerous if run twice (it is its own inverse, matching the swap matrix S\mathbf{S} above). Fix the ingestion normaliser rather than sprinkling flips through queries.
  • Deprecated proj4 strings. A legacy +proj=longlat +datum=WGS84 PROJ.4 string implies traditional lon,lat order and will not trigger the authority-order behaviour, so mixing proj4-string sources with EPSG-code sources in the same pipeline yields inconsistent ordering. Normalise every source to an EPSG code (or full WKT2) at ingestion so axis handling is uniform.

Authoritative axis-order behaviour is documented by the pyproj gotchas guide and the PROJ axis-order documentation.

Up: CRS Alignment for Coastal GIS Projects