Python Script: Convert WGS84 to Local CRS for Fleet Routing

This page addresses a focused implementation question within Coordinate Reference System Mapping for Fleet Data: how to write the actual Python code that converts a batch of raw GPS telemetry records from WGS84 (EPSG:4326) into a local projected coordinate system, ready for routing solvers and spatial indexing. The conversion is a mandatory first step inside any GPS data preprocessing pipeline — routing engines, Euclidean distance matrices, and KD-tree nearest-neighbour lookups all assume planar, meter-based coordinates and produce meaningless results when fed decimal degrees. Getting always_xy, datum handling, and output validation right the first time prevents a class of bugs that are notoriously hard to detect at runtime because the coordinates look plausible but are geometrically wrong.

Coordinate-flow overview

The diagram below shows where the WGS84-to-local-CRS step sits relative to the broader ingestion and matching pipeline, and which artefacts it produces.

Compatibility and configuration requirements

Install into a clean virtual environment:

pip install "pyproj>=3.0.0" "pandas>=1.3.0"

Production-ready conversion function

The function below is self-contained and handles empty DataFrames, NaN coordinates, and wrong-hemisphere zone selection. Copy it verbatim into your ingestion layer.

import math
import pandas as pd
import numpy as np
from pyproj import Transformer


def utm_epsg_for_centroid(lat: float, lon: float) -> int:
    """
    Return the EPSG code for the UTM zone that best covers a geographic centroid.

    Parameters
    ----------
    lat : float
        Centroid latitude in decimal degrees (WGS84).
    lon : float
        Centroid longitude in decimal degrees (WGS84).

    Returns
    -------
    int
        EPSG code: 326xx for northern hemisphere, 327xx for southern.
    """
    zone_number = math.floor((lon + 180.0) / 6.0) + 1
    # Norway/Svalbard special cases (UTM standard exceptions)
    if 56.0 <= lat < 64.0 and 3.0 <= lon < 12.0:
        zone_number = 32
    if 72.0 <= lat <= 84.0:
        if 0.0 <= lon < 9.0:
            zone_number = 31
        elif 9.0 <= lon < 21.0:
            zone_number = 33
        elif 21.0 <= lon < 33.0:
            zone_number = 35
        elif 33.0 <= lon < 42.0:
            zone_number = 37
    base = 32600 if lat >= 0 else 32700
    return base + zone_number


