Transforming Point Coordinates with pyproj

To convert coordinates between two CRSs, build one Transformer with always_xy=True and reuse it:

from pyproj import Transformer

# always_xy=True → inputs and outputs are (x, y) = (lon, lat) or (easting, northing)
to_utm = Transformer.from_crs("EPSG:4326", "EPSG:32636", always_xy=True)

easting, northing = to_utm.transform(36.82, -1.29)      # lon, lat
print(round(easting, 2), round(northing, 2))

Coordinate transformation is the point-level counterpart of the raster warping covered in Mastering CRS Transformations in rasterio.

Why This Arises in Remote Sensing Workflows

Points cross CRS boundaries constantly. A bounding box arrives in degrees and has to become a pixel window in UTM metres. Field plots recorded on a handheld GPS have to be sampled against a projected scene. A STAC search takes geographic coordinates while the imagery it returns is projected. Tile boundaries in Web Mercator have to be mapped back onto native grids.

Two properties make this error-prone. The first is axis order: the EPSG registry defines EPSG:4326 with latitude first, while almost every software convention, GeoJSON included, writes longitude first. A transformer built without always_xy follows the authority, so passing longitude first produces a coordinate that is not merely wrong but confidently wrong — often somewhere in the ocean.

The second is that transformation is not a fixed function. Between two datums there may be several defined operations with different accuracies, and PROJ picks the best available given the grids installed on the machine. That is usually what you want, and it means the same code can produce slightly different numbers on two machines — a difference that matters when centimetres are being claimed.

The axis-order trap, drawn The pair 36.82 and minus 1.29 means longitude 36.82 east and latitude 1.29 south — a point in Kenya. Interpreted latitude-first it becomes latitude 36.82 north and longitude 1.29 west, a point in Spain. Both transform without error, and only one is right. The pair (36.82, −1.29) — two readings, two continents always_xy=True x = 36.82° E, y = 1.29° S → EPSG:32636 E 258,000 · N 9,857,000 inside the scene, as intended authority order (default) lat = 36.82° N, lon = 1.29° W → EPSG:32636 far outside zone 36 no error raised — just a wrong place A magnitude assertion catches this instantly: zone-36 eastings live between 166,021 and 833,978.

Environment & Setup

Package Version Why
pyproj ≥3.4 Transformer, CRS, area-of-interest selection
numpy ≥1.23 Vectorised transformation of coordinate arrays
rasterio ≥1.3.0 Supplies the raster CRS the points are transformed into 
pip install "pyproj>=3.4" "numpy>=1.23" "rasterio>=1.3.0"

Complete Working Example

This helper builds one transformer per CRS pair, caches it, transforms arrays in bulk, and asserts that the results are plausible for the target CRS.

Cost of building versus using a Transformer Constructing a Transformer resolves a coordinate operation from the PROJ database and may consult grid files, costing milliseconds. Transforming an array of coordinates costs microseconds. Building one per point inverts that ratio and dominates the runtime of any large job. Build once, transform in bulk pattern time for 100,000 points why Transformer per point ~46 s 100,000 database lookups one Transformer, point loop ~1.9 s Python call overhead per point one Transformer, arrays ~0.04 s a single vectorised call cached per CRS pair ~0.04 s same, and safe across call sites The last row is what the lru_cache in the example buys: correctness by default at call sites. 
from functools import lru_cache

import numpy as np
from pyproj import CRS, Transformer


@lru_cache(maxsize=32)
def transformer(src_crs: str, dst_crs: str) -> Transformer:
    """One Transformer per CRS pair, built once and reused."""
    return Transformer.from_crs(src_crs, dst_crs, always_xy=True)


def to_crs(xs, ys, src_crs: str, dst_crs: str) -> tuple[np.ndarray, np.ndarray]:
    """Transform coordinate arrays, returning float64 arrays in the target CRS."""
    tf = transformer(src_crs, dst_crs)
    out_x, out_y = tf.transform(np.asarray(xs, dtype="float64"),
                                np.asarray(ys, dtype="float64"))
    out_x, out_y = np.asarray(out_x), np.asarray(out_y)

    if not np.isfinite(out_x).all() or not np.isfinite(out_y).all():
        n_bad = int((~np.isfinite(out_x) | ~np.isfinite(out_y)).sum())
        raise ValueError(f"{n_bad} coordinate(s) fell outside the transformation's domain")
    return out_x, out_y


def assert_plausible(xs: np.ndarray, ys: np.ndarray, crs: str) -> None:
    """Cheap sanity check that the numbers belong to the CRS they claim."""
    c = CRS.from_user_input(crs)
    if c.is_geographic:
        assert np.all(np.abs(xs) <= 180) and np.all(np.abs(ys) <= 90), "not degrees"
    elif "UTM" in (c.name or "").upper():
        assert np.all((xs > 100_000) & (xs < 900_000)), "eastings outside the UTM band"
        assert np.all((ys >= 0) & (ys <= 10_000_000)), "northings outside the UTM range"


