Spatial Epi Public Health GIS

Picking a CRS for Multi-State Surveillance Regions

This guide, part of Coordinate Reference Systems for Public Health, addresses the projection choice for surveillance that crosses state lines — a regional outbreak response, a multi-state vector program, a health information exchange — where every participant’s native system is the wrong one for the whole.

Problem Context & Constraints

Within one state the choice is usually made for you: the state plane zone or a state-wide equal-area system is standard, the boundary files arrive in it, and distortion is small across the state’s extent. A multi-state region breaks all three. No state plane zone is designed for a neighbouring state; UTM zones split most regions down the middle; and a national equal-area system, while valid, distorts more at the region’s edges than a projection centred on the region would.

The distortion matters for the operations regional surveillance actually performs. Kernel density bandwidths, buffer radii, scan statistic radii and nearest-facility distances are all specified in metres, and they are only metres where the projection preserves distance. Rates by area depend on area, which only an equal-area projection preserves. A projection whose scale error reaches half a per cent at the region’s edges shrinks or stretches a 10 km scan radius by 50 m there — small, but systematic and geographically structured, which is the worst kind of error for a spatial statistic.

Three constraints shape the choice. Equal area is usually non-negotiable, because rates and densities depend on it. Distortion should be measured where the population is, not over the bounding box, since an error in an empty desert corner affects nothing. And the handoff between each state’s native system and the regional one must be explicit and lossless, because participating states will keep producing data in their own systems.

Distortion Across the Region Under Three Choices Scale error plotted along a west-to-east transect crossing three states. One state's plane zone has almost no error inside its own state and grows steeply across the neighbours. A national equal-area system has a modest error that varies gently across the whole region. A custom equal-area projection with standard parallels chosen for the region keeps the error near zero across all three states. A shaded band marks where the region's population is concentrated, which is where the comparison should be judged. Judge the projection where the people are, not over the bounding box population concentrated here one state's plane zone national system region-centred state A state B state C absolute scale error along a west-to-east transect scale error

Prerequisites

  • python 3.11
  • pyproj 3.6.1, geopandas 0.14.4, numpy 1.26.4
  • The region’s boundary and a population-weighted point sample or block centroids with population
  • The native coordinate systems and realizations each participating state uses
  • Agreement among participants on a single datum realization for the regional system

Step-by-Step Solution

Define a region-centred Albers equal-area projection with standard parallels at roughly one sixth and five sixths of the region’s latitude span, measure its scale error weighted by population, compare it with the candidates, and publish the definition for all participants.

# Region-centred equal-area CRS and population-weighted distortion comparison.
# Pinned: pyproj==3.6.1, geopandas==0.14.4, numpy==1.26.4
import logging

import numpy as np
import geopandas as gpd
from pyproj import CRS, Proj

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


def regional_albers(region: gpd.GeoDataFrame, datum: str = "NAD83(2011)") -> CRS:
    """Albers equal-area with standard parallels placed for the region's latitude span."""
    geo = region.to_crs("EPSG:6318")
    minx, miny, maxx, maxy = geo.total_bounds
    span = maxy - miny
    lat1 = miny + span / 6.0
    lat2 = maxy - span / 6.0
    lon0 = (minx + maxx) / 2.0
    lat0 = (miny + maxy) / 2.0
    wkt = (f"+proj=aea +lat_1={lat1:.4f} +lat_2={lat2:.4f} +lat_0={lat0:.4f} "
           f"+lon_0={lon0:.4f} +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs")
    crs = CRS.from_proj4(wkt)
    log.info("regional Albers: lat1=%.3f lat2=%.3f lon0=%.3f (datum intent %s)",
             lat1, lat2, lon0, datum)
    return crs


def population_weighted_distortion(points: gpd.GeoDataFrame, crs: CRS,
                                   weight_col: str = "population") -> dict:
    """Scale error at each point, weighted by where people live."""
    geo = points.to_crs("EPSG:6318")
    p = Proj(crs)
    f = p.get_factors(geo.geometry.x.to_numpy(), geo.geometry.y.to_numpy())
    # For an equal-area projection areal scale is 1; linear scale varies by direction.
    lin_err = np.maximum(np.abs(np.asarray(f.meridional_scale) - 1.0),
                         np.abs(np.asarray(f.parallel_scale) - 1.0))
    area_err = np.abs(np.asarray(f.areal_scale) - 1.0)
    w = points[weight_col].to_numpy(float)
    w = w / w.sum()
    out = {"linear_err_weighted_mean": float(w @ lin_err),
           "linear_err_p99_by_pop": float(np.quantile(lin_err, 0.99)),
           "areal_err_max": float(area_err.max())}
    log.info("%s: pop-weighted linear error %.5f, p99 %.5f, max areal error %.2e",
             crs.name or "custom", out["linear_err_weighted_mean"],
             out["linear_err_p99_by_pop"], out["areal_err_max"])
    return out


