Fixing PyProj CRS Transformation Errors: Axis Order & EPSG Validation

A pyproj transform that raises CRSError on the first call is annoying but honest; the one that returns numbers off by thousands of kilometres without complaint is the one that ships bad data. This guide isolates both failure modes — hard CRSError/ProjError exceptions and silent latitude/longitude swaps — and gives a deterministic fix for each. It is for anyone wiring Coordinate Systems with PyProj into a production pipeline within Mastering Core Geospatial Python Libraries, and it pairs with Choosing the Right UTM Zone Automatically in Python once you know which projected CRS you should be transforming into.

Why This Approach / What Goes Wrong

PROJ 6 rewrote how coordinate reference systems are described. The old behaviour — treating every CRS as if it emitted (x, y) in the order you happened to want — is gone. Modern PROJ (bundled with pyproj>=3.0) honours the axis order declared in the EPSG registry, and for geographic CRSs like EPSG:4326 that order is latitude first, longitude second. Two distinct problems fall out of this.

The first is the hard failure. Passing a raw integer, a hand-written +proj= string with a typo, or a deprecated datum alias to Transformer.from_crs() raises CRSError (the CRS could not be built) or ProjError (the transformation pipeline could not be constructed). These are loud and easy to fix once you know to validate the CRS objects before you build the transformer.

The second is the silent swap, and it is far more dangerous. A Transformer built from EPSG:4326 expects its input as (lat, lon) by default. If your data — like almost every GeoJSON, shapefile, or API response in the wild — is ordered (lon, lat), PyProj will happily transform it, feeding your longitude in where it expects latitude. Nothing raises. You just get an easting and northing for the wrong place. The fix is always_xy=True, which forces the traditional GIS (lon, lat) / (x, y) convention on both ends regardless of what the EPSG registry says. This single flag is the most common cause of "my coordinates are in the ocean" bug reports on this site; the same trap is covered from the reprojection angle in Coordinate Reference System Transformations.

There is a third class, and it sits between the two. PROJ chooses an operation from a ranked candidate list, and when the highest-accuracy candidate needs a datum-shift grid that is not present on disk it quietly falls back to a lower-ranked one — a Helmert approximation, or in the worst case a "ballpark" null offset that pretends the two datums coincide. Nothing in the returned tuple distinguishes a 0.01 m grid-based result from a 2 m ballpark. The tell-tale symptom is a round-trip that closes to nine decimal places while the coordinates disagree with a surveyed control point, because the same wrong operation ran in both directions and its error cancelled itself out. Chasing that class down is a separate exercise, covered in datum shifts and transformation grids in PyProj; what matters here is knowing it exists, so that you never read a clean round-trip as proof of accuracy — only as proof of self-consistency.

always_xy decides where the Berlin sensor lands One WGS 84 sensor at longitude 13.404954, latitude 52.520008 feeds two transformers. The top path with always_xy=False (the default) makes PROJ interpret the first value as latitude and the second as longitude, producing 4847457, 1852284 — roughly 4000 km off, with no exception raised. The bottom path with always_xy=True keeps the traditional (lon, lat) order and returns 391776, 5820073, the correct ETRS89 / UTM zone 33N easting and northing for Berlin. One flag decides where the point lands Berlin sensor input order (lon, lat) 13.404954, 52.520008 default — silently wrong Transformer.from_crs(...) always_xy=False (default) PROJ reads as (lat, lon) 13.40 → latitude 52.52 → longitude 4 847 457, 1 852 284 ✗ ~4000 km off explicit — correct Transformer.from_crs( …, always_xy=True) keeps (lon, lat) order 13.40 → longitude 52.52 → latitude 391 776, 5 820 073 ✓ UTM 33N · Berlin
always_xy decides where the Berlin sensor lands

Prerequisites

conda install -c conda-forge "pyproj=3.4.*" "geopandas=0.14.*"