if __name__ == "__main__":
    lon = np.array([36.82, 36.90, 37.01])
    lat = np.array([-1.29, -1.35, -1.41])

    east, north = to_crs(lon, lat, "EPSG:4326", "EPSG:32636")
    assert_plausible(east, north, "EPSG:32636")

    back_lon, back_lat = to_crs(east, north, "EPSG:32636", "EPSG:4326")
    print("round-trip error (m-scale degrees):",
          float(np.max(np.abs(back_lon - lon))), float(np.max(np.abs(back_lat - lat))))

The round trip is the strongest cheap test available: a transformation that is wrong in either direction rarely returns to its starting point, and the residual should be at the level of floating-point noise rather than metres.

Variant Patterns

1. Transforming a bounding box, not a pair of corners

A projected bounding box is not the transform of its two corners: the edges curve, so the true extent is larger than the corner-to-corner rectangle. Use the dedicated helper, which densifies the edges.

from pyproj import Transformer

tf = Transformer.from_crs("EPSG:4326", "EPSG:32636", always_xy=True)

# Wrong: corners only — under-covers the true footprint
naive = tf.transform([36.5, 37.2], [-1.6, -1.0])

# Right: densified edges
minx, miny, maxx, maxy = tf.transform_bounds(36.5, -1.6, 37.2, -1.0, densify_pts=21)

The difference matters most for wide extents and high latitudes, where a naive box can miss several kilometres of ground — and it is the same reasoning behind calculate_default_transform densifying its own footprint, described in Reprojecting a Raster from UTM to WGS84.

2. Pinning the transformation for reproducibility

Why two machines can disagree by a metre PROJ chooses between available coordinate operations. A grid-based operation gives sub-decimetre accuracy but requires the grid file to be installed. A Helmert approximation is always available and accurate to about a metre. The same code picks different paths on machines with different grid availability. Same source and target CRS, two operations source CRS EPSG:4326 grid-based operation accuracy ≈ 0.05 m · needs the grid file Helmert approximation accuracy ≈ 1 m · always available target CRS EPSG:32636 List the candidates with TransformerGroup and pin one when reproducibility matters more than accuracy. 
from pyproj.transformer import TransformerGroup

group = TransformerGroup("EPSG:4326", "EPSG:32636", always_xy=True)
for t in group.transformers[:3]:
    print(t.description, "| accuracy:", t.accuracy)

# Pin the first available operation so every machine produces identical numbers
pinned = group.transformers[0]

For most remote sensing work the metre-level difference is irrelevant next to a 10 m pixel. For survey-grade or change-detection work at sub-pixel scale, it is not, and pinning is the difference between reproducible and approximately reproducible.

3. Points into pixel indices

The common end goal is not a coordinate but a pixel. Combine the transformer with the raster’s inverse transform in one step.

import rasterio
from pyproj import Transformer

with rasterio.open("scene_utm.tif") as src:
    tf = Transformer.from_crs("EPSG:4326", src.crs, always_xy=True)
    xs, ys = tf.transform(lon, lat)
    rows, cols = src.index(xs, ys)          # map coords → (row, col)
    inside = [(0 <= r < src.height) and (0 <= c < src.width) for r, c in zip(rows, cols)]

Checking inside explicitly matters because index happily returns out-of-range indices for points outside the scene, and reading with them either raises or silently wraps depending on how the array is sliced — the same trap covered in Extracting Pixel Values at Point Locations.

Performance Notes

Transformation itself is fast; the surrounding code usually is not. Three habits keep it that way.

Build the transformer outside the loop. Constructing one resolves the operation from the PROJ database and may consult grid files, which is orders of magnitude more expensive than transforming a point. The lru_cache in the example is there precisely because call sites forget.

Pass arrays, not points. transform accepts NumPy arrays and processes them in one call; transforming a million points one at a time spends its life in Python overhead rather than in PROJ.

Reuse a single transformer across threads with care. A Transformer is not guaranteed thread-safe in every PROJ build, so give each worker its own instance — the cache keyed per pair makes that a one-line change — rather than sharing one across a pool.

Common Errors

The transformed points are in the wrong hemisphere

Axis order. Add always_xy=True and confirm you are passing longitude first.

inf values come back for some points

Those coordinates fall outside the transformation’s domain — commonly a UTM zone that does not cover them. Check the points against the zone’s valid range, or transform to a CRS with global coverage.

The transformation is slow inside a loop over features

A Transformer is being constructed per feature. Hoist it out, or cache it as shown above.

Frequently Asked Questions

Q: What does always_xy actually change? It forces longitude-latitude and easting-northing ordering regardless of what the CRS authority declares. EPSG defines EPSG:4326 as latitude first, so without always_xy a transformer expects (lat, lon) and silently returns nonsense when handed (lon, lat).

Q: Is creating a Transformer expensive? Yes, relative to using one. It resolves the coordinate operation from the PROJ database, sometimes downloading a grid. Build it once outside any loop and reuse it; building one per point is the most common performance mistake here.

Q: Why do my transformed coordinates differ by a metre from another tool? Different datum transformation paths. PROJ may choose a grid-based operation when the grid is available and a coarser Helmert approximation when it is not, and the two differ by up to a few metres. Pin the operation explicitly when the difference matters.