Spatial Epidemiology Fundamentals & Data Standards

Production spatial epidemiology fails not at the statistics but at the data layer — an undefined datum, a silently truncated attribute, or a census-tract revision that manufactures a “cluster” out of a denominator change. This section establishes the engineering standards a Python/GIS surveillance pipeline must enforce so that every case rate, exposure surface, and cluster declaration it emits is reproducible, privacy-compliant, and defensible under regulatory audit.

The reference architecture below shows how raw feeds become defensible intelligence through sequential validation gates:

Each gate in this diagram is a hard stop: a record that fails it is quarantined, not passed downstream with a warning. The four implementation areas that own these gates each have a dedicated guide — Coordinate Reference Systems for Public Health, Spatial Data Types & Formats, Precision Standards in Epi-Mapping, and Compliance Mapping Frameworks — and the sections below walk the pipeline in execution order, linking each to the standard that governs it.

Data Governance & Compliance Architecture

Public health GIS operates inside a lattice of overlapping legal constraints: the HIPAA Safe Harbor and Expert Determination de-identification standards at the federal level, GDPR Article 9 special-category protections for any EU-resident data, and state-level health information privacy statutes that frequently impose stricter small-count suppression thresholds than HIPAA alone. Treating these as a post-hoc reporting filter is the most common architectural error — by the time data reaches the reporting tier, raw coordinates have already been written to intermediate tables, logs, and spatial indexes. Compliance must therefore be an ingestion-layer property, enforced before any patient-level geometry is persisted at full precision. The full statutory mapping, suppression-threshold logic, and lineage requirements are covered in Compliance Mapping Frameworks; the architectural commitments below are the ones every pipeline in this section assumes.

De-identification at ingestion means three things run before the first spatial join. Direct identifiers (name, MRN, full date of birth) are hashed or dropped; small-count cells are suppressed against a documented threshold so that no published areal unit can be back-solved to an individual; and geocoded coordinates that fall below an analytical confidence interval are masked or aggregated rather than stored rooftop-precise. The suppression threshold itself is a logged parameter, not a constant buried in code — auditors must be able to read it from the run record.

Audit logging is the other half of governance. Every pipeline stage emits structured logs capturing input/output record counts, validation failure rates, and a SHA-256 checksum of the input artifact, so that any published result can be traced back to the exact bytes it was computed from. Immutable lineage tracking records transformation timestamps, library and algorithm versions, and operator credentials. This is what makes a result defensible: a reviewer can re-run the documented configuration and obtain byte-identical output. The following pattern hashes the input and stamps an ISO 19115-style lineage record at the moment of ingestion:

# Versions: geopandas==0.14.4, pyproj==3.6.1, pandas==2.2.2 (Python 3.11)
import hashlib
import json
import logging
from datetime import datetime, timezone

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

def attest_ingest(path: str, suppression_threshold: int, operator: str) -> dict:
    """Hash the source artifact and emit an immutable lineage record at ingestion."""
    with open(path, "rb") as fh:
        digest = hashlib.sha256(fh.read()).hexdigest()
    record = {
        "source_path": path,
        "sha256": digest,
        "ingested_utc": datetime.now(timezone.utc).isoformat(),
        "suppression_threshold": suppression_threshold,  # logged, not hard-coded silently
        "operator": operator,
        "schema_version": "epi-ingest/1.4.0",
    }
    logging.info("INGEST lineage: %s", json.dumps(record))
    return record  # persist to an append-only provenance store

Spatial Data Preparation

Once an artifact is attested, preparation enforces the geometric preconditions that every downstream statistic silently assumes. The first is coordinate reference governance. Projection mismatch is the primary failure vector in multi-source epidemiological mapping: surveillance points, environmental rasters, and administrative boundaries routinely arrive with conflicting datums or no defined CRS at all. Distance and adjacency computed on unprojected WGS84 degrees are geometrically meaningless, and the distortion grows with latitude — a fixed-distance band that behaves at 30°N silently mis-neighbours features at 60°N. Pipelines must declare a single canonical CRS for analytical operations, reject geometries lacking a valid EPSG code, and choose a projection by analytical intent: equal-area projections preserve the rate denominators in continental incidence modeling, while localized outbreak investigations demand high-precision UTM or State Plane zones. The selection rationale, datum-shift handling, and tolerance checks are detailed in Coordinate Reference Systems for Public Health; the official pyproj documentation is the authoritative reference for transformation-pipeline configuration.