Install PyProj from conda-forge, not a bare pip wheel layered on top of a conda GDAL — a mismatched PROJ database is the classic source of CRSError: Invalid projection on codes that are perfectly valid. Confirm the versions with python -c "import pyproj; print(pyproj.__version__, pyproj.proj_version_str)"; you want PROJ 8+ under PyProj 3.4+.

Step-by-Step Implementation

The worked example transforms an air-quality sensor in Berlin from WGS 84 lon/lat into a projected, metre-based CRS suitable for distance and area work. Berlin sits in UTM zone 33N, so the correct target is EPSG:25833 (ETRS89 / UTM zone 33N) — a proper projected CRS, not Web Mercator, which distorts metric measurement badly at high latitudes.

1. Build and validate the CRS objects before you build the transformer. CRS.from_epsg() hits the PROJ registry and rejects bad codes immediately, instead of deferring the error to transform time.

from pyproj import CRS

# Sensor metadata says coordinates are WGS 84 lon/lat; target is metric UTM 33N.
source_crs = CRS.from_epsg(4326)      # WGS 84 (geographic, degrees)
target_crs = CRS.from_epsg(25833)     # ETRS89 / UTM zone 33N (metres)

# Fail fast on a malformed or unknown definition
assert source_crs.is_valid and target_crs.is_valid, "Invalid CRS definition"
assert target_crs.axis_info[0].unit_name == "metre", "Target CRS is not metric"

Prefer CRS.from_epsg(25833) over the shorthand integer 25833 you might pass elsewhere: it forces validation now and gives you a CRS object you can inspect (.axis_info, .is_geographic, .to_wkt()).

2. Construct the Transformer with an explicit axis order. This is the line that fixes the silent swap.

from pyproj import Transformer

# always_xy=True => inputs and outputs are (lon, lat) / (easting, northing),
# regardless of the lat-first axis order EPSG:4326 declares.
transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)

3. Transform using the (lon, lat) -> (x, y) convention.

sensor_lon, sensor_lat = 13.404954, 52.520008   # Berlin, WGS 84
easting, northing = transformer.transform(sensor_lon, sensor_lat)
print(f"{easting:.2f}, {northing:.2f}")
# 391776.16, 5820073.12

4. See the failure for yourself. Rebuild the transformer without always_xy and feed it the same (lon, lat) tuple your data actually carries. PROJ reads the first value as latitude, and the result lands nowhere near Berlin — with no exception raised.

wrong = Transformer.from_crs(source_crs, target_crs)   # always_xy defaults to False
bad_e, bad_n = wrong.transform(sensor_lon, sensor_lat) # lon fed where lat is expected
print(f"{bad_e:.2f}, {bad_n:.2f}")
# 4847457.95, 1852284.37   <- silently wrong, ~4000 km off

5. Reuse the transformer; never rebuild it per point. Each Transformer.from_crs() call triggers a PROJ database lookup and can trigger grid-file loading. Instantiate once and cache it — module-level, or via functools.lru_cache keyed on the EPSG pair.

Rebuilding the Transformer per point versus caching it once Two work timelines for the same batch of sensor coordinates. The top lane rebuilds the Transformer inside the loop, so a PROJ database lookup precedes every single transform call and repeats for every point. The bottom lane pays one lookup and pipeline build, then reuses the cached Transformer for every remaining point, so only cheap transform steps follow. Build the Transformer once, then reuse it every from_crs() call re-queries the PROJ database and may reload datum grids per-point rebuild from_crs() in the loop PROJ lookup tf PROJ lookup tf PROJ lookup tf … repeated per point N × lookup build once, cache lru_cache on the EPSG pair PROJ lookup + pipeline build tf tf tf tf tf tf tf tf … reused 1 × lookup Cache on the EPSG pair, then hand whole arrays to transform() — one pipeline, many points
The database lookup, not the arithmetic, is what makes a per-point rebuild slow — caching on the EPSG pair collapses N lookups to one.
from functools import lru_cache
from pyproj import Transformer

