Spatial Epi Public Health GIS

Reconciling ESRI and OSRM CRS in Isochrone Generation

When you overlay an ESRI Network Analyst service area on an open-source OSRM isochrone during drive-time isochrone generation, the two polygons that should describe the same catchment can drift apart by tens to hundreds of metres, silently inflating or deflating the population a facility appears to cover. This guide resolves that misalignment in Python by normalizing both sources to one projected equal-area coordinate system, handling axis order explicitly, and area-checking the result.

Problem Context & Constraints

The failure is subtle because both files often look correct in isolation. An ESRI Network Analyst polygon frequently arrives tagged as WKID 102100 (ESRI’s label for Web Mercator) or in a State Plane zone measured in US survey feet, while a self-hosted OSRM isochrone comes back as WGS84 longitude/latitude, exactly the EPSG:4326 the engine consumes. Load both without thinking and a GIS library reprojects each according to its declared authority — but three classic traps mean the declared authority is not always the effective one.

The first trap is the 102100↔3857 identity. ESRI’s WKID 102100 and the modern EPSG:3857 describe the same spherical Web Mercator, but older toolchains registered 102100 with slightly different WKT and, historically, a different treatment of the datum. If your pyproj install resolves 102100 through an ESRI WKT that pins a sphere while another layer resolves 3857 through the EPSG definition, a transformation pipeline can be inserted where none is needed, nudging vertices by a few metres. The safe move is to recognize 102100 as an alias and force both to a single authoritative code before any transform runs.

The second trap is axis order. EPSG:4326 is formally defined as latitude, longitude — the opposite of the longitude, latitude ordering that OSRM emits and that most desktop GIS displays. If an isochrone’s coordinates were serialized lon/lat but the file (or a downstream reader) asserts the CRS is EPSG:4326 with authority-compliant axis order, the polygon lands with its X and Y swapped, dropping the catchment in the wrong hemisphere or, near 40°N, shearing it by hundreds of metres. A wall of implausible coordinates is the visible symptom; a modest shift that still plots “near” the facility is the dangerous one.

The third trap is unit and datum mismatch on State Plane sources. A Network Analyst output in a State Plane zone reported in US survey feet, reprojected by a reader that assumes international feet or metres, scales the whole polygon around the projection origin. Close to the facility the error is small; at the catchment edge — exactly where the coverage count is decided — it is largest.

The naive reconciliation, gdf_b.to_crs(gdf_a.crs), trusts every one of those declarations. Because coverage counts feed resource allocation, an undetected shift is not a cosmetic bug: it moves the boundary between “served” and “underserved” populations. The constraint for public health work is therefore that reconciliation must be explicit, logged, and verified by an independent area check — never assumed from metadata alone.

The reconciliation flow below normalizes both sources to a common equal-area target, repairs axis order, and gates on an area-delta check before the polygons are allowed to merge:

Reconciling ESRI and OSRM Isochrones to a Common CRS Three panels left to right. Panel one shows two mis-registered catchment blobs: an ESRI Network Analyst polygon tagged WKID 102100 and an OSRM polygon in WGS84 lon/lat, offset from each other. A reconcile stage in the middle normalizes authority codes, repairs lat/lon axis order, and reprojects both to EPSG:5070 equal-area. Panel three shows the two polygons aligned, passing an area-delta check before they merge. Mis-registered same catchment, two sources ESRI 102100 OSRM WGS84 offset by tens–hundreds of m Reconcile map 102100 → 3857 repair lat/lon axis order to_crs(EPSG:5070) equal-area, metres Aligned area-delta check passes one boundary |ΔA| / A < tolerance Explicit authority mapping and axis repair, then an area check, turn two drifting polygons into one catchment

Prerequisites

Pin the stack so the transformation pipeline is reproducible across analysts and reviewers:

  • python 3.11
  • geopandas 0.14.4 (pulls shapely 2.0.x)
  • pyproj 3.6.1 — the CRS authority resolver and transformer
  • numpy 1.26.4

Input state assumed by the code below:

  • ESRI isochrone: a polygon layer exported from Network Analyst, commonly tagged WKID 102100 or a State Plane zone (e.g. EPSG:2229, California Zone V in US survey feet). Confirm the declared authority in the layer metadata, not by eye.
  • OSRM isochrone: a polygon in EPSG:4326, coordinates serialized in (longitude, latitude) order as OSRM emits them.
  • A single target CRS chosen for metric, equal-area comparison. For a multi-state or national footprint, use EPSG:5070 (CONUS Albers Equal Area); for a single metropolitan study area a local UTM zone is acceptable. The canonical selection rules live in coordinate reference systems for public health, and getting the WGS84-to-projected step exactly right is covered in how to align WGS84 and UTM for county health data.

Step-by-Step Solution

The routine below resolves each source’s authority code explicitly (mapping the 102100 alias), detects and repairs a swapped-axis WGS84 layer, reprojects both to the equal-area target, and returns the reconciled pair with the authority codes logged for audit.

# python 3.11
# geopandas==0.14.4  shapely==2.0.4  pyproj==3.6.1  numpy==1.26.4
import logging
import numpy as np
import geopandas as gpd
from pyproj import CRS

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("crs_reconcile")

# ESRI WKID 102100 is a documented alias of EPSG:3857 (spherical Web Mercator).
ESRI_ALIAS = {102100: 3857, 102113: 3857}

def canonical_crs(gdf: gpd.GeoDataFrame) -> CRS:
    """Resolve a layer's declared CRS to an authoritative EPSG code, mapping ESRI aliases."""
    if gdf.crs is None:
        raise ValueError("Layer has no CRS; refusing to guess. Set the source authority explicitly.")
    auth = gdf.crs.to_authority()          # e.g. ('ESRI', '102100') or ('EPSG', '4326')
    if auth and auth[0] == "ESRI" and int(auth[1]) in ESRI_ALIAS:
        code = ESRI_ALIAS[int(auth[1])]
        log.info("mapped ESRI:%s -> EPSG:%d (Web Mercator identity)", auth[1], code)
        return CRS.from_epsg(code)
    return gdf.crs

