Detecting Axis Order Issues in CRS Definitions

An axis-order mistake moves a point across the world, and it is the one CRS error that is trivially detectable and endemic anyway — because most code assumes longitude comes first and most authoritative definitions say latitude does. This guide, part of working with EPSG and WKT2 CRS definitions, reads the declared order from the definition instead of assuming it, detects a transposition from the coordinates themselves, and states the one convention a pipeline should adopt internally.

Why There Are Two Orders #

The authority definition of a geographic CRS states its axis order, and for most geographic systems — including EPSG:4326 and EPSG:6318 — that order is latitude then longitude. Most software, most file formats and most APIs use longitude then latitude, because that matches the x-then-y convention of everything else in graphics and GIS. Neither is wrong; they are different conventions, and the mismatch lives at every boundary between an authority definition and a piece of code.

Figure — the authority says latitude first; nearly all code says longitude first.

The result is that “EPSG:4326 coordinates” is ambiguous unless the order is stated, and that the two readings put a point in completely different places. A position at latitude 45, longitude −122 read in the other order is latitude −122, which is not a latitude at all — the useful case, because it is detectable — while a position at latitude 40, longitude 20 read the other way is a perfectly valid location in a different hemisphere.

Complete Runnable Implementation #

from __future__ import annotations

from dataclasses import dataclass

import numpy as np
from pyproj import CRS


@dataclass(frozen=True)
class AxisOrder:
    """The declared axis order of a CRS, read rather than assumed."""

    first: str                # abbreviation, e.g. "Lat" or "E"
    second: str
    is_lat_first: bool
    units: tuple[str, str]

    @classmethod
    def of(cls, crs: CRS) -> "AxisOrder":
        axes = crs.axis_info
        if len(axes) < 2:
            raise ValueError(f"{crs.name} declares fewer than two axes")
        a, b = axes[0], axes[1]
        return cls(
            first=a.abbrev, second=b.abbrev,
            is_lat_first=a.direction.lower() in ("north", "south"),
            units=(a.unit_name, b.unit_name),
        )


def looks_transposed(crs: CRS, first: np.ndarray, second: np.ndarray) -> bool:
    """Heuristic transposition test from the coordinate values themselves.

    Two independent signals: a latitude outside +/-90 is impossible, and for a
    projected CRS a northing and an easting differ in magnitude in a way that is
    usually decisive. Neither is proof; both are worth acting on.
    """
    if crs.is_geographic:
        order = AxisOrder.of(crs)
        lat = first if order.is_lat_first else second
        return bool(np.any(np.abs(lat) > 90.0))
    # Projected: in a zone-based system the northing is far larger than the easting.
    med_first, med_second = float(np.median(np.abs(first))), float(np.median(np.abs(second)))
    return med_first > 10.0 * max(med_second, 1.0)


