Converting OSM coordinates to local CRS with pyproj Jump to heading#

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 Jump to heading#

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.

• pyproj ≥ 3.4.0 installed (pip install "pyproj>=3.4") — it bundles PROJ ≥ 9.0 and the Transformer operation database.
• numpy ≥ 1.23 installed for vectorized array transforms.
• A streaming PBF reader available (pip install osmium for pyosmium) so coordinates feed the buffer without per-object Python overhead.
• The target EPSG code chosen for your study area (e.g. 32633 for UTM Zone 33N over central Europe; 3035 for Europe-wide equal-area; 3857 for slippy-map tiles).
• Required datum-shift grids provisioned in PROJ_DATA; set PROJ_NETWORK=OFF in air-gapped runs so PROJ never silently fetches or falls back to a lower-accuracy path.
• Input coordinates already validated against the WGS 84 envelope (this is the output of the reprojection stage’s upstream bounds check).

Conceptual minimum Jump to heading#

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:

$$\text{zone} = \left\lfloor \frac{\lambda + 180}{6} \right\rfloor + 1$$

with $\lambda$ the centroid longitude in decimal degrees.

Runnable solution Jump to heading#

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.

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 Jump to heading#

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.

Verification Jump to heading#

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 Jump to heading#

Specification reference Jump to heading#

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.

Related Jump to heading#

• Coordinate Reference Systems in OSM — how OSM encodes coordinates and why the CRS is implicit.
• Spatial Indexing for OSM Extracts — index the projected (x, y) arrays this page produces.
• PBF File Structure Deep Dive — the delta/granularity decode that yields the coordinates you reproject.
• Node-Way-Relation Data Model — which primitives carry coordinates versus inherit them.
• Memory-Efficient Chunk Processing — bound memory while streaming large extracts through the transformer.
• OSM Data Fundamentals & Architecture — the foundation this stage sits within.

Up one level: Coordinate Reference Systems in OSM.