Converting OSM coordinates to local CRS with pyproj

Take decoded OpenStreetMap node arrays stored in implicit WGS 84 (EPSG:4326) and reproject them into a local projected CRS — UTM, LAEA, or Web Mercator — using pyproj’s Transformer API without silently swapping the latitude and longitude axes.

Prerequisites

Confirm each item before running the code below; an unmet prerequisite is the usual cause of a “works on my laptop, NaNs in production” reprojection bug.

Conceptual minimum

OpenStreetMap persists every node, way, and relation in unprojected WGS 84 geographic coordinates, and — as Coordinate Reference Systems in OSM explains — that CRS is implicit: no projection string is stored against any primitive. Angular degrees are the wrong unit for buffering, distance, and topology work, so an analytics pipeline must reproject into a metric Cartesian system. Because only nodes carry coordinates in the Node-Way-Relation data model, you transform raw node arrays first and assemble way and relation geometries afterward.

The single rule that prevents most corruption is axis order. PROJ follows the EPSG registry, which defines EPSG:4326 as latitude-first, while OSM tooling, GeoJSON, and Shapely all expect (longitude, latitude). Passing always_xy=True forces pyproj to treat the X argument as longitude and Y as latitude regardless of the CRS pair or PROJ version, removing a brittle implicit dependency. For local accuracy, pick the UTM zone covering the extract centroid, where the zone number follows:

Cost of transformer construction against transformer reuse A horizontal bar chart of wall-clock seconds to reproject one million nodes. Constructing a Transformer inside the loop takes 612 seconds. Constructing it once and calling it per point takes 44 seconds. Constructing it once and passing the whole coordinate array in a single call takes 1.1 second. Where the time goes: constructing the transformer, not using it 1 M nodes reprojected, wall-clock seconds Transformer per point 612 s Built once, called per point 44 s Built once, whole array at once 1.1 s Construction parses the CRS definition and builds a pipeline; the projection maths itself is a handful of floating-point operations.
Two independent wins stack here: hoisting construction out of the loop, then handing pyproj whole arrays so the per-call overhead is amortised across a million points instead of paid a million times.

zone=⌊λ+1806⌋+1\text{zone} = \left\lfloor \frac{\lambda + 180}{6} \right\rfloor + 1

with λ\lambda the centroid longitude in decimal degrees.

Axis-order pitfall: always_xy=True versus the default Two parallel pipelines start from one OSM node with latitude 52.52 and longitude 13.40. The top lane passes always_xy=True, so the pyproj Transformer maps X to longitude and Y to latitude, yielding valid UTM Zone 33N coordinates x≈392,440 and y≈5,820,080 over Berlin. The bottom lane omits always_xy, so PROJ follows EPSG authority order and reads the first argument as latitude; the axes are swapped and the point projects outside the zone extent, landing in the wrong place. OSM node lat = 52.52 lon = 13.40 always_xy=True → X is longitude Transformer in: (13.40, 52.52) UTM 33N (x, y) 392440, 5820080 ✓ Valid — lands on Berlin inside zone extent default → X read as latitude Transformer reads lat = 13.40 axes swapped lon 52.52 → off-zone ✗ Out of bounds ocean / inf / nan

Runnable solution

The snippet builds a cached Transformer, then streams (lat, lon) tuples through it in memory-bounded NumPy chunks, yielding projected (x, y) arrays in the target CRS.

python 
import logging
from typing import Iterable, Iterator, Tuple

import numpy as np
from pyproj import CRS, Transformer

logger = logging.getLogger("osm.reproject")

SOURCE_CRS = CRS.from_epsg(4326)            # OSM's implicit WGS 84
TARGET_CRS = CRS.from_epsg(32633)           # UTM Zone 33N — set per study area

# Build ONE transformer per process. Initialization queries the PROJ
# operation database and loads any datum-shift grids, so never rebuild
# it inside a loop or per worker task.
TRANSFORMER = Transformer.from_crs(
    SOURCE_CRS,
    TARGET_CRS,
    always_xy=True,                         # map X<-lon, Y<-lat explicitly
)


def chunk_transform(
    lat_lon_iter: Iterable[Tuple[float, float]],
    transformer: Transformer = TRANSFORMER,
    chunk_size: int = 1_000_000,
) -> Iterator[np.ndarray]:
    """Yield projected (N, 2) float64 arrays in the target CRS.

    Input tuples are (lat, lon) as OSM stores them; they are reordered to
    (lon, lat) before transforming. Memory per chunk is ~16 MB at 1M points.
    """
    buffer: list[Tuple[float, float]] = []

    def _flush(rows: list[Tuple[float, float]]) -> np.ndarray:
        arr = np.asarray(rows, dtype=np.float64)        # columns: lon, lat
        x, y = transformer.transform(arr[:, 0], arr[:, 1])
        out = np.column_stack((x, y))
        finite = np.isfinite(out).all(axis=1)
        if not finite.all():
            logger.warning(
                "dropped %d of %d points outside target CRS extent",
                int((~finite).sum()), out.shape[0],
            )
        return out[finite]

    for lat, lon in lat_lon_iter:
        buffer.append((lon, lat))                        # enforce (x=lon, y=lat)
        if len(buffer) >= chunk_size:
            yield _flush(buffer)
            buffer.clear()
    if buffer:
        yield _flush(buffer)


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    nodes = [(52.5200, 13.4050), (48.1372, 11.5755)]     # Berlin, Munich (lat, lon)
    for projected in chunk_transform(nodes):
        logger.info("projected chunk shape=%s sample=%s", projected.shape, projected[0])