def looks_axis_swapped(gdf: gpd.GeoDataFrame) -> bool:
    """Heuristic: a lon/lat layer mislabeled as lat/lon lands with |X| <= 90 and |Y| > 90."""
    b = gdf.total_bounds  # (minx, miny, maxx, maxy)
    x_in_lat_range = np.abs([b[0], b[2]]).max() <= 90.0
    y_out_of_lat_range = np.abs([b[1], b[3]]).max() > 90.0
    return bool(x_in_lat_range and y_out_of_lat_range)

def reconcile_isochrones(esri: gpd.GeoDataFrame, osrm: gpd.GeoDataFrame,
                         target_epsg: int = 5070) -> tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]:
    """Normalize two isochrone sources to one projected equal-area CRS with explicit handling."""
    esri = esri.set_crs(canonical_crs(esri), allow_override=True)
    osrm = osrm.set_crs(canonical_crs(osrm), allow_override=True)

    # Axis-order repair: OSRM serializes (lon, lat); a reader asserting authority-compliant
    # EPSG:4326 axis order can swap them. Detect and flip before any reprojection.
    if osrm.crs.to_epsg() == 4326 and looks_axis_swapped(osrm):
        log.warning("OSRM layer appears axis-swapped (lat/lon); flipping to (lon, lat)")
        osrm = osrm.set_geometry(osrm.geometry.apply(_swap_xy))

    target = CRS.from_epsg(target_epsg)
    esri_p = esri.to_crs(target)   # equal-area, metres
    osrm_p = osrm.to_crs(target)
    log.info("reconciled: esri=%s osrm=%s -> EPSG:%d",
             esri.crs.to_authority(), osrm.crs.to_authority(), target_epsg)
    return esri_p, osrm_p

def _swap_xy(geom):
    """Swap X/Y ordinates of a (Multi)Polygon using shapely's transform."""
    from shapely.ops import transform
    return transform(lambda x, y, z=None: (y, x), geom)

Two design choices make the routine auditable. First, canonical_crs never guesses a missing CRS — it raises — so a layer with no declared authority cannot slip through with a silent default. Second, the axis-swap repair is a logged, conditional branch keyed on an objective bounds test, not a blanket flip that would corrupt an already-correct layer.

Validation & Edge Cases

A reprojected pair is not yet a reconciled pair. Validate the geometry with an independent area check before merging, and watch for these three failure modes.

1. Area delta between the two sources. After both isochrones sit in the same equal-area CRS, the same catchment should enclose nearly the same area. A large relative delta is the signature of an unresolved CRS problem — a missed unit conversion or a residual axis swap — not a genuine modeling difference:

def area_delta(esri_p, osrm_p):
    """Relative area difference in the shared equal-area CRS. Flag > 5%."""
    a_esri = float(esri_p.geometry.area.sum())
    a_osrm = float(osrm_p.geometry.area.sum())
    rel = abs(a_esri - a_osrm) / max(a_esri, a_osrm)
    log.info("area_check esri=%.1f m2 osrm=%.1f m2 rel_delta=%.4f", a_esri, a_osrm, rel)
    if rel > 0.05:
        log.warning("area delta %.1f%% exceeds tolerance; suspect CRS/unit/axis error", rel * 100)
    return rel
WARNING area delta 41.3% exceeds tolerance; suspect CRS/unit/axis error

A ~40% delta is the classic US-survey-foot-read-as-metre scaling error; a delta that resolves only after _swap_xy confirms an axis problem.

2. The 102100 pipeline that shifts by metres. If canonical_crs did not fire — because the layer was tagged EPSG-side but carried an ESRI WKT internally — a transformation pipeline can insert a null datum shift that displaces vertices by 2–5 m. Diagnose by comparing pyproj’s chosen transform explicitly:

from pyproj import Transformer
t = Transformer.from_crs(esri.crs, CRS.from_epsg(5070), always_xy=True)
log.info("pipeline=%s", t.description)   # inspect for an unexpected datum step

Force both Web Mercator variants to EPSG:3857 before the projected transform so the pipeline is a clean planar reprojection, not a datum hop.

3. Snap tolerance after alignment. Even correctly reconciled polygons carry sub-metre vertex drift from independent routing engines. A light snap to a common grid removes spurious slivers when you later intersect the two boundaries, so an overlay does not fracture a shared edge into thousands of micro-polygons:

from shapely import set_precision
osrm_p["geometry"] = set_precision(osrm_p.geometry.values, grid_size=0.5)  # 0.5 m grid
esri_p["geometry"] = set_precision(esri_p.geometry.values, grid_size=0.5)

Compliance Notes

Three artifacts make a reconciled isochrone defensible in a public health review:

  • Explicit authority provenance. Log the resolved authority code of every input layer and the target EPSG, including any ESRI-to-EPSG alias mapping applied. The record proves which CRS each catchment actually used, so a reviewer can distinguish a genuine coverage difference from a projection artifact.
  • The area-delta metric. Persist the relative area check with each reconciled pair. Because coverage counts drive allocation, the area check is the auditable evidence that the two sources describe the same catchment rather than two silently misregistered ones.
  • Deterministic ordering and pinned versions. pyproj transformation pipelines can change between grid-shift releases; record the pyproj and PROJ versions alongside the output so a re-run reproduces the same vertices. Combined with a stable feature sort, this makes the reconciled catchment byte-reproducible for interagency handoff.