Spatial Epi Public Health GIS

How to Align WGS84 and UTM for County Health Data

This guide solves one specific failure mode: case-level coordinates arriving in WGS84 (EPSG:4326) that must be deterministically joined to county boundaries projected in UTM, without introducing silent metric distortion. It is part of the Coordinate Reference Systems for Public Health cluster within the broader Spatial Epidemiology Fundamentals & Data Standards framework.

Problem Context & Constraints

The naive approach — loading WGS84 points and UTM polygons into the same GeoDataFrame and running a spatial join — fails in ways that are invisible until an audit catches them. GPS-derived telemetry, EHR geocodes, and legacy survey exports almost always land in WGS84, while jurisdictional layers and analytical grids are distributed in UTM (e.g., EPSG:32617 for the eastern United States). The two systems are not interchangeable:

  • WGS84 is a geographic coordinate system (GCS). It stores positions as angular degrees on an ellipsoid. A degree of longitude is not a fixed ground distance — it shrinks from ~111 km at the equator toward zero at the poles.
  • UTM is a conformal projected coordinate system (PCS). It flattens the ellipsoid into metric grids within 6° longitudinal zones, so x/y are meters and distance/area operations are valid.

Overlaying degrees onto a metric grid without an explicit transformation makes distance buffers elliptical, skews area-normalized incidence rates, and causes within predicates to return false negatives — cases silently drop out of their own county. Because many desktop GIS tools reproject on the fly for display only, the rendered map can look correct while the underlying join is wrong. For production epidemiology you must perform the transformation explicitly, log it, and validate the result. This is the same enforcement gate described in the parent CRS guide; here we apply it to the exact WGS84-to-UTM county case.

The specific constraints this scenario imposes:

  1. The transformation must be deterministic — identical inputs always produce identical UTM coordinates and identical join assignments.
  2. The correct UTM zone must be derived from the data, not hardcoded, or cases near a zone seam (e.g., 78°W on the Zone 17/18 boundary) project into the wrong grid.
  3. Coordinate precision must be masked in meters for de-identification before outputs leave the pipeline, matching the thresholds in Setting Decimal Precision for Disease Coordinate Mapping.
WGS84-to-UTM Alignment Pipeline vs. the Naive Same-Table Join Two WGS84 case points and a UTM county boundary layer enter a CRS validation gate. The correct path derives or locks the UTM zone, reprojects the points to metric UTM, masks coordinate precision in meters, runs a strict within spatial join, and writes an audit log so every case is assigned to its county. A contrasting lower branch shows the naive route: degrees are overlaid directly on the metric grid in one table, so within predicates return false negatives and cases silently drop out of their own county. WGS84 case points EPSG:4326 · degrees UTM county polygons EPSG:326xx · meters CRS gate explicit CRS? explicit transform — every case assigned Derive / lock UTM zone from data extent Reproject to UTM (meters) to_crs(epsg=…) Precision mask round meters de-identify Spatial join within 1:1 to county Audit log CRS · zone · counts · unmatched naive route — degrees + meters in one table no transform Same-table join degrees over metric grid within → false negatives Silent dropouts cases fall outside their own county · map looks fine

Prerequisites

Pin the GIS stack so the transformation is reproducible across environments. PROJ ships the datum grids that drive the WGS84/UTM math, so its version is part of your provenance record.

# requirements.txt (pinned for deterministic geodesy)
geopandas==0.14.4
pyproj==3.6.1     # bundles PROJ 9.x datum grids
shapely==2.0.4
pandas==2.2.2

Input data requirements:

  • Case layer: point geometries with an explicit WGS84 CRS, or a documented assumption that they are WGS84. Coordinates that are silently in some other GCS are the single most common source of misaligned joins.
  • County boundary layer: clean, valid polygons with an explicit CRS. A stable, unique county identifier (e.g., GEOID/FIPS) is required so join cardinality and ordering are reproducible.

Step-by-Step Solution

The pipeline below validates CRS declarations, repairs geometry, derives the UTM zone from the data, reprojects deterministically, applies metric precision masking, runs the spatial join, and writes a structured audit trail. It uses only the pinned stack above.