def to_xy(crs: CRS, a: np.ndarray, b: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Normalise a coordinate pair to the pipeline's internal (x, y) convention.

    Internally this codebase is always longitude/easting first. The conversion
    happens exactly once, at ingress, and is recorded — the alternative is every
    function having its own opinion.
    """
    order = AxisOrder.of(crs)
    return (b, a) if order.is_lat_first else (a, b)

Parameter Reference #

Worked Example #

import numpy as np
from pyproj import CRS

crs = CRS.from_epsg(6318)                    # NAD83(2011) geographic
order = AxisOrder.of(crs)
print(order.first, order.second, order.is_lat_first)
# Lat Lon True

a = np.array([45.5, 45.6])                   # as stored, authority order
b = np.array([-122.6, -122.7])
x, y = to_xy(crs, a, b)
print(x[:1], y[:1])
# [-122.6] [45.5]

print(looks_transposed(crs, np.array([-122.6]), np.array([45.5])))
# True   -- a "latitude" of -122.6 cannot be right

Figure — three detection tests, from the cheapest heuristic to the definitive one.

Validation Check #

def assert_within_area_of_use(crs: CRS, x: np.ndarray, y: np.ndarray) -> None:
    """Range-check against the CRS area of use — the definitive transposition test."""
    area = crs.area_of_use
    if area is None:
        return
    inside = ((x >= area.west) & (x <= area.east)
              & (y >= area.south) & (y <= area.north))
    assert inside.all(), (
        f"{int((~inside).sum())} coordinate(s) fall outside the area of use of "
        f"{crs.name}; check the axis order before anything else"
    )

The area-of-use check is stronger than either heuristic and costs four comparisons per point. It also catches wrong-zone and wrong-CRS errors, which is why it belongs in the ingress guard rather than in a debugging session.

Common Mistakes #

Assuming longitude first because the file “looks like” it. A file whose first column happens to be in the range of a valid latitude tells you nothing, because longitudes in that range exist too. Read the declared order from the CRS, and where the file format has no CRS, treat the order as an explicit, recorded assumption.

Normalising in more than one place. Two functions that each swap the order produce coordinates that are correct in one and transposed in the other, and the bug appears only in whichever path calls both. Normalise once, at ingress, and let everything downstream assume the internal convention.

Trusting a heuristic as proof. Both heuristics above are one-sided: they detect an obvious transposition and stay silent on a subtle one. Use them to catch the obvious case early, and the area-of-use check as the real gate — the ordering argument made in handling CRS mismatches in cadastral datasets.

Where the Convention Is Decided in a Pipeline #

Axis order causes trouble in proportion to how many places decide it, so the design goal is to have exactly one. Four boundaries are where the decision gets made, usually implicitly.

Figure — four boundaries where the convention gets decided, usually implicitly.

At ingress, when a file is read. This is the right place: read the declared order from the CRS, normalise to the internal convention, and record what was done. Every function downstream then has one assumption instead of a question.

In the transformer, via the axis-order flag. Setting it consistently for every transformer the pipeline builds makes the library agree with the internal convention rather than with the authority definition, which is what you want once ingress has normalised.

At egress, when a file is written. The output format may require the authority order — a WKT2-carrying container will be read according to what it declares — so the conversion back happens once, at the boundary, and is the mirror of ingress.

In visualisation and debugging code, which is where the convention is most often decided by accident. A plotting call that takes x and y will happily accept latitude and longitude in that order and produce a map that is transposed, and because it is only a diagnostic view nobody treats it as a defect — until a decision is made from it.

Writing the convention down in one sentence at the top of the module that does the normalising is worth more than any amount of defensive checking further in.

Frequently Asked Questions #

Should I use always_xy=True on every transformer? #

Yes, if the pipeline’s internal convention is longitude first — and it should be, because that is what NumPy arrays, file formats and plotting libraries expect. The flag makes the transformer speak that convention regardless of what the CRS declares, which turns the axis order from a per-CRS property into a single project-wide decision.

Does axis order affect projected CRSs too? #

It does, and it is less often noticed because easting and northing usually differ enough in magnitude to look obviously wrong. Some national grids declare northing first, though, and in a system where both values are of similar size a transposition is entirely plausible-looking. Read the declaration for projected systems as well.

What about the third axis? #

Height, where present, is declared with its own direction — up or down — and a downward-positive axis exists in some engineering and bathymetric systems. It is rarer than the horizontal problem and produces a sign error rather than a transposition, which is easier to spot and just as wrong.

How does this interact with GeoJSON? #

GeoJSON fixes the order as longitude then latitude on WGS84 by specification, which is one of the few places where the ambiguity does not exist. That makes it a useful interchange format for the specific case it covers and a poor one for anything projected, as geospatial file formats and CRS metadata sets out.

Does axis order affect a transformation’s accuracy? #

Not its accuracy — its correctness. A transposed pair is not a slightly worse coordinate, it is a different place, usually thousands of kilometres away. That distinction matters when triaging a residual report: anything at the scale of a real transformation error is not an axis problem, and anything at continental scale almost certainly is.

The summary is short: read the declared order rather than assuming it, normalise once at ingress, and range-check against the area of use — three habits that between them retire the entire class of problem.

Related #

• Working with EPSG and WKT2 CRS definitions — the parent topic
• Parsing WKT2:2019 CRS strings in Python
• Handling CRS mismatches in cadastral datasets
• Reading and writing .prj files with Python
• Setting up high-precision coordinate reference systems