@lru_cache(maxsize=32)
def get_transformer(src_epsg: int, dst_epsg: int) -> Transformer:
    return Transformer.from_crs(src_epsg, dst_epsg, always_xy=True)

tf = get_transformer(4326, 25833)
# Vectorised: pass arrays, not a Python loop, for large sensor batches
lons = [13.404954, 13.388860, 13.428555]
lats = [52.520008, 52.517037, 52.523430]
xs, ys = tf.transform(lons, lats)

If you are moving whole layers rather than coordinate arrays, GeoPandas .to_crs() delegates to PyProj and applies always_xy semantics internally, so you only manage axis order by hand when you extract raw coordinates. For layers too large to reproject in memory, see Reprojecting Large Datasets Without Memory Errors.

6. Decide whether a bad coordinate should raise or be flagged. By default PROJ writes an infinity sentinel into any output slot it could not compute and returns normally, which is exactly the behaviour that lets a corrupt row travel three stages downstream before anyone notices. errcheck=True converts that sentinel into a ProjError at the point of failure.

import numpy as np
from pyproj import Transformer

tf = Transformer.from_crs(4326, 25833, always_xy=True)

# One corrupt row: a latitude of 152.52 cannot exist
lons = np.array([13.404954, 13.388860, 13.428555])
lats = np.array([52.520008, 52.517037, 152.523430])

# Default behaviour — the bad row comes back as inf, the good rows are untouched
easting, northing = tf.transform(lons, lats)
print(np.isfinite(easting))
# [ True  True False]

# errcheck=True raises on the first coordinate PROJ cannot handle
try:
    tf.transform(lons, lats, errcheck=True)
except Exception as exc:                      # pyproj.exceptions.ProjError
    print(type(exc).__name__, "->", exc)
# ProjError -> Error ... latitude or longitude exceeded limits

Pick the mode that matches the caller. A request handler transforming a single user-supplied coordinate wants errcheck=True, because an exception maps cleanly onto a 400 response and there is nothing to salvage. A nightly batch of ten million sensor fixes wants the default, because raising discards the 9,999,999 rows that were fine; there you transform the whole array, mask with np.isfinite, quarantine the offending rows, and carry on. The check itself is a single comparison per coordinate, so the cost is irrelevant next to the projection maths — the choice is about error semantics, not speed.

7. Pin the operation to your study area when several candidates exist. For CRS pairs served by regional grids — NAD27 to NAD83, OSGB36 to ETRS89, many national realizations — PROJ may know a dozen operations with different footprints and accuracies. Without a hint it ranks them across the whole overlap of the two systems, which can select a nationwide low-accuracy operation when a precise regional one covers your data. area_of_interest (pyproj 3.0+) narrows the ranking.

from pyproj import Transformer
from pyproj.aoi import AreaOfInterest

# The AOI is ALWAYS given in lon/lat degrees, regardless of always_xy
berlin_extent = AreaOfInterest(
    west_lon_degree=13.0, south_lat_degree=52.3,
    east_lon_degree=13.8, north_lat_degree=52.7,
)

tf_scoped = Transformer.from_crs(
    "EPSG:4326", "EPSG:25833", always_xy=True, area_of_interest=berlin_extent
)

tf_scoped.transform(13.404954, 52.520008)     # operation is chosen on first use
print(tf_scoped.description)
# Inverse of ETRS89 to WGS 84 (1) + UTM zone 33N

Two details bite people here. The AOI is expressed in longitude/latitude degrees even when always_xy=True is set — that flag governs coordinate arguments, not the AOI constructor. And transformer.description reports unavailable until proj_trans is called until you have pushed at least one coordinate through, because operation selection is deferred to the first transform; print it after a call, never before.

Verification

Prove the fix with a round-trip and a bounds check. Expected console output is shown inline as comments.

import math
from pyproj import Transformer

