Reprojecting GeoJSON Between EPSG:4326 and EPSG:3857 with pyproj

TL;DR: Build one pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) and reuse it, then run each feature’s geometry through shapely.ops.transform(transformer.transform, geom) to convert every coordinate from geographic degrees to Web Mercator metres. The always_xy=True flag is non-negotiable: GeoJSON is always longitude-first per RFC 7946, but pyproj otherwise honours EPSG:4326’s latitude-first authority axis order and will silently swap your coordinates.

The Core Challenge: Axis Order and Two Very Different Units

Reprojecting a FeatureCollection looks trivial until two subtleties collide. The first is axis order. GeoJSON fixes coordinate order as [longitude, latitude] for every geometry, in every CRS, by specification — RFC 7946 leaves no room for interpretation. pyproj, by contrast, respects the axis order declared by the EPSG authority, and EPSG:4326 officially lists latitude as axis 1 and longitude as axis 2. Feed a GeoJSON coordinate straight into a default Transformer and it reads your longitude as a latitude. The output is not an error — it is a plausible-looking point in the wrong hemisphere. This is the same class of defect that plagues OGC request handling, covered from the service-boundary angle in SRS and Coordinate Reference System Handling and diagnosed in detail in the Handling Spatial Reference Mismatches in OGC Requests guide.

The second subtlety is units and domain. EPSG:4326 is geographic: coordinates are decimal degrees on the WGS 84 ellipsoid. EPSG:3857 — “Web Mercator”, the projection behind virtually every slippy map tile — is planar: coordinates are metres east and north of the origin, spanning roughly ±20,037,508 metres. Web Mercator is only defined between about 85.06°N and 85.06°S; push a latitude to the poles and the northing runs to infinity. Any reprojection routine that ignores this clipping boundary will emit inf values that corrupt the whole file.

GeoJSON Reprojection Data Flow A GeoJSON FeatureCollection whose coordinates are longitude and latitude in degrees flows into a cached pyproj Transformer built with always_xy set to true, which is applied per geometry by shapely.ops.transform, producing a new FeatureCollection whose coordinates are easting and northing in metres in EPSG:3857. EPSG:4326 FeatureCollection [lon, lat] degrees cached Transformer from_crs(4326, 3857) always_xy = True shapely.ops.transform per geometry EPSG:3857 FeatureCollection [x, y] metres RFC 7946: GeoJSON is always lon, lat — always_xy keeps pyproj in the same order

Production-Ready Code

The routine below reprojects a whole FeatureCollection in either direction. It caches the Transformer per CRS pair, loads each geometry with shapely.geometry.shape (which handles Point, LineString, Polygon, and their Multi* and GeometryCollection variants without any per-type branching), applies the transform with shapely.ops.transform, and serialises the result back to GeoJSON with mapping. The only third-party dependencies are pyproj and shapely.

"""
reproject_geojson.py — reproject a GeoJSON FeatureCollection with pyproj + shapely

Dependencies (pip install):
    pyproj>=3.6, shapely>=2.0
"""
from __future__ import annotations

import json
from functools import lru_cache
from typing import Any

from pyproj import Transformer
from shapely.geometry import mapping, shape
from shapely.ops import transform as shapely_transform

# Web Mercator (EPSG:3857) is only defined to ~85.06 degrees of latitude.
WEB_MERCATOR_MAX_LAT = 85.06


@lru_cache(maxsize=16)
def get_transformer(src_crs: str, dst_crs: str) -> Transformer:
    """
    Build and cache a pyproj Transformer for a source -> target CRS pair.

    always_xy=True forces (longitude, latitude) / (easting, northing) input
    and output order, matching GeoJSON's RFC 7946 lon-first convention and
    overriding EPSG:4326's authority (latitude-first) axis order.
    """
    return Transformer.from_crs(src_crs, dst_crs, always_xy=True)


def _round_coords(obj: Any, ndigits: int) -> Any:
    """Recursively round every numeric coordinate in a GeoJSON geometry dict."""
    if isinstance(obj, float):
        return round(obj, ndigits)
    if isinstance(obj, (list, tuple)):
        return [_round_coords(item, ndigits) for item in obj]
    return obj


def _clamp_to_mercator(lon: float, lat: float) -> tuple[float, float]:
    """Clamp a lon/lat pair into Web Mercator's valid latitude domain."""
    clamped_lat = max(-WEB_MERCATOR_MAX_LAT, min(WEB_MERCATOR_MAX_LAT, lat))
    return lon, clamped_lat