The second precondition is topological validity. Self-intersecting polygons, unclosed rings, and mixed geometry types cause geopandas operations to fail or, worse, return wrong areas without raising. Validation gates must assert validity and heal slivers before any weights matrix or buffer is constructed, routing irreparable geometries to the quarantine table rather than dropping them silently:

# Versions: geopandas==0.14.4, shapely==2.0.4, pyproj==3.6.1
import geopandas as gpd
from shapely.validation import make_valid
import logging

CANONICAL_CRS = "EPSG:5070"  # CONUS Albers equal-area; rate-preserving for incidence

def prepare(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Enforce CRS, repair topology, and quarantine irreparable geometry."""
    if gdf.crs is None:
        raise ValueError("Geometry has no CRS; reject to prevent silent misalignment.")
    gdf = gdf.to_crs(CANONICAL_CRS)
    bad = ~gdf.geometry.is_valid
    if bad.any():
        logging.warning("Repairing %d invalid geometries", int(bad.sum()))
        gdf.loc[bad, "geometry"] = gdf.loc[bad, "geometry"].apply(make_valid)
    still_bad = ~gdf.geometry.is_valid
    quarantine = gdf[still_bad].assign(error_code="TOPOLOGY_IRREPARABLE")
    return gdf[~still_bad], quarantine  # valid frame, quarantined frame

Schema validation belongs in this same layer. Mandatory fields — case identifier, onset date, diagnostic code, geocoded coordinates — are enforced declaratively at parse time with pydantic or great_expectations, and records that fail are quarantined with a structured error code rather than propagating into spatial joins or kernel density estimates. Validation failure rate is itself a monitored metric: a sudden spike usually signals an upstream feed schema change, not bad luck.

Core Methods Overview

Prepared, governed geometry feeds the detection layer, where the correct method is dictated by data geometry rather than analyst preference. The full implementations live in the Disease Clustering & Spatial Statistical Modeling section; the decision order below is what routes a dataset to the right one.

Aggregated counts per administrative unit — cases per census tract, normalized to a population denominator — go first to global autocorrelation. Global & Local Moran’s I Implementation answers a single yes/no question about whether spatial structure exists at all before any local map is drawn. The global statistic is

$$I = \frac{n}{\sum_{i}\sum_{j} w_{ij}} \cdot \frac{\sum_{i}\sum_{j} w_{ij}(x_i - \bar{x})(x_j - \bar{x})}{\sum_{i}(x_i - \bar{x})^2}$$

where $w_{ij}$ is the row-standardized spatial weight between units $i$ and $j$. Only when global $I$ is significant does the local decomposition (LISA) earn the right to flag individual units.

When the operational need is to zone an intervention — to draw the boundary of a hot area for resource targeting — Getis-Ord Gi* Hotspot Detection produces a signed z-score surface in which high-value and low-value clusters are distinguishable. Its statistic is

$$G_i^{*} = \frac{\sum_{j} w_{ij} x_j - \bar{x}\sum_{j} w_{ij}}{s\sqrt{\dfrac{n\sum_{j} w_{ij}^2 - \left(\sum_{j} w_{ij}\right)^2}{n-1}}}$$

For precise event points rather than areal aggregates — geocoded vector-trap captures or individual case residences where the events are the data — K-Function & Point Pattern Analysis measures clustering as a function of distance scale via Ripley’s estimator,

$$\hat{K}(r) = \frac{|A|}{n(n-1)} \sum_{i \neq j} \mathbf{1}\{d_{ij} \le r\}\, e_{ij}$$

where $e_{ij}$ is an edge-correction weight. Finally, when where and when must be answered together — a timestamped event stream in active surveillance — Spatial Scan Statistics Configuration applies a likelihood-ratio scan across overlapping space-time windows. Every one of these methods consumes the same adjacency object, so the choices made in Spatial Weights Matrix Construction — contiguity versus distance versus k-nearest-neighbor, and whether the matrix is row-standardized — propagate into every statistic downstream. Routing every dataset through this ordered decision keeps the entire section operating on one consistent weights object and one consistent CRS.

Threshold Tuning & Validation

Local statistics test one hypothesis per spatial unit, so a raw $\alpha = 0.05$ cutoff applied across a thousand units yields roughly fifty false hotspots by construction. Every method above must therefore correct for multiplicity. The recommended control is the Benjamini–Hochberg false-discovery-rate procedure, which bounds the expected proportion of false positives among declared clusters while retaining far more power than Bonferroni at surveillance scale:

# Versions: numpy==1.26.4
import numpy as np

def bh_fdr(pvals: np.ndarray, alpha: float = 0.05) -> np.ndarray:
    """Benjamini-Hochberg FDR mask over local p-values; returns survivors."""
    p = np.asarray(pvals, dtype=float)
    n = p.size
    order = np.argsort(p)
    thresh = (np.arange(1, n + 1) / n) * alpha
    passed = p[order] <= thresh
    k = np.max(np.where(passed)[0]) + 1 if passed.any() else 0
    keep = np.zeros(n, dtype=bool)
    if k:
        keep[order[:k]] = True
    return keep  # boolean mask of units surviving FDR control at alpha

Per-run correction is necessary but not sufficient. Drift detection monitors shifts in the denominator population and covariate distributions between releases — a quietly revised census tract will manufacture a “signal” that is really a reclassification artifact, the same failure mode that motivates the versioned-boundary registry below. Cross-validation against historical baselines separates a genuine emerging signal from recurring seasonal structure, and in near-real-time architectures, reporting-delay correction through nowcasting and adaptive windowing prevents both premature declarations and missed onsets. The alpha level, FDR method, drift thresholds, and permutation seed are all logged parameters, because a clustering declaration is only defensible if a reviewer can reproduce the exact tuning that produced it.

Where These Pipelines Actually Break

The order above describes how a governed pipeline is meant to run. It is worth knowing, separately, where such pipelines fail once they are running unattended — because the distribution of real incidents is not the distribution a first-time reader would guess. The stages that consume the most engineering attention are rarely the ones that generate the most re-runs.

Three things follow from that shape, and they are the reason this section is ordered the way it is.

Preparation failures outnumber analysis failures roughly four to one. No statistic in the catalogue — not Moran’s I, not a scan window, not a 2SFCA decay function — appears in the top six causes. That is not because the statistics are simple; it is because a wrong statistic usually produces a suspicious number that somebody questions, whereas a wrong coordinate reference produces a confident number that nobody does. The asymmetry is worth stating plainly: the cheapest place to catch an error is the stage where the error is still visible as a type violation rather than as an epidemiological claim.

The top two causes are both mismatches between two datasets, not defects in one. A datum mismatch requires a case layer and a boundary layer that disagree; a match-rate problem requires an address file and a reference range that disagree. Single-file validation — is this geometry valid, does this file parse — catches neither. The gates that matter are therefore relational: they compare the authority code on one frame against the authority code on the other, and they compare the geocoded distribution against the population distribution it should resemble. Any validation suite built only from per-file assertions will pass a pipeline that is about to produce a systematically displaced map.

Boundary vintage drift is the slow one. It sits fourth by count but first by time-to-detect, because nothing about a 2020 tract boundary joined to a 2010 population denominator throws an error — the join succeeds, the rates come out, and the trend line bends for a reason nobody can name until someone plots the denominators. Pin the vintage of every boundary file in the same place you pin the CRS, and make a vintage mismatch as loud as a CRS mismatch, because epidemiologically it is the same class of mistake: two layers that describe different worlds being treated as though they describe one.

The remainder of this section takes those failure modes in the order a pipeline meets them. Coordinate handling and precision come first because they are upstream of everything and because their failures are silent. Formats and types come next, since a format decision fixes which metadata can survive at all. Weights construction comes last of the preparation topics because it is the first place the geometry stops being a map and becomes a model — and the point at which an undetected island quietly changes the sample size every later statistic reports.

Operationalization & Output Standards

The final stage serializes results for handoff, and format choice is itself a compliance decision. Legacy shapefiles introduce attribute truncation, encoding inconsistency, and topology errors that compromise reproducibility; production pipelines deprecate them. Prefer GeoPackage for transactional vector storage, GeoParquet for analytical handoff — it preserves column types, CRS, and bounding box, and reads column-selectively at scale — and Cloud-Optimized GeoTIFF for environmental covariate surfaces. The format-by-format selection criteria, spatial indexing, and conversion patterns are detailed in Spatial Data Types & Formats, with the OGC GeoPackage specification as the authoritative baseline.

Every serialized artifact carries the same provenance block recorded at ingestion: the CRS authority code, an ISO 19115 lineage record, and the SHA-256 configuration hash, so the audit loop opened at ingestion closes at output. Positional accuracy metadata travels with the geometry as well — geocoder confidence, the precision tier, and any uncertainty buffer — following Precision Standards in Epi-Mapping so that downstream consumers can propagate locational error into regression and exposure models rather than treating coordinates as exact.

Administrative boundaries are not static, which is why output must be stamped with the boundary vintage it was computed against. Census tracts, voting districts, and health service areas undergo periodic revision, so a time series spanning multiple vintages requires deterministic alignment — areal interpolation or dasymetric mapping to normalize historical counts onto a contemporary geography — before aggregation. A versioned boundary registry and automated drift detection for slivers, gaps, and overlaps guarantee longitudinal consistency. When a confirmed surface is handed to a downstream accessibility analysis — pairing a hotspot with Healthcare Access & Network Analysis Automation to size the response — the shared provenance metadata is what lets the two pipelines be reconciled after the fact.

The matrix below scores the four serialization targets a pipeline actually emits against the properties that decide which one a given handoff requires:

Where this section ends and modelling begins

The boundary is worth naming precisely, because work that belongs on one side of it is routinely attempted on the other. Everything covered here is about making two or more datasets describe the same world: the same coordinate reference, the same vintage, the same units, the same precision, the same neighbour relations. None of it involves a hypothesis. The moment a statistic is computed that could come out one way or another, the work has moved into the modelling sections, and the inputs should already be frozen.

Keeping that boundary sharp has a practical payoff. When a result looks wrong, the first question is always whether it is a modelling problem or a preparation problem, and the answer is much faster to find when preparation has its own gates, its own logs, and its own pass or fail. A pipeline that interleaves the two — reprojecting inside the loop that fits the model, say — makes every wrong answer expensive to diagnose.

Production Implementation Checklist

• Source artifact hashed (SHA-256) and an immutable lineage record written before any geometry is persisted
• Direct identifiers hashed or dropped and small-count cells suppressed at ingestion, with the suppression threshold logged
• Mandatory fields enforced declaratively at parse time; failed records quarantined with a structured error code
• Every layer reprojected to a single documented equal-area or zone CRS; geometries lacking a valid EPSG code rejected
• Geometry validity asserted and slivers healed before any weights matrix or buffer is built; irreparable geometry quarantined
• Projection chosen by analytical intent (equal-area for rates, UTM/State Plane for local precision) and recorded
• Positional accuracy metadata and uncertainty buffers attached to all geocoded outputs
• Local p-values corrected with Benjamini–Hochberg FDR at a documented alpha, not a raw cutoff
• Drift checks run on denominator population and covariate distributions between releases
• Administrative boundaries version-controlled; historical counts areally interpolated to a single vintage before aggregation
• Output written as GeoPackage / GeoParquet / COG with CRS, ISO 19115 lineage, and the configuration hash embedded
• Pipelines containerized with pinned dependency manifests so results are byte-reproducible

Related Topics

• Coordinate Reference Systems for Public Health — projection selection and datum governance the entire pipeline depends on.
• Spatial Data Types & Formats — vector, raster, and columnar serialization standards for audit-ready handoff.
• Precision Standards in Epi-Mapping — geocoder accuracy tiers and uncertainty propagation into spatial models.
• Compliance Mapping Frameworks — HIPAA, GDPR, and state-law constraints encoded as ingestion-layer gates.
• Disease Clustering & Spatial Statistical Modeling — the detection methods this prepared, governed data feeds.