fwd = Transformer.from_crs(4326, 25833, always_xy=True)
inv = Transformer.from_crs(25833, 4326, always_xy=True)

sensor_lon, sensor_lat = 13.404954, 52.520008
easting, northing = fwd.transform(sensor_lon, sensor_lat)

# 1. Forward result matches the known-good UTM 33N easting/northing for Berlin
assert math.isclose(easting, 391776.16, abs_tol=0.5), f"Unexpected easting: {easting}"
assert math.isclose(northing, 5820073.12, abs_tol=0.5), f"Unexpected northing: {northing}"

# 2. Round-trip back to lon/lat returns the original within floating-point tolerance
lon_rt, lat_rt = inv.transform(easting, northing)
assert math.isclose(lon_rt, sensor_lon, abs_tol=1e-6)
assert math.isclose(lat_rt, sensor_lat, abs_tol=1e-6)
print(f"round-trip OK: {lon_rt:.6f}, {lat_rt:.6f}")
# round-trip OK: 13.404954, 52.520008

A clean round-trip that lands back on the input is the strongest single signal that both the axis order and the CRS pair are correct. If the round-trip drifts by whole metres or the forward result misses the assertion, you almost certainly swapped an axis or targeted the wrong zone.

A single point, though, only proves that one point worked. Batch pipelines need an assertion over the whole array, and the two cheapest ones are a finiteness mask and a zone-extent check. UTM eastings are defined with a 500,000 m false easting and a 6° wide zone, so every legitimate easting in any UTM zone falls between roughly 166,000 m and 834,000 m; a value outside that band means the point is not in the zone you targeted, whatever the transform returned.

import numpy as np

lons = np.array([13.404954, 13.388860, 13.428555])
lats = np.array([52.520008, 52.517037, 152.523430])   # row 2 is corrupt

easting, northing = fwd.transform(lons, lats)

# 1. Isolate the rows PROJ could not compute instead of trusting a spot check
failed = ~(np.isfinite(easting) & np.isfinite(northing))
print("failed row indices:", np.flatnonzero(failed))
# failed row indices: [2]

# 2. Every surviving easting must sit inside the UTM zone's defined band
ok = ~failed
assert np.all((easting[ok] > 166_000) & (easting[ok] < 834_000)), "easting outside zone 33"

# 3. Northings in the northern hemisphere run 0 to ~9,330,000 m
assert np.all((northing[ok] >= 0) & (northing[ok] < 9_400_000)), "northing outside hemisphere"
print(f"{ok.sum()} of {len(lons)} rows transformed cleanly")
# 2 of 3 rows transformed cleanly

Run those three checks as a fixture in CI rather than as a one-off. They cost microseconds on a million-row array and they catch the two failures that unit tests on a single hard-coded point cannot: a partially corrupt input batch, and a target zone that was correct for last quarter's extract but not for this one.

Triage table from symptom to root cause to fix Four rows, each pairing a symptom with its cause and its one-line fix. A CRSError on a valid EPSG code means a mismatched PROJ database from mixing pip and conda, fixed by reinstalling from conda-forge. A silent result about four thousand kilometres off means the lat-first axis order of EPSG 4326, fixed with always_xy equals True. Infinity or not-a-number output means the point falls outside the CRS area of use, fixed by bounds-checking the input. Coordinates off by about a hundred metres mean a bare PROJ string with no datum shift, fixed by rebuilding the CRS from EPSG or WKT2. From symptom to fix: the four PyProj failures SYMPTOM ROOT CAUSE ONE-LINE FIX CRSError: Invalid projection raised on a valid EPSG code PROJ database mismatch pip wheel over a conda GDAL reinstall from conda-forge check pyproj.proj_version_str no exception, ~4000 km off the numbers still look plausible EPSG:4326 declares lat first lon fed where lat is expected always_xy=True on every from_crs() call inf / nan in the output only some rows affected outside the area of use point sits beyond the zone bounds-check lon and lat then confirm the target zone off by roughly 100 m round-trip drifts by metres +proj= with no datum shift PROJ.4 text names no grids rebuild from EPSG or WKT2 CRS.from_epsg(25833)
Each failure has one signature and one fix: validate the CRS pair, force the axis order, bounds-check the input, and never trust a bare PROJ string.