def convert_wgs84_to_local_crs(
    df: pd.DataFrame,
    lat_col: str = "latitude",
    lon_col: str = "longitude",
    target_epsg: int | None = None,
    auto_utm: bool = True,
) -> tuple[pd.DataFrame, int]:
    """
    Convert WGS84 (EPSG:4326) coordinates to a local projected CRS.

    The function adds two columns — ``x_meters`` and ``y_meters`` — to a copy
    of the input DataFrame and returns both the enriched DataFrame and the EPSG
    code that was used, so callers can cache and reuse the transformer.

    Parameters
    ----------
    df : pd.DataFrame
        Input telemetry frame. Must contain ``lat_col`` and ``lon_col``.
    lat_col : str
        Name of the latitude column.  Default: ``"latitude"``.
    lon_col : str
        Name of the longitude column. Default: ``"longitude"``.
    target_epsg : int | None
        Explicit EPSG code for the target projection (e.g. 27700 for British
        National Grid, 32632 for UTM Zone 32N).  When ``None`` and
        ``auto_utm=True`` the function computes the optimal UTM zone from the
        centroid of the non-null coordinate pairs.
    auto_utm : bool
        When ``True`` and ``target_epsg`` is not supplied, auto-select the UTM
        zone from the data centroid.  Set to ``False`` to force an explicit
        ``target_epsg``.

    Returns
    -------
    tuple[pd.DataFrame, int]
        A copy of ``df`` with ``x_meters`` and ``y_meters`` columns appended,
        plus the EPSG code used for the transform (useful for caching).

    Raises
    ------
    ValueError
        If ``target_epsg`` is None and ``auto_utm`` is False, or if the
        coordinate columns are missing from the DataFrame.
    """
    for col in (lat_col, lon_col):
        if col not in df.columns:
            raise ValueError(f"Column '{col}' not found in DataFrame.")

    if df.empty:
        out = df.copy()
        out["x_meters"] = np.nan
        out["y_meters"] = np.nan
        epsg = target_epsg or 32632
        return out, epsg

    if target_epsg is None:
        if not auto_utm:
            raise ValueError(
                "Supply a target_epsg or set auto_utm=True to derive the UTM zone automatically."
            )
        valid = df[[lat_col, lon_col]].dropna()
        if valid.empty:
            raise ValueError("All coordinate rows are NaN — cannot determine UTM zone.")
        centroid_lat = valid[lat_col].mean()
        centroid_lon = valid[lon_col].mean()
        target_epsg = utm_epsg_for_centroid(centroid_lat, centroid_lon)

    # always_xy=True: enforce (longitude, latitude) input order regardless of
    # the CRS authority's official axis definition.  Without this flag, some
    # projected CRS definitions swap the axes and the output grid is silently
    # transposed — vehicles appear in the wrong country.
    transformer = Transformer.from_crs(
        "EPSG:4326",
        f"EPSG:{target_epsg}",
        always_xy=True,
    )

    # Vectorized transform: pyproj delegates to C-level PROJ routines.
    # Passing .values bypasses pandas index overhead; NaN coordinates pass
    # through unchanged (pyproj 3.0+ behaviour).
    x, y = transformer.transform(
        df[lon_col].values,
        df[lat_col].values,
    )

    out = df.copy()
    out["x_meters"] = x
    out["y_meters"] = y
    return out, target_epsg


def validate_projection(df: pd.DataFrame, epsg: int) -> None:
    """
    Assert that projected coordinates fall within expected UTM bounding boxes.

    Routing solvers that receive out-of-range coordinates produce unbounded
    distance matrices.  Call this immediately after ``convert_wgs84_to_local_crs``.

    Parameters
    ----------
    df : pd.DataFrame
        DataFrame containing ``x_meters`` and ``y_meters`` columns.
    epsg : int
        The EPSG code used for the projection.
    """
    # UTM easting is always in [100 000, 900 000] m; northing [0, 10 000 000] m.
    # State Plane and national grids have wider ranges so we use a broad check.
    x_ok = df["x_meters"].dropna().between(-10_000_000, 10_000_000).all()
    y_ok = df["y_meters"].dropna().between(-10_000_000, 10_000_000).all()
    if not (x_ok and y_ok):
        raise ValueError(
            f"Projected coordinates for EPSG:{epsg} contain out-of-range values. "
            "Check that lat/lon columns are not accidentally swapped in the source data."
        )


# ---------------------------------------------------------------------------
# Example: end-to-end usage
# ---------------------------------------------------------------------------
if __name__ == "__main__":
    sample = pd.DataFrame({
        "vehicle_id": ["V001", "V001", "V002"],
        "latitude":  [51.5074, 51.5085, 51.4901],
        "longitude": [-0.1278, -0.1260, -0.1430],
    })

    # Auto-detect UTM zone from data centroid (UTM Zone 30N → EPSG:32630)
    projected, used_epsg = convert_wgs84_to_local_crs(sample)
    validate_projection(projected, used_epsg)
    print(f"Used EPSG:{used_epsg}")
    print(projected[["vehicle_id", "x_meters", "y_meters"]])

Execution and tuning guidelines

Running the function

Call convert_wgs84_to_local_crs(df) at the ingestion boundary — the single place in your pipeline where raw OBD-II or mobile GPS records first enter your system. Returning the used epsg alongside the DataFrame lets you cache the Transformer instance for subsequent batches from the same fleet region, eliminating the PROJ CRS-lookup overhead on repeated calls.

Key parameter knobs