Step-by-step walkthrough

  1. CRS construction — CRS.from_epsg(4326) and CRS.from_epsg(32633) resolve full WKT definitions from the EPSG registry. Using from_epsg rather than a raw proj-string guarantees the datum and ellipsoid are unambiguous.
  2. Transformer caching — Transformer.from_crs(...) is assigned to a module-level TRANSFORMER. The constructor performs a database lookup and may load grid-shift files, so it is built once and reused across calls and worker threads.
  3. always_xy=True — this pins argument order to (longitude, latitude), the convention every downstream tool expects, so axis order can never silently flip between PROJ versions.
  4. Axis reordering — each incoming (lat, lon) tuple is appended as (lon, lat), matching the X/Y contract enforced above.
  5. Vectorized flush — _flush builds a float64 array and calls transformer.transform on whole columns, letting PROJ’s C routines run without Python per-point overhead.
  6. Finite masking — np.isfinite(...).all(axis=1) drops any row that transformed to inf/nan (input outside the target CRS extent), and the count is logged rather than silently swallowed.
  7. Bounded streaming — chunk_size caps live memory; at 1,000,000 points a chunk holds roughly 16 MB of float64, leaving headroom on a standard worker even when several stages run concurrently. Pair this with memory-efficient chunk processing to keep the whole ingestion deterministic on multi-gigabyte extracts.
The three stages a pyproj Transformer resolves at construction time A left-to-right chain: source CRS EPSG:4326, then a datum shift from WGS 84 to the target datum, then the projection itself such as transverse Mercator, Lambert azimuthal equal-area or Mercator, ending in projected x and y in metres. A panel below notes that all three stages are resolved once at construction and cached, making the object expensive to build, cheap to call, and something to build per worker after a fork. What Transformer.from_crs actually assembles source CRS EPSG:4326 datum shift WGS 84 → target datum projection tmerc · laea · merc projected (x, y) metres All three stages are resolved once, at construction — and cached on the object A Transformer is therefore expensive to build, cheap to call, and safe to reuse for every node in the extract. It is not thread-safe to share across processes, so build one per worker after the fork, not before.
Reuse is safe because the object is immutable once built — but build it inside each worker process. A Transformer created before a fork carries a PROJ context that does not survive the copy cleanly.

Verification

Confirm the reprojection is correct before wiring it into the next stage:

  • Range check. For UTM Zone 33N, easting (x) should sit near 166,000–834,000 m and northing (y) be positive in the northern hemisphere. Berlin (52.52, 13.405) projects to roughly x ≈ 392,440, y ≈ 5,820,080.
  • Round-trip residual. Build the inverse transformer (Transformer.from_crs(TARGET_CRS, SOURCE_CRS, always_xy=True)), reproject the output back, and assert the residual is below your tolerance (< 1e-6 degrees for a clean grid path).
  • Log lines. A healthy run logs projected chunk shape=... and emits no dropped N of M points warnings; any drop warning means inputs fell outside the target CRS extent.
  • Sample audit. Compare a 1% random sample against known control points to confirm sub-meter agreement before trusting downstream metric joins.

Common errors and fixes

Symptom Root cause One-line fix
Coordinates land in the ocean / wrong hemisphere Axis order swapped (lat passed as X) Set always_xy=True and feed (lon, lat).
Output columns are all inf Point outside the target CRS valid extent (e.g. wrong UTM zone) Select the UTM zone for the extract centroid, or use EPSG:3035/3857.
CRSError: Invalid projection Source CRS never assigned to raw OSM input Construct with CRS.from_epsg(4326) at the ingestion boundary.
Throughput collapses on large extracts Per-node transform() calls in a Python loop Batch into NumPy arrays and transform whole columns.
Sub-meter drift vs. control points Datum-shift grid missing; PROJ fell back Provision grids in PROJ_DATA; set PROJ_NETWORK=OFF to fail loudly.
DeprecationWarning on Proj/transform Legacy pyproj 1.x API Migrate to the Transformer API shown above.

Specification reference

OpenStreetMap stores all geometry in WGS 84 (EPSG:4326); the datum is fixed by convention and is not encoded in the data. See the OSM Wiki on Node coordinates and the EPSG:4326 and EPSG:32633 definitions for axis order and valid extents. In PBF, raw integers are reconstructed via granularity and lat_offset/lon_offset before any reprojection — the PBF File Structure Deep Dive covers that decode step.

Projected node arrays from this procedure feed directly into the metric stages that follow — most often spatial indexing for OSM extracts, where R-tree, H3, or Quadkey structures accelerate proximity queries and boundary clipping.

Frequently Asked Questions

Is a Transformer safe to share between threads?

Reading from one concurrently is safe in current pyproj releases, because the underlying PROJ context is thread-local and the transformation itself does not mutate the object. Sharing across a fork, however, is not: the child inherits a context that was created in the parent, and behaviour ranges from a silent slow path to a crash. Build the transformer inside each worker after the fork, which is cheap enough when it happens once per worker and catastrophic when it happens once per point.

Why does pyproj sometimes download grid files?

Because an accurate datum transformation between some pairs of systems needs a correction grid that is too large to ship with the package. When the grid is missing, PROJ falls back to a less accurate transformation rather than failing, so results differ by metres between a machine with the grid and one without. If reproducibility across environments matters, either pre-fetch the grids with pyproj sync and bake them into the image, or pin the transformation explicitly so no fallback is possible.

Should I reproject before or after filtering?

After. Filtering is usually cheaper in the source system — a bounding-box test in degrees is a comparison, and the same test after projection needs the projection first — and reprojecting fewer points is strictly less work. The exception is a filter expressed in metres, such as “within 500 m of a line”, which cannot be evaluated correctly in degrees at all.

Up one level: Coordinate Reference Systems in OSM.