Picking a UTM Zone for an OSM Extract Jump to heading#

#!/usr/bin/env python3
"""Derive the right projected CRS for an OSM extract from its own bounds."""
from __future__ import annotations

import json
import logging
import math
import subprocess
from dataclasses import dataclass

from pyproj import CRS, Transformer

logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
logger = logging.getLogger(__name__)

#: Beyond this span a single transverse Mercator zone cannot cover the extract
#: without unacceptable edge distortion.
MAX_UTM_SPAN_DEGREES = 4.0


@dataclass(frozen=True)
class Bounds:
    west: float
    south: float
    east: float
    north: float

    @property
    def centre_lon(self) -> float:
        return (self.west + self.east) / 2

    @property
    def centre_lat(self) -> float:
        return (self.south + self.north) / 2

    @property
    def span_lon(self) -> float:
        return self.east - self.west


def bounds_from_pbf(path: str) -> Bounds:
    """The header carries the bbox — one block read, no full pass."""
    info = json.loads(subprocess.run(
        ["osmium", "fileinfo", "--extended", "--json", path],
        capture_output=True, text=True, check=True).stdout)
    box = info["data"]["bbox"]
    return Bounds(west=box[0], south=box[1], east=box[2], north=box[3])


def utm_zone(lon: float) -> int:
    """Zone 1 starts at 180°W; each zone spans 6°."""
    return int(math.floor((lon + 180.0) / 6.0) % 60) + 1


def utm_epsg(lon: float, lat: float) -> int:
    """326xx for the northern hemisphere, 327xx for the southern."""
    return (32600 if lat >= 0 else 32700) + utm_zone(lon)


def choose_crs(bounds: Bounds) -> CRS:
    """UTM when the extract fits a zone; an equal-area projection when it does not."""
    if bounds.span_lon <= MAX_UTM_SPAN_DEGREES:
        epsg = utm_epsg(bounds.centre_lon, bounds.centre_lat)
        crs = CRS.from_epsg(epsg)
        logger.info("extract spans %.2f° — using %s (%s)",
                    bounds.span_lon, crs.name, f"EPSG:{epsg}")
        return crs

    # Too wide for one zone: a Lambert azimuthal equal-area projection centred on
    # the extract distorts shape a little everywhere instead of a lot at the edges.
    crs = CRS.from_proj4(
        f"+proj=laea +lat_0={bounds.centre_lat:.4f} +lon_0={bounds.centre_lon:.4f} "
        f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
    logger.warning("extract spans %.2f° — too wide for UTM; using LAEA centred on "
                   "%.3f, %.3f", bounds.span_lon, bounds.centre_lat, bounds.centre_lon)
    return crs


def make_transformer(bounds: Bounds) -> Transformer:
    """always_xy keeps the argument order (lon, lat) on both sides."""
    return Transformer.from_crs(CRS.from_epsg(4326), choose_crs(bounds), always_xy=True)


def scale_error_ppm(crs: CRS, lon: float, lat: float) -> float:
    """Measured, not assumed: project a 1 km east-west line and compare its length."""
    fwd = Transformer.from_crs(CRS.from_epsg(4326), crs, always_xy=True)
    geod = CRS.from_epsg(4326).get_geod()
    delta = 0.01                                         # ~700 m at 52°N
    x1, y1 = fwd.transform(lon, lat)
    x2, y2 = fwd.transform(lon + delta, lat)
    projected = math.hypot(x2 - x1, y2 - y1)
    _, _, true_distance = geod.inv(lon, lat, lon + delta, lat)
    return (projected / true_distance - 1.0) * 1e6


if __name__ == "__main__":
    b = bounds_from_pbf("ireland.osm.pbf")
    crs = choose_crs(b)
    for lon in (b.west, b.centre_lon, b.east):
        logger.info("scale error at lon %.2f: %+.0f ppm", lon,
                    scale_error_ppm(crs, lon, b.centre_lat))

scale error at lon -10.48: +112 ppm
scale error at lon  -8.00: -398 ppm
scale error at lon  -5.30: +121 ppm

fwd = make_transformer(b)
inv = Transformer.from_crs(choose_crs(b), CRS.from_epsg(4326), always_xy=True)
for lon, lat in ((b.west, b.south), (b.centre_lon, b.centre_lat), (b.east, b.north)):
    x, y = fwd.transform(lon, lat)
    back_lon, back_lat = inv.transform(x, y)
    assert abs(back_lon - lon) < 1e-9 and abs(back_lat - lat) < 1e-9