Edge Cases & Debugging

Frequently Asked Questions

Should always_xy=True just be the default in every project? In practice, yes — set it on every Transformer.from_crs() call and treat a transformer built without it as a code-review failure. The reason is that essentially every data format you will actually meet stores coordinates x-first: GeoJSON mandates longitude-then-latitude, shapefiles and GeoPackages store x-then-y, and both Shapely and GeoPandas expose .x and .y in that order. The authority-declared latitude-first ordering for EPSG:4326 is correct as metadata and almost never matches the bytes on disk. The narrow exception is when you are deliberately consuming a service that honours the registry order — some WFS 1.1.0 and WMS 1.3.0 endpoints do — in which case leave the flag off, but comment the line so the next reader knows it was a decision rather than an omission.

Why does my round-trip close perfectly when the coordinates are still wrong? Because a round-trip only tests that the forward and inverse operations are mutual inverses, and a wrong operation is still its own inverse. If PROJ selected a ballpark datum shift, the same ballpark runs in reverse and its error cancels exactly, so the check returns to the input to within floating-point noise. The same is true of an axis swap applied consistently in both directions. A round-trip proves internal consistency; only an external reference proves accuracy. Add at least one known control point — a published survey mark, or a coordinate you have verified in another tool — and assert against its true value, not against your own input.

Is it safe to share one Transformer across threads or processes? Across threads, yes on pyproj 3.7 and later: the library maintains one PROJ context per thread automatically, which is why pyproj.set_use_global_context() is now deprecated. Across processes, no. A transformer created before a fork carries a context that is not valid in the child, and the failure mode ranges from a segfault to silently wrong output rather than a clean exception. The safe pattern is a module-level lru_cached factory that workers call after starting, so each process compiles its own pipeline once. On Windows and macOS, where spawn is the default start method, the problem does not arise because nothing is inherited — which is exactly why this bug tends to appear only when a pipeline moves to a Linux host.

My code broke after upgrading from pyproj 2.x — what actually changed? Three removals account for almost all of it. pyproj.transform(p1, p2, x, y) and the Proj-object calling convention are gone, replaced by Transformer. +init=epsg:4326 style strings are rejected; use EPSG:4326 or CRS.from_epsg(4326). And PROJ 6 introduced authority axis order, so code that previously got away with passing (lon, lat) now needs always_xy=True to keep doing so. The upgrade is worth it — the Transformer object is what makes caching and vectorized array transforms possible — but do it with a regression fixture of known coordinate pairs rather than by inspection, because the axis change is the one that produces no error at all.

How do I catch an axis swap in a test instead of in production? Assert on the extent rather than on individual points. After building geometry, compare total_bounds against the bounding box you expect the dataset to occupy; a swapped pair puts longitude values where latitudes belong and the bounds immediately exceed ±90 in the y slot, which is impossible for real latitudes. For projected output, the UTM easting band check shown under Verification does the same job in metres. Both are one-line assertions that run in microseconds and fail loudly on the exact class of bug that produces no exception, which makes them far more valuable in CI than a broader but slower geometry validity sweep.

Does Transformer.from_crs() hit the network? Only if grid downloads are enabled and the chosen operation needs a grid that is not on disk. Building the transformer itself is a local proj.db lookup; the network fetch, when it happens, is triggered lazily by the first coordinate and pulls only the byte ranges of the GeoTIFF grid that cover your extent. In a container that means a first-call latency spike and a hard dependency on the PROJ CDN being reachable. If either is unacceptable, bake the grids you need into the image and leave PROJ_NETWORK off, then assert the operation you expected was actually used before serving traffic.

Related