Reprojecting UTM-Zone Data to Web Mercator with pyproj

Part of the CRS & Projection Management guide.

Operative rule: never feed raw UTM easting/northing meters (EPSG:326xx / EPSG:327xx) to a web tile renderer — reproject to EPSG:3857 (or EPSG:4326 for GeoJSON) with geopandas .to_crs() first, or every feature lands hundreds of kilometres off-target.

How UTM-to-Web-Mercator Reprojection Works

The Universal Transverse Mercator system divides the globe into 60 six-degree longitude zones, each with its own projected coordinate system measured in meters. Survey data, government cadastral exports, and drone-derived rasters commonly arrive in a single UTM zone because it minimises local distortion — an ideal choice for area and distance calculations. Web maps, however, live in a different world: tile renderers like Leaflet and maplibregl draw on the spherical Web Mercator grid, EPSG:3857, and accept geographic degrees (EPSG:4326) on their public API. UTM meters and Web Mercator meters are both meter-based but are not interchangeable — an easting of 500000 in EPSG:32633 is a completely different point from 500000 in EPSG:3857.

Reprojection bridges the two. pyproj builds a Transformer from the PROJ database that maps each source coordinate through the source zone’s inverse projection back to the ellipsoid, then forward into the target projection. In practice you rarely call the transformer directly: geopandas wraps it in GeoDataFrame.to_crs(), which reprojects an entire geometry column in one vectorised pass. The one subtlety that bites everyone is axis order — some CRS definitions are latitude-first — which is why always_xy=True matters when you do reach for the raw transformer. This is the same axis-order trap discussed in implementing EPSG:3857 vs EPSG:4326 in Folium, and getting the CRS right upstream is what lets downstream tile vs vector rendering strategies align features to the basemap.

UTM to Web Mercator Reprojection Pipeline Stage one detects the UTM zone EPSG code from the data extent. Stage two reprojects with geopandas to_crs to EPSG:3857 for tiles or EPSG:4326 for GeoJSON. Stage three validates the resulting bounds fall in the legal Web Mercator range. Stage four hands clean coordinates to the web renderer. UTM source EPSG:326xx / 327xx easting / northing (m) to_crs() pyproj transform always_xy=True validate bounds within ±20037508 m (EPSG:3857 extent) web renderer Leaflet / MapLibre tiles aligned ⚠ Skipping stage 2 sends UTM meters to the renderer — features land off-map

Production-Ready Implementation

The module below detects the correct UTM EPSG code from a data sample, reprojects a GeoDataFrame to EPSG:3857, and validates that the output bounds fall inside the legal Web Mercator extent. It uses geopandas for the bulk transform and a raw pyproj.Transformer only for the single-point zone-detection helper, where always_xy=True is set explicitly.

from __future__ import annotations

import math
from pathlib import Path

import geopandas as gpd
import pyproj

# EPSG:3857 valid extent in meters (square Web Mercator world).
WEB_MERCATOR_LIMIT = 20_037_508.342789244


def utm_epsg_for_lonlat(lon: float, lat: float) -> int:
    """Return the EPSG code of the UTM zone containing (lon, lat).

    Northern hemisphere -> 326zz, southern -> 327zz, zz = zone number.
    """
    zone = int(math.floor((lon + 180) / 6) % 60) + 1
    return (32600 if lat >= 0 else 32700) + zone


def reproject_utm_to_web_mercator(
    gdf: gpd.GeoDataFrame,
    target_epsg: int = 3857,
) -> gpd.GeoDataFrame:
    """Reproject a UTM GeoDataFrame to Web Mercator (or 4326) and validate.

    Args:
        gdf: Source frame; must carry a defined CRS (a UTM zone).
        target_epsg: 3857 for tile rendering, 4326 for GeoJSON export.

    Returns:
        A reprojected copy with validated bounds.
    """
    if gdf.crs is None:
        raise ValueError(
            "Source CRS is undefined. Set it first, e.g. "
            "gdf.set_crs('EPSG:32633') — reprojecting without a known "
            "source zone silently produces garbage."
        )

    # Single vectorised pass over the geometry column. geopandas builds the
    # pyproj Transformer internally with correct axis handling.
    out = gdf.to_crs(epsg=target_epsg)

    # Validate: for 3857 the extent is a fixed square; for 4326 it is degrees.
    minx, miny, maxx, maxy = out.total_bounds
    if target_epsg == 3857:
        limit = WEB_MERCATOR_LIMIT
        ok = (-limit <= minx <= limit) and (-limit <= maxy <= limit)
    else:  # 4326 degrees
        ok = (-180 <= minx <= 180) and (-90 <= miny <= 90) and (maxy <= 90)

    if not ok:
        raise ValueError(
            f"Reprojected bounds {out.total_bounds} fall outside the "
            f"EPSG:{target_epsg} extent — the source CRS was likely wrong."
        )
    return out