def compare_candidates(points: gpd.GeoDataFrame, candidates: dict[str, CRS]) -> list:
    rows = [(name, population_weighted_distortion(points, crs))
            for name, crs in candidates.items()]
    rows.sort(key=lambda r: r[1]["linear_err_weighted_mean"])
    log.info("best by population-weighted linear error: %s", rows[0][0])
    return rows

The standard-parallel rule of thumb — one sixth in from each edge of the latitude span — minimises the maximum scale error for a region of moderate north-south extent, and it is a good starting point rather than a law. Weighting the distortion by population then judges the result where it matters. A region whose population is concentrated in a southern metropolitan corridor may be better served by parallels shifted south, and the comparison function will show it.

Checking areal error as well as linear error is a correctness guard: an equal-area projection should have an areal scale of one everywhere, and any departure indicates the definition is not what was intended.

The Handoff Between State and Regional Systems Three participating states each hold data in their own native projected system and realization. Each state's data is transformed once into the shared regional equal-area system, after first unifying to the agreed datum realization, and all regional analysis happens there. Products for an individual state are transformed back out to that state's native system as a one-way output. The regional system definition is published as a single shared file so every participant uses byte-identical parameters. In once, analyse once, out once — with one published definition state A native state B native state C native regional equal-area all analysis here state products one-way out to native unify the datum realization before entering; publish the regional definition as one shared file

Validation & Edge Cases

A region with a long north-south extent. Albers equal-area handles moderate latitude spans well; a region spanning twenty degrees or more accumulates linear distortion even with optimised parallels. Consider a Lambert azimuthal equal-area projection centred on the region, and compare both with the population-weighted measure.

Participants on different realizations. A shared regional projection defined on one datum realization, fed with data in another, carries the realization shift into every regional analysis. Agree the realization first, as covered in Handling Datum Shifts Between NAD83 Realizations.

Custom definitions that drift between participants. A projection defined by parameters typed into each participant’s software will differ in a decimal somewhere. Publish the definition as a single well-known-text file and have everyone load that file.

Tools that do not accept custom projections. Some web mapping and routing tools accept only registered codes. Keep the regional system for analysis and produce web outputs in a standard web projection as a final, one-way step.

Distance operations that should be geodesic. For very long distances — cross-region routing, air travel — projected distance in any system is an approximation. Compute those geodesically and reserve projected distance for the local operations it is designed for.

Boundary files clipped at state lines. Each participant’s tract file ends at its own state border, and the digitised borders rarely coincide exactly. Once transformed into the regional system, adjacent states’ tracts overlap or leave slivers along the shared line, which breaks contiguity weights and scan windows that cross it. Reconcile the border — snapping both sides to one authoritative state-line geometry — before building any regional weights.

Region membership that changes. A fourth state joining the program extends the region and moves its centre. Re-optimising the projection for the larger region would change every existing coordinate; usually the better choice is to keep the published definition and accept slightly higher distortion in the new state, recording that decision. Measure the new state’s population-weighted error before deciding.

Distance thresholds tuned in one system and applied in another. A scan radius or buffer calibrated in a state plane zone behaves slightly differently in the regional system near its edges. Re-validate distance-based parameters after the move rather than carrying them over, particularly for detection thresholds whose false-alarm rate was tuned empirically.

The border reconciliation step in that list is the one most often skipped, and its effect on contiguity is easy to picture:

Reconcile the State Line Before Building Regional Weights Two adjacent states' tract layers meeting at their shared border. Before reconciliation, each state's digitised border differs slightly, leaving a thin sliver between them in some stretches and an overlap in others; tracts on opposite sides of the sliver are not detected as neighbours, so cross-border contiguity is lost. After both sides are snapped to one authoritative state-line geometry, the tracts share an edge and cross-border neighbours are found, which regional cluster detection depends on. A sliver on the border hides every cross-border neighbour Before state A tract state B tract 8 m sliver: not neighbours After state A tract state B tract shared edge: neighbours 
INFO crs.region regional Albers: lat1=35.217 lat2=40.883 lon0=-86.450 (datum intent NAD83(2011))
INFO crs.region custom: pop-weighted linear error 0.00031, p99 0.00094, max areal error 2.20e-15
INFO crs.region best by population-weighted linear error: regional_albers

Compliance Notes

Record the regional projection’s full definition, the method used to choose its parameters, the population-weighted distortion comparison against the alternatives, the agreed datum realization, and the transformations used to bring each participant’s data in.

In a multi-agency program the projection definition is a shared standard, and it belongs in the data-sharing agreement or its technical annex rather than in one agency’s code. Changing it later moves every coordinate in every participant’s regional archive, so it should be changed only by agreement, versioned like any other shared reference, and accompanied by the transformation that maps old regional coordinates to new ones.