• target_epsg — hardcode this for deterministic multi-batch pipelines that must guarantee consistent coordinate origins. For a UK urban fleet use 27700 (British National Grid); for a Central European fleet use 32633 (UTM Zone 33N). Leaving it None with auto_utm=True is safe for single-region deployments but can produce different zone selections if a vehicle crosses a UTM zone boundary.

• auto_utm — when True, the function averages all non-null coordinate pairs to find a centroid and derives the optimal UTM zone. This is appropriate for fleets whose stop-set centroid is stable day-to-day. For dynamically growing networks that span multiple UTM zones, override this by setting the zone to the one that minimises worst-case distortion across your operational bounding box.

• always_xy=True inside Transformer.from_crs — this parameter is not a tunable knob; it is always required. Removing it is the single most common source of silent projection errors. The parameter tells pyproj to ignore the CRS authority’s axis-order declaration and always expect (longitude, latitude) input, which is the convention used by GPS hardware, GeoJSON, and most telematics APIs.

Chunked processing for large exports

For ingestion jobs that process raw telematics exports exceeding 10 million rows, instantiate the transformer once outside the loop and process the file in chunks:

from pyproj import Transformer

EPSG_TARGET = 32633
transformer = Transformer.from_crs("EPSG:4326", f"EPSG:{EPSG_TARGET}", always_xy=True)

chunks = []
for chunk in pd.read_csv("telematics_export.csv", chunksize=500_000):
    x, y = transformer.transform(chunk["longitude"].values, chunk["latitude"].values)
    chunk["x_meters"] = x
    chunk["y_meters"] = y
    chunks.append(chunk)

projected_df = pd.concat(chunks, ignore_index=True)

This keeps peak RAM usage flat and avoids repeated CRS registry lookups, which account for most of the latency in naive per-row implementations.

Common pitfalls

• Silent axis transposition without always_xy=True. If you omit the always_xy flag and your target EPSG uses (Northing, Easting) as its official axis order (common in national grids such as EPSG:27700 or EPSG:4258), pyproj will silently swap longitude and latitude. The resulting easting and northing values are numerically plausible but spatially wrong — vehicles appear tens or hundreds of kilometres from their true positions. The validation function above catches this in most cases because the transposed output falls outside the expected coordinate range for the zone.

• Hardcoding a single UTM zone for a multi-region fleet. A VRP solver fed coordinates from a UTM zone that does not cover its operational area will accumulate scale distortion of up to 1:400 near zone edges. This inflates computed distances for vehicles near the zone boundary relative to those near the central meridian, causing the solver to underestimate travel times for edge-area routes. Use auto_utm to let the centroid logic select the appropriate zone, or project into a national grid that covers the full operational area.

• Re-projecting on every query rather than at ingestion. Calling Transformer.transform inside a per-request path (e.g. inside a routing API handler) adds 10–50 ms of overhead per batch, and — more critically — means different requests may use different EPSG codes if the fleet’s active centroid drifts. Project coordinates once at ingestion, persist x_meters/y_meters in your telemetry store alongside the raw WGS84 columns, and expose only the projected values to downstream stop detection and map matching services.

Up: Coordinate Reference System Mapping for Fleet Data | GPS Data Preprocessing & Cleaning Fundamentals

Related

• Coordinate Reference System Mapping for Fleet Data — broader CRS strategy, multi-zone fleet workflows, and storage schema recommendations
• Kalman Filtering for GPS Noise Reduction — smooth noisy positional traces before projecting, so outlier fixes don’t accumulate projection errors
• Outlier Removal in Raw Telematics Streams — strip implausible coordinates prior to projection to prevent transformer NaN bleed-through
• Stop Detection and Dwell Time Analytics — the first major downstream consumer of projected coordinates; clustering algorithms require meter-based inputs
• Trajectory Analysis and Map Matching Techniques — map matching engines expect projected or WGS84 coordinates in a documented axis order; confirm before routing