def reproject_feature_collection(
    collection: dict[str, Any],
    src_crs: str = "EPSG:4326",
    dst_crs: str = "EPSG:3857",
    precision: int | None = None,
    clamp_poles: bool = True,
) -> dict[str, Any]:
    """
    Return a new FeatureCollection with every geometry reprojected to dst_crs.

    Features are copied, not mutated. Features with a null geometry pass
    through untouched. When reprojecting geographic degrees into EPSG:3857,
    clamp_poles guards against infinite northings beyond +/-85.06 degrees.
    """
    transformer = get_transformer(src_crs, dst_crs)
    project_to_mercator = src_crs == "EPSG:4326" and dst_crs == "EPSG:3857"

    def _project(x: Any, y: Any, z: Any = None) -> tuple[Any, ...]:
        # shapely calls this with flat coordinate sequences; pyproj accepts
        # scalars or iterables and returns them in the same shape.
        if clamp_poles and project_to_mercator:
            x, y = _clamp_to_mercator(x, y)
        if z is None:
            return transformer.transform(x, y)
        return transformer.transform(x, y, z)

    out_features: list[dict[str, Any]] = []
    for feature in collection.get("features", []):
        geometry = feature.get("geometry")
        if geometry is None:
            out_features.append({**feature, "geometry": None})
            continue

        reprojected = shapely_transform(_project, shape(geometry))
        new_geometry = mapping(reprojected)
        if precision is not None:
            new_geometry["coordinates"] = _round_coords(
                new_geometry["coordinates"], precision
            )
        out_features.append({**feature, "geometry": new_geometry})

    return {"type": "FeatureCollection", "features": out_features}


if __name__ == "__main__":
    sample = {
        "type": "FeatureCollection",
        "features": [
            {
                "type": "Feature",
                "properties": {"name": "Greenwich Observatory"},
                "geometry": {"type": "Point", "coordinates": [-0.0014, 51.4779]},
            },
            {
                "type": "Feature",
                "properties": {"name": "Thames reach"},
                "geometry": {
                    "type": "LineString",
                    "coordinates": [[-0.02, 51.50], [0.00, 51.51], [0.02, 51.49]],
                },
            },
        ],
    }

    projected = reproject_feature_collection(sample, precision=2)
    print(json.dumps(projected, indent=2))

    # Round trip back to EPSG:4326 to confirm the transform is reversible.
    restored = reproject_feature_collection(
        projected, src_crs="EPSG:3857", dst_crs="EPSG:4326", precision=6
    )
    print(json.dumps(restored["features"][0]["geometry"], indent=2))

Step-by-Step Walkthrough

Building and caching the transformer. get_transformer wraps Transformer.from_crs behind functools.lru_cache. The first call for a given (src_crs, dst_crs) pair does the expensive work — resolving both CRS definitions and compiling a PROJ transformation pipeline — and every subsequent call for the same pair returns the identical, thread-safe Transformer instance. On a collection of thousands of features this turns thousands of pipeline builds into one. The cache is keyed on the CRS strings, so reprojecting in both directions (4326→3857 and 3857→4326) stores two distinct transformers.

The always_xy=True flag. This is the single most important argument in the file. Without it, Transformer.from_crs("EPSG:4326", "EPSG:3857") expects input in EPSG:4326’s authority order — latitude first. Your GeoJSON supplies longitude first, so the transformer would treat -0.0014 as a latitude and 51.4779 as a longitude, projecting a point near Greenwich to somewhere off the coast of Somalia. always_xy=True tells pyproj to use the “traditional GIS” order, (x, y) / (longitude, latitude), on both input and output. That order is exactly what GeoJSON guarantees under RFC 7946, so no manual axis swapping is ever needed.

Applying the transform per geometry. shape(geometry) parses a GeoJSON geometry dict into the corresponding shapely object. shapely.ops.transform then walks that object’s coordinate sequences, calling the supplied _project function and rebuilding a geometry of the same type with the returned coordinates. Because shapely handles the traversal, one code path covers Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon, and GeometryCollection — ring structure, holes, and part counts are all preserved automatically. The _project callable accepts an optional z so that geometries carrying elevation survive the round trip.

