VehicleMapping HD AV Spatial Data

Converting WGS84 to UTM for AV Pipelines

Implement a zone-aware, batch WGS84 (EPSG:4326) → Universal Transverse Mercator (UTM) transformation for the ingestion stage of an autonomous-vehicle mapping pipeline, where global GNSS telemetry must become a metric planning grid before any topology, occupancy, or perception code runs.

This conversion is the first hard gate in the projection workflow for coordinate reference systems for AVs: naive implementations trigger zone-boundary discontinuities, introduce scale-factor distortion, and accumulate floating-point drift that corrupts lane-level topology modeling downstream. In production it must execute within strict latency budgets while preserving sub-centimeter fidelity across terabyte-scale map-tile datasets, holding the ≤0.05 m RMSE projection budget the parent stage commits to.

Batch projection with a round-trip validation gate that quarantines drifting tiles:

Batch WGS84 → UTM projection with a round-trip validation gate A WGS84 longitude/latitude batch resolves a UTM zone from its trajectory centroid, is batch-transformed with pyproj over NumPy or memmap arrays, then inverse-transformed for a round trip. A gate tests whether the residual is within 1e-4 metres: passing tiles flow to the topology and occupancy grid, failing tiles are quarantined for forensic analysis. INPUT WGS84 lon/lat batch EPSG:4326 · decimal degrees STEP 1 Resolve UTM zone per trajectory centroid · 326xx / 327xx STEP 2–3 Batch transform pyproj.Transformer · NumPy / memmap STEP 4 · GATE round trip Δ ≤ 1e-4 m? pass Metric coords → topology / occupancy local origin subtracted fail Quarantine tile flag for forensic analysis

Prerequisites #

  • Python ≥ 3.10 (the type hints below use tuple[...] builtins).
  • pyproj ≥ 3.6 with a bundled PROJ ≥ 9.2. Pin both in the container image and ship the PROJ data directory read-only.
  • NumPy ≥ 1.24 for contiguous-array transforms and numpy.memmap out-of-core access.
  • Input format: decimal-degree longitude/latitude arrays in WGS84, axis order lon-then-lat, already de-duplicated and time-ordered.
  • Upstream stage: this step runs immediately after GNSS/IMU trajectory ingestion and before topology graph construction, semantic labelling, and occupancy grid generation. It consumes the smoothed pose stream produced upstream and is the metric foundation for everything documented under HD mapping architecture and spatial data standards.
  • Environment: set PROJ_LIB to a read-only cache and run with allow_network=False so transforms are byte-for-byte reproducible inside CI and across worker pools.

Step-by-step #

Step 1 — Resolve the UTM zone per trajectory centroid #

The longitudinal heuristic zone = int((lon + 180) / 6) + 1 is insufficient for safety-critical routing. It fails near zone edges, in high-latitude regions like Svalbard, and across Norway's non-standard UTM zones (32V/34V/36V replace 31V/33V/35V). Resolve the EPSG code from the trajectory centroid, not per point, so a whole tile lands in one projection, and handle the special cases explicitly:

python
def resolve_utm_epsg(lon: float, lat: float) -> int:
    """Return the correct EPSG code for a WGS84 lon/lat point.

    Handles the Svalbard (zone 33X) and Norway (32V) special cases
    defined in the UTM standard.
    """
    zone = int((lon + 180) / 6) + 1

    # Norway special case: zone 32V extends over zone 31V for latitudes 56–64°N
    if 56.0 <= lat < 64.0 and 3.0 <= lon < 12.0:
        zone = 32

    # Svalbard special cases (latitudes 72–84°N)
    if 72.0 <= lat < 84.0:
        if 0.0 <= lon < 9.0:
            zone = 31
        elif 9.0 <= lon < 21.0:
            zone = 33
        elif 21.0 <= lon < 33.0:
            zone = 35
        elif 33.0 <= lon < 42.0:
            zone = 37

    # EPSG: 326xx for the northern hemisphere, 327xx for the southern
    return 32600 + zone if lat >= 0 else 32700 + zone

Feed resolve_utm_epsg(lons.mean(), lats.mean()) to pick target_epsg. Expected output: an integer EPSG code such as 32632 (UTM 32N) or 32733 (UTM 33S).

Step 2 — Batch-transform on contiguous NumPy arrays #

Row-wise coordinate transformations are computationally prohibitive at scale. Python's GIL and per-row overhead make DataFrame apply() unusable in a perception-critical loop. Call pyproj.Transformer.transform() once on whole arrays, with always_xy=True to enforce deterministic lon/lat axis order:

python
import numpy as np
from pyproj import Transformer

def wgs84_to_utm_batch(
    lons: np.ndarray,
    lats: np.ndarray,
    target_epsg: int,
) -> tuple[np.ndarray, np.ndarray]:
    """Vectorized WGS84 → UTM transformation with explicit axis order.

    Args:
        lons: 1-D array of longitudes in decimal degrees.
        lats: 1-D array of latitudes in decimal degrees.
        target_epsg: EPSG code for the target UTM zone (e.g. 32632).

    Returns:
        Tuple of (easting, northing) arrays in metres.
    """
    transformer = Transformer.from_crs(
        "EPSG:4326",
        f"EPSG:{target_epsg}",
        always_xy=True,   # lon first, lat second — never rely on implicit axis order
        allow_network=False,  # deterministic in containerised environments
    )
    return transformer.transform(lons, lats)