# geopandas==0.14.4, pyproj==3.6.1, shapely==2.0.4, pandas==2.2.2
import logging
import os

import geopandas as gpd
import pandas as pd
from pyproj import CRS
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info
from shapely.validation import make_valid

# Structured audit logging — every run records the CRS transform for provenance
logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(message)s",
    handlers=[logging.FileHandler("crs_alignment_audit.log"), logging.StreamHandler()],
)


def derive_utm_epsg(cases_wgs84: gpd.GeoDataFrame) -> int:
    """Derive the UTM zone from the data extent in WGS84, instead of hardcoding it.

    Locking the zone to the data centroid prevents cases near a 6-degree zone seam
    from being projected into an adjacent grid.
    """
    minx, miny, maxx, maxy = cases_wgs84.total_bounds
    candidates = query_utm_crs_info(
        datum_name="WGS 84",
        area_of_interest=AreaOfInterest(minx, miny, maxx, maxy),
    )
    if not candidates:
        raise ValueError("No UTM zone resolved for the case extent; inspect coordinates.")
    epsg = int(candidates[0].code)
    logging.info("Derived UTM zone: EPSG:%s for bounds %s", epsg, (minx, miny, maxx, maxy))
    return epsg


def align_and_join_county_health_data(
    cases_path: str,
    boundaries_path: str,
    target_utm_epsg: int | None = None,
    precision_decimals: int = 2,
    id_column: str = "GEOID",
) -> gpd.GeoDataFrame:
    """Align WGS84 case coordinates to UTM county boundaries and assign each case to a county.

    target_utm_epsg=None derives the zone from the data; pass an int to lock an agency standard.
    precision_decimals applies to METERS in UTM, so 2 == 0.01 m resolution.
    """
    if not os.path.exists(cases_path) or not os.path.exists(boundaries_path):
        raise FileNotFoundError("Input paths must resolve to valid vector files.")

    cases = gpd.read_file(cases_path)
    counties = gpd.read_file(boundaries_path)

    # 1. CRS validation & explicit WGS84 assignment
    if cases.crs is None or counties.crs is None:
        raise ValueError("Input datasets must contain explicit CRS definitions.")
    if not cases.crs.equals(CRS.from_epsg(4326)):
        logging.warning("Case CRS is %s, not EPSG:4326; reprojecting to WGS84 first.", cases.crs.to_string())
        cases = cases.to_crs(epsg=4326)

    # 2. Geometry repair — prevents sjoin from raising GEOSException on bad polygons
    cases["geometry"] = cases.geometry.apply(lambda g: make_valid(g) if g is not None else g)
    counties["geometry"] = counties.geometry.apply(lambda g: make_valid(g) if g is not None else g)

    # 3. Deterministic zone selection + projection to UTM (meters)
    if target_utm_epsg is None:
        target_utm_epsg = derive_utm_epsg(cases)
    cases_utm = cases.to_crs(epsg=target_utm_epsg)
    counties_utm = counties.to_crs(epsg=target_utm_epsg)

    # 4. De-identification: round Point coordinates in METERS before join/output
    point_mask = cases_utm.geom_type == "Point"
    if point_mask.any():
        cases_utm.loc[point_mask, "geometry"] = cases_utm.loc[point_mask, "geometry"].apply(
            lambda geom: type(geom)(round(geom.x, precision_decimals), round(geom.y, precision_decimals))
        )

    # 5. Spatial join — strict 'within' for jurisdictional assignment
    joined = gpd.sjoin(cases_utm, counties_utm, how="left", predicate="within")

    # 6. Validation & audit reporting
    unjoined = int(joined["index_right"].isna().sum())
    if unjoined:
        logging.warning("Unmatched cases: %d. Review zone seams, datum drift, or boundary topology.", unjoined)
    logging.info(
        "Alignment complete | target=EPSG:%s | records=%d | unmatched=%d",
        target_utm_epsg, len(joined), unjoined,
    )

    # Deterministic ordering for reproducible downstream aggregation
    if id_column in joined.columns:
        joined = joined.sort_values(id_column, kind="mergesort").reset_index(drop=True)
    return joined

Validation & Edge Cases

Three failure modes account for nearly all production incidents on this exact join. Each has a concrete diagnostic.