Pole clamping. When the source is EPSG:4326 and the target is EPSG:3857, _project first passes each coordinate through _clamp_to_mercator, which pins latitude to the ±85.06° band. Without this guard, a feature touching 90°N would project to an infinite northing and poison the entire output document. Clamping is disabled automatically for any other CRS pair and can be turned off with clamp_poles=False when you would rather drop out-of-domain features yourself.

Rewriting the FeatureCollection. Each reprojected geometry is serialised back to a plain dict with mapping, and features are rebuilt with {**feature, "geometry": new_geometry} so properties and IDs are preserved and the original input is never mutated. The optional precision argument rounds coordinates recursively — a small step that removes noisy floating-point digits and can shrink a serialised file by a third with no visible loss of accuracy.

Verification

Run the script directly. It reprojects the sample collection to EPSG:3857, then reverses the transform to confirm the original degrees are recovered.

python reproject_geojson.py

Expected output (truncated):

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": { "name": "Greenwich Observatory" },
      "geometry": {
        "type": "Point",
        "coordinates": [ -155.84, 6708225.7 ]
      }
    }
  ]
}

The recovered geometry after the round trip prints as:

{
  "type": "Point",
  "coordinates": [ -0.0014, 51.4779 ]
}

Confirm two things: the EPSG:3857 easting for a point near 0° longitude is close to zero metres, and the round trip returns the original longitude and latitude to six decimal places. A quick numeric assertion belongs in your test suite:

restored = reproject_feature_collection(
    reproject_feature_collection(sample),
    src_crs="EPSG:3857",
    dst_crs="EPSG:4326",
)
lon, lat = restored["features"][0]["geometry"]["coordinates"]
assert abs(lon - (-0.0014)) < 1e-6 and abs(lat - 51.4779) < 1e-6

Gotchas & Edge Cases

Why is always_xy=True required when reprojecting GeoJSON?

GeoJSON coordinates are always longitude-first, latitude-second per RFC 7946, whatever the CRS. pyproj instead honours the authority axis order of EPSG:4326, which is latitude-first. Passing always_xy=True forces the Transformer to accept and emit coordinates in (longitude, latitude) / (x, y) order, matching GeoJSON exactly. Omit it and the output is silently swapped — no exception is raised, so the bug only surfaces when a point lands in the wrong hemisphere. The same axis-order trap appears when building OGC request parameters, as the SRS and Coordinate Reference System Handling guide explains.

Why do points near the poles produce infinite EPSG:3857 coordinates?

Web Mercator (EPSG:3857) is only defined between roughly 85.06°N and 85.06°S. As latitude approaches ±90° the Mercator projection tends to infinity, so pyproj returns inf for the northing. Clamp latitudes into the valid band before reprojecting — as _clamp_to_mercator does above — or filter out features that fall outside the projection’s domain. Never serialise an inf into GeoJSON: it is not valid JSON and most parsers reject the whole document.

Should I build a new Transformer for every feature?

No. Transformer.from_crs performs a CRS lookup and compiles a PROJ pipeline on first construction, which is expensive relative to a single coordinate transform. Build the Transformer once per source-to-target CRS pair and reuse it across every feature — here via functools.lru_cache. Reusing a single instance across thousands of geometries is dramatically faster and, because pyproj Transformer objects are immutable and thread-safe, safe to share across worker threads.

How much precision should reprojected coordinates keep?

EPSG:3857 coordinates are metres, so one decimal place already gives centimetre precision — more than enough for web mapping. In EPSG:4326, six decimal places of longitude is about 0.1 metres at the equator. Rounding EPSG:3857 output to one or two decimals and EPSG:4326 output to six or seven decimals via the precision argument strips meaningless floating-point noise and shrinks the serialised file without any practical loss of accuracy.

Do I need to update the GeoJSON crs member after reprojecting?

RFC 7946 removed the named crs member and mandates that GeoJSON coordinates be WGS 84 longitude/latitude — that is, EPSG:4326. Strictly speaking, a FeatureCollection in EPSG:3857 is no longer conformant GeoJSON, so a compliant consumer will assume degrees. When you reproject to Web Mercator for an internal pipeline or a tiling step, keep the target CRS out of band (a filename convention, an HTTP header, or a wrapper object) rather than relying on a deprecated in-document crs member, and reproject back to EPSG:4326 before publishing the file to any standards-based client such as a WMS or WFS.

Back to SRS and Coordinate Reference System Handling

Related