always_xy=True is non-negotiable: without it, axis order is CRS-dependent and a silent inversion corrupts every downstream spatial index. Keeping allow_network=False and pinning PROJ_LIB to a read-only cache eliminates cold-start write-lock races in multi-threaded worker pools — see the PROJ environment-variable reference for the containerised deployment knobs. Expected output: two float64 arrays of eastings (~166 000–834 000 m) and northings (0–10 000 000 m).

Step 3 — Stream tiles larger than RAM with numpy.memmap #

For tiles that exceed available heap, operate on memory-mapped coordinate arrays in chunks instead of loading the full dataset:

python
import numpy as np
from pyproj import Transformer

def transform_memmap_tile(
    lon_mmap: np.memmap,
    lat_mmap: np.memmap,
    target_epsg: int,
    chunk_size: int = 500_000,
) -> tuple[np.ndarray, np.ndarray]:
    """Out-of-core batch transform for large tile datasets."""
    transformer = Transformer.from_crs(
        "EPSG:4326", f"EPSG:{target_epsg}", always_xy=True, allow_network=False
    )
    n = len(lon_mmap)
    eastings = np.empty(n, dtype=np.float64)
    northings = np.empty(n, dtype=np.float64)

    for start in range(0, n, chunk_size):
        sl = slice(start, start + chunk_size)
        eastings[sl], northings[sl] = transformer.transform(lon_mmap[sl], lat_mmap[sl])

    return eastings, northings

Construct the Transformer once outside the loop — re-instantiating it per chunk reloads PROJ grids and dominates wall-clock time. A 500 000-point chunk keeps the transient buffer near 4 MB per coordinate axis, well under a 500 MB per-worker ceiling.

Step 4 — Subtract a local origin for downstream consumers #

UTM northings near 5–6 million metres burn precision in IEEE 754 double arithmetic once perception code multiplies them through transform matrices. Subtract a fixed per-tile origin (the zone centroid) before handing coordinates to planners or to the spatial index, and record the offset in the tile metadata so the lane-level topology model can reconstruct global coordinates on demand. This local-origin convention is also what keeps multi-sensor coordinate alignment numerically stable when LiDAR and camera frames are fused into the same metric grid during sensor fusion and spatial data alignment.

Verification and acceptance criteria #

Repeated forward and inverse projections accumulate floating-point error that can exceed the 0.1 m tolerance required for lane-centerline alignment. After projecting, inverse-transform back to WGS84 and measure the residual. Tiles whose maximum residual exceeds 1e-4 m must be quarantined, not silently merged:

python
def validate_round_trip(
    lons_orig: np.ndarray,
    lats_orig: np.ndarray,
    eastings: np.ndarray,
    northings: np.ndarray,
    target_epsg: int,
    tolerance_m: float = 1e-4,
) -> tuple[bool, float]:
    """Inverse-project back to WGS84 and check residuals.

    Returns:
        (passed, max_error_m) — max_error_m is in metres.
    """
    inv_transformer = Transformer.from_crs(
        f"EPSG:{target_epsg}", "EPSG:4326", always_xy=True, allow_network=False
    )
    lons_rt, lats_rt = inv_transformer.transform(eastings, northings)

    # Approximate metre-scale residual in the geodetic domain
    # 1° lat ≈ 111 320 m; cos(lat) factor for longitude
    dlat_m = (lats_rt - lats_orig) * 111_320.0
    dlon_m = (lons_rt - lons_orig) * 111_320.0 * np.cos(np.radians(lats_orig))
    errors_m = np.hypot(dlat_m, dlon_m)
    max_error = float(errors_m.max())

    return max_error <= tolerance_m, max_error

Acceptance gate, wired into CI on a fixed survey-point fixture:

  • passed is True and max_error_m ≤ 1e-4 for every tile in the regression set.
  • Resolved EPSG code matches the surveyed zone for each fixture trajectory (catches off-by-one zone resolution).
  • Easting falls within 166_000–834_000 m; out-of-band values signal a wrong-zone projection.

Cross-zone trajectory stitching #

When a trajectory straddles a UTM zone boundary, projecting the whole path into either adjacent zone introduces growing scale-factor distortion in the far half. For cross-zone workloads, project overlapping segments into a shared local tangent plane (ENU centred on the boundary crossing) or a custom transverse Mercator CRS spanning both zones. This prevents topological fractures when merging tiles across longitudinal zone boundaries — a recurring failure mode in fleet-scale HD mapping architecture and spatial data standards.

Common errors and fixes #

Silent axis inversion (always_xy omitted). Eastings and northings come out swapped, and the spatial index appears to work until tiles "rotate" 90°. There is no exception. Diagnosis: assert that eastings sit in 166_000–834_000 m. Fix: always construct the transformer with always_xy=True.

pyproj.exceptions.CRSError: Invalid projection on the resolved EPSG. Caused by passing a southern-hemisphere latitude through northern-zone logic, or an out-of-range zone number. Fix: route every code through resolve_utm_epsg() and assert 32601 <= epsg <= 32660 or 32701 <= epsg <= 32760 before building the transformer.

RuntimeError/write-lock races from PROJ grid cache under threads. Manifests as intermittent transform failures only under load in a worker pool. Diagnosis: workers race to write the PROJ cache. Fix: set allow_network=False, pin PROJ_LIB to a pre-populated read-only directory, and warm the cache during image build.

Round-trip residual above 1e-4 m on a "clean" tile. Usually a datum mismatch — coordinates were not actually WGS84, or a geoid/ellipsoid height crept into the planar transform. Fix: confirm the source CRS is EPSG:4326, strip any height channel before the planar transform, and re-run validate_round_trip; if it still fails, quarantine and inspect with the Open Geospatial Consortium coordinate-transformation standard.

Up one level: Coordinate reference systems for AVs.