1. Zone-boundary splits. A case at 78°W sits on the Zone 17/18 seam. If the EPSG is hardcoded to the wrong zone, the point projects 100+ km off-grid and falls outside every county polygon. Derive the zone from the data (as above) or, when locking an agency standard, flag coordinates outside the nominal zone:

# Quarantine cases whose longitude falls outside the locked UTM zone's +/-3 deg span
zone_central_lon = -81.0  # Zone 17 central meridian
out_of_zone = cases[(cases.geometry.x < zone_central_lon - 3) |
                    (cases.geometry.x > zone_central_lon + 3)]
# log output:
# WARNING | 412 cases outside locked zone span; routed to quarantine table 'cases_zone_review'

2. Invalid geometries crashing the join. Self-intersecting boundary polygons or degenerate points make sjoin raise a GEOSException. The make_valid() calls prevent the crash, but you should still assert validity and quarantine survivors:

assert cases_utm.geometry.is_valid.all(), "Invalid case geometries remain after repair"
bad = counties_utm[~counties_utm.geometry.is_valid]
# log output:
# ERROR | shapely.errors.GEOSException: TopologyException: side location conflict
#   -> resolved by make_valid(); 0 invalid county polygons remain

3. Predicate mismatch and join cardinality. County assignment must be 1:1. A many-to-one result means overlapping administrative layers or boundary topology errors, not a coordinate problem. Verify cardinality and choose the predicate deliberately — use within for strict assignment, and switch to intersects only for intentionally jittered de-identified points (and document the tolerance):

dupes = joined.index[joined.index.duplicated()]
assert dupes.empty, f"{len(dupes)} cases matched multiple counties — inspect overlapping boundaries"
# log output:
# WARNING | 7 cases matched 2 counties each; boundary layer has overlapping slivers at FIPS 37119/37025

Compliance Notes

This join is a coordinate transformation on protected health data, so it must be defensible in an audit. Record, per run: the input CRS of each layer, the resolved target_utm_epsg (and whether it was derived or locked), the transformation method, the precision_decimals applied, and input/output record counts plus the unmatched count. The pipeline writes all of these to crs_alignment_audit.log.

Metric precision masking (rounding to 0.01 m at precision_decimals=2, or coarser for public dashboards) reduces re-identification risk while preserving county-level analytical utility; align the chosen threshold with Setting Decimal Precision for Disease Coordinate Mapping. Never persist raw WGS84 coordinates beside transformed UTM outputs in the same analytical table without an explicit de-identification flag. Maintain version-controlled boundary files with SHA-256 checksums so unauthorized edits to jurisdictional layers are detectable, and sort outputs by a stable identifier (above) so aggregation is reproducible across runs.

A note on datums: WGS84 and NAD83 differ by under ~2 m across most of the continental US, but cross-border or legacy survey data may need an explicit datum shift. Verify whether additional parameters are required against the PROJ projection documentation before treating the two as equivalent.

Frequently Asked Questions

Should I project to UTM or just compute everything in WGS84? Any operation that produces distances, buffers, or areas — incidence per km², exposure rings, distance-based spatial weights — must run in a metric CRS. Projecting points to UTM before the join keeps those downstream operations valid. Use WGS84 only as an interchange/storage format.

Why derive the UTM zone instead of hardcoding EPSG:32617? Hardcoding is safe only when every case is guaranteed to fall inside one zone. County datasets that straddle a 6° seam (or agencies that ingest multi-state data) will misproject edge cases into the wrong grid. Deriving the zone from the data extent, then logging it, gives you both correctness and a provenance record.

Cases are dropping out of the join — what should I check first? In order: (1) the case CRS is actually WGS84 and not a different GCS; (2) the resolved UTM zone matches the data, not a stale hardcoded value; (3) county polygons are valid after make_valid(); (4) the predicate is appropriate (within vs intersects for jittered points). The audit log’s unmatched count plus the zone-quarantine output isolate which of these is responsible.

What precision should I round UTM coordinates to? Two decimal places (0.01 m) preserves full analytical fidelity for county-level work; coarser rounding is appropriate for public-facing outputs. Match the threshold to the analysis scale and the governing data use agreement rather than defaulting to maximum precision.