def batch_reproject(src_dir: Path, dst_dir: Path, target_epsg: int = 3857) -> list[Path]:
    """Reproject every vector file in src_dir, writing results to dst_dir."""
    dst_dir.mkdir(parents=True, exist_ok=True)
    written: list[Path] = []
    for path in sorted(src_dir.glob("*.geojson")):
        gdf = gpd.read_file(path)

        # If the file lacks a CRS, infer the UTM zone from its centroid.
        if gdf.crs is None:
            c = gdf.to_crs(4326).geometry.union_all().centroid  # rough guess
            guessed = utm_epsg_for_lonlat(c.x, c.y)
            gdf = gdf.set_crs(epsg=guessed, allow_override=True)

        out = reproject_utm_to_web_mercator(gdf, target_epsg)
        out_path = dst_dir / path.name
        out.to_file(out_path, driver="GeoJSON")
        written.append(out_path)
    return written


if __name__ == "__main__":
    # Example: a survey file known to be UTM zone 33N (EPSG:32633).
    gdf = gpd.read_file("survey_utm33n.geojson").set_crs(
        "EPSG:32633", allow_override=True
    )
    web = reproject_utm_to_web_mercator(gdf, target_epsg=4326)  # for GeoJSON
    web.to_file("survey_wgs84.geojson", driver="GeoJSON")
    print("bounds:", web.total_bounds)

Alternative Variants

Raw pyproj.Transformer for loose coordinate pairs

When data is not in a GeoDataFrame — for instance a list of easting/northing tuples from a CSV — build the transformer directly. always_xy=True guarantees (x, y) input and output ordering, sidestepping the latitude-first axis surprise. Reuse the transformer across the whole batch rather than reconstructing it per point.

import pyproj

# UTM zone 33N (EPSG:32633) -> WGS84 lon/lat (EPSG:4326).
transformer = pyproj.Transformer.from_crs(
    "EPSG:32633", "EPSG:4326", always_xy=True
)

utm_points = [(500000.0, 5200000.0), (512345.0, 5211987.0)]
lonlat = [transformer.transform(e, n) for e, n in utm_points]
# -> [(lon, lat), ...] ready for GeoJSON [lon, lat] order

UTM zone and EPSG reference

Region example Longitude band Zone North EPSG South EPSG
UK / Portugal −6° to 0° 30 EPSG:32630 EPSG:32730
Central Europe 12° to 18° 33 EPSG:32633 EPSG:32733
US East Coast −78° to −72° 18 EPSG:32618 EPSG:32718
Japan (Tokyo) 138° to 144° 54 EPSG:32654 EPSG:32754
Web target (all) — — EPSG:3857 tiles EPSG:4326 GeoJSON

Verification Steps

  • Source CRS is defined — assert gdf.crs is not None before reprojecting; an undefined CRS silently produces wrong output rather than an error.
  • Zone sanity — cross-check utm_epsg_for_lonlat() against a known point in the dataset; a mismatch means the file was labelled with the wrong zone.
  • Bounds in range — after a 3857 reproject, assert total_bounds sits within ±20037508; after a 4326 reproject, assert it sits within (-180, -90, 180, 90).
  • Visual alignment — export a sample to GeoJSON, drop it on a Leaflet or maplibregl map, and confirm features land on the expected streets, not offset by kilometres.
  • Round-trip check — reproject back to the source UTM zone and confirm coordinates match the originals to sub-meter precision, proving the transform chain is lossless.

Common Errors & Fixes

Features land in the middle of the ocean or on another continent

Raw UTM easting/northing meters were passed to the web renderer without reprojection, or the wrong UTM zone was assigned. A UTM easting near 500000 interpreted as a longitude or as Web Mercator meters lands far off-target. Reproject with .to_crs(epsg=3857) (or 4326 for GeoJSON) and verify the source zone with utm_epsg_for_lonlat() before trusting the output.

ValueError: Cannot transform naive geometries from geopandas

The GeoDataFrame has no CRS attached, so geopandas refuses to guess. Set it explicitly with gdf.set_crs('EPSG:32633', allow_override=True) using the zone the data was actually captured in. Never call to_crs on a naive frame hoping it defaults sensibly — it does not, and the result is silently misaligned.

Coordinates come out swapped (latitude where longitude should be)

A latitude-first CRS axis definition was honoured by the raw pyproj.Transformer. Reconstruct it with always_xy=True to force (x, y) / (longitude, latitude) ordering. When using geopandas .to_crs() this is handled for you, so prefer the GeoDataFrame path whenever the data is already a frame.

Reprojected bounds are correct but tiles still misalign at high zoom

The geometry was reprojected to EPSG:4326 and then handed to a component expecting EPSG:3857 meters, or vice versa. Match the target CRS to the consumer: GeoJSON overlays and folium want EPSG:4326 degrees, while pre-tiled pipelines feeding a tool like tippecanoe — see generating vector tiles from PostGIS with Tippecanoe — also consume EPSG:4326 input and tile into EPSG:3857 internally.