Spatial Epi Public Health GIS

Testing Geospatial Analysis Code

This guide is part of the Surveillance Pipeline Automation & Spatial Data Quality section. It covers the discipline that makes a scheduled pipeline safe to change: tests that catch a wrong answer before it reaches a published map.

Concept & Epidemiological Alignment

Conventional unit testing compares a function’s output to an expected value written by hand. Spatial analysis code resists this in three ways. The correct output of a reprojection, an overlay or a buffer is a geometry with thousands of coordinates that nobody can write down. The correct output of a clustering statistic depends on a random permutation and differs from run to run. And the correct output of many functions depends on library versions — a projection engine update can legitimately shift a coordinate by a metre — so a hard-coded expected value becomes wrong for reasons unrelated to the code under test.

The discipline that works instead rests on four ideas.

Invariants rather than values. Many spatial operations have properties that must hold regardless of input: a reprojection round-trip returns the original within tolerance; an areal interpolation conserves the total; a buffer contains its source geometry; a dissolve preserves total area. These are writable when the expected output is not, and they catch the errors that matter.

Synthetic data with a known answer. Generate case points with a planted cluster at a known location and size, run the detection code, and assert it finds that cluster. The expected answer is known because you put it there.

Pinned dependencies. Projection and geometry engines move answers. Pin them, and treat an upgrade as a deliberate change that must pass a regression comparison.

Regression against a stored result. For a fixed input and fixed seeds, the current output of an end-to-end analysis is compared to a stored output from a known-good version. A difference is not automatically a failure, but it is always a decision.

The epidemiological reason this matters is that the errors spatial code produces are usually plausible. A tract assignment wrong by one polygon, a weight matrix with one missing neighbour, a projection off by the datum shift — none look wrong on a map, and all change a cluster boundary or a rate.

A second reason is institutional. Surveillance code is maintained by rotating staff — a fellow for two years, a contractor for six months — and the test suite is the only durable description of what the code is supposed to do. A pipeline without tests can be run by its author and nobody else; a pipeline with them can be changed by a newcomer who can see, within minutes, whether they broke anything that mattered.

Testing Strategies by What Is Being Tested

What is tested Strategy Example assertion
Coordinate transformations Round-trip invariant Transform and inverse return input within 1 cm
Areal interpolation Conservation invariant Source total equals target total
Buffers and overlays Containment and area invariants Buffer contains source; intersection area ≤ both inputs
Spatial weights Structural properties Symmetric for contiguity; no self-neighbours; islands declared
Cluster statistics Planted signal Detects a planted cluster at a known location
Whole pipeline Regression on fixed input Output matches stored known-good within tolerance
Library upgrades Differential comparison Old and new versions agree on the regression fixture
Three Layers of Geospatial Tests A three-layer arrangement. The broad base is invariant and property tests: many fast tests checking round-trips, conservation, containment and structural properties, run on every commit. The middle layer is planted-signal tests on synthetic data, fewer and slower, checking that statistics find what was put there, run on every commit to analysis code. The narrow top is whole-pipeline regression on a fixed fixture with stored known-good output, run on every deployment and every dependency upgrade. Each layer names what it catches that the one below cannot. Many cheap invariants, fewer planted signals, one end-to-end regression Pipeline regression fixed fixture, stored output every deploy and upgrade Planted-signal tests synthetic cases with a known cluster, known rate, known exposure every change to analysis code Invariant and property tests round-trips, conservation, containment, symmetry, generated inputs every commit, seconds to run catches: composition catches: logic Each layer catches a class the layer below structurally cannot

Spatial Data Prerequisites

  • A fixed test coordinate reference system declared in the test configuration, identical to production.
  • Synthetic fixtures generated from a seed, covering the geometry types and edge cases the code must handle: islands, multipart polygons, points on boundaries, empty inputs.
  • A small real-shaped fixture — for example a few dozen synthetic tracts with realistic topology — for regression tests.
  • Pinned versions of every geometry and projection library, with the projection data version pinned alongside the library.
  • Tolerances declared per assertion, derived from the precision the analysis actually needs rather than from floating-point defaults.

Production Implementation

The tests below illustrate each layer: invariants over generated inputs, a planted-cluster test, and a regression comparison with a declared tolerance.

# Geospatial tests: invariants, a planted signal, and a pipeline regression.
# Pinned: pytest==8.2.2, hypothesis==6.103.1, geopandas==0.14.4, shapely==2.0.4,
#         pyproj==3.6.1, numpy==1.26.4, libpysal==4.10, esda==2.5.1
import numpy as np
import geopandas as gpd
import pandas as pd
import pytest
from hypothesis import given, settings, strategies as st
from pyproj import Transformer
from shapely.geometry import Point, box

WORK_CRS = "EPSG:5070"

# ---------- Layer 1: invariants over generated inputs ----------

lon = st.floats(min_value=-124.0, max_value=-67.0, allow_nan=False)
lat = st.floats(min_value=25.0, max_value=49.0, allow_nan=False)


@given(lon, lat)
@settings(max_examples=500, deadline=None)
def test_reprojection_round_trip(x, y):
    fwd = Transformer.from_crs("EPSG:4326", WORK_CRS, always_xy=True)
    inv = Transformer.from_crs(WORK_CRS, "EPSG:4326", always_xy=True)
    ex, ny = fwd.transform(x, y)
    bx, by = inv.transform(ex, ny)
    # 1e-7 degrees is about a centimetre: far tighter than any epi use needs.
    assert abs(bx - x) < 1e-7 and abs(by - y) < 1e-7


@given(st.floats(min_value=10, max_value=5000))
def test_buffer_contains_source(radius):
    p = Point(1_000_000, 2_000_000)
    assert p.buffer(radius).contains(p)


def test_areal_interpolation_conserves_total(interpolate_fn, source_zones, target_zones):
    """Whatever the method, the total must survive the move between geographies."""
    src_total = float(source_zones["cases"].sum())
    out = interpolate_fn(source_zones, target_zones, value="cases")
    assert out["cases"].sum() == pytest.approx(src_total, rel=1e-9)


# ---------- Layer 2: planted signal ----------

def synthetic_cases(n_background: int, cluster_center, cluster_n: int,
                    cluster_radius: float, extent, seed: int) -> gpd.GeoDataFrame:
    rng = np.random.default_rng(seed)
    minx, miny, maxx, maxy = extent
    bx = rng.uniform(minx, maxx, n_background)
    by = rng.uniform(miny, maxy, n_background)
    r = cluster_radius * np.sqrt(rng.uniform(0, 1, cluster_n))
    t = rng.uniform(0, 2 * np.pi, cluster_n)
    cx = cluster_center[0] + r * np.cos(t)
    cy = cluster_center[1] + r * np.sin(t)
    pts = [Point(x, y) for x, y in zip(np.r_[bx, cx], np.r_[by, cy])]
    return gpd.GeoDataFrame({"is_planted": [False] * n_background + [True] * cluster_n},
                            geometry=pts, crs=WORK_CRS)


def test_hotspot_finds_planted_cluster(detect_fn):
    extent = (0, 0, 50_000, 50_000)
    center = (32_000, 18_000)
    cases = synthetic_cases(2_000, center, 180, 1_500, extent, seed=20260918)
    found = detect_fn(cases)                       # returns cluster polygons
    assert len(found) >= 1, "planted cluster not detected"
    hit = found.geometry.contains(Point(center)).any()
    assert hit, "detected cluster does not contain the planted centre"


def test_hotspot_quiet_on_null(detect_fn):
    """Equally important: no cluster when none was planted (at the stated alpha)."""
    extent = (0, 0, 50_000, 50_000)
    cases = synthetic_cases(2_180, (0, 0), 0, 1, extent, seed=20260919)
    assert len(detect_fn(cases)) == 0


# ---------- Layer 3: pipeline regression ----------

def test_pipeline_regression(run_pipeline, fixture_dir, stored_output):
    """Fixed fixture, fixed seeds: output must match the stored known-good."""
    out = run_pipeline(fixture_dir, seed=20260918)
    ref = pd.read_parquet(stored_output)
    merged = ref.merge(out, on="geoid", suffixes=("_ref", "_new"), how="outer")
    assert merged["cases_ref"].equals(merged["cases_new"]), "case counts changed"
    assert np.allclose(merged["rate_ref"], merged["rate_new"], atol=1e-9, equal_nan=True)
    assert (merged["cluster_ref"] == merged["cluster_new"]).all(), "cluster membership changed"

The null test beside the planted-cluster test is the half most test suites omit. A detector that finds the planted cluster and also finds three others in pure noise passes the first test and is useless. Testing both directions — sensitivity on a planted signal, specificity on a null — is what makes the pair meaningful.

The regression comparison treats counts as exact and rates as floating-point with a tight tolerance, and cluster membership as exact. That split reflects what can legitimately vary between runs: nothing should change a count, arithmetic ordering can change the last bit of a rate, and a change in cluster membership is always a finding.

Parameter Selection & Tuning

Tolerances from requirements, not defaults. A centimetre tolerance on a reprojection round-trip is appropriate; a millimetre tolerance fails on legitimate engine differences, and a metre tolerance misses a datum error. Derive each tolerance from the precision the analysis needs.

Number of generated examples. A few hundred examples per invariant run in seconds and find most edge cases. Property-based frameworks shrink failures to a minimal example, which is often the most useful debugging artifact the suite produces.

Planted-signal strength. A planted cluster so strong that any method finds it tests nothing. Plant signals at the effect sizes the pipeline is expected to detect in practice, and include one just below the detection limit to confirm the method’s sensitivity is where you think it is.

Regression fixture size. Large enough to exercise every code path — islands, multipart features, empty areas, suppressed cells — and small enough to run in a minute. A few dozen areas and a few thousand synthetic cases usually suffice.

Seeds. Every stochastic test uses a fixed seed. A test that passes on some seeds and fails on others is reporting a real property of the method, and should be rewritten to test that property explicitly across a set of seeds.

Edge Cases & Failure Modes

Tests that pass because the fixture is too easy. Synthetic tracts that are all squares never exercise the island, multipart and sliver cases where spatial code actually breaks. Include deliberately awkward geometry in every fixture.

Stored regression output updated to make a test pass. When a regression test fails, updating the stored output is sometimes correct and always a decision. Require the update to be a separate, reviewed change with a stated reason.

Projection data outside version control. The projection library’s grid files and database can change independently of the library version. Pin the data version alongside the library, as covered in Pinning PROJ and GDAL Versions for Reproducible Results.

Tests that need protected data. A test suite that only runs inside the protected environment is a test suite most developers cannot run. Synthetic fixtures exist precisely so that the tests run everywhere.

Flaky permutation-based tests. A significance test using permutation inference varies with the permutation count and seed. Fix both, and choose planted effects comfortably above or below the threshold rather than at it.

Sensitivity and Specificity Need Separate Tests Two synthetic study areas. The left one contains uniform background cases plus a planted cluster of known centre and radius; the detector's output is overlaid and correctly contains the planted centre. The right one contains uniform background cases only, with the same total count; the detector correctly returns no cluster. A detector that passed the left test and failed the right would find clusters everywhere, and a suite containing only the left test would not notice. A detector must be tested for what it finds and for what it does not Planted cluster found, contains the centre: pass Null: no cluster planted nothing found: pass Only the pair together shows the detector distinguishes signal from noise

The fixture-difficulty point in that list is easiest to see by listing which code paths an easy fixture never reaches:

An Easy Fixture Leaves the Risky Paths Untested Six code paths in a spatial pipeline listed down the side: ordinary contiguity, island handling in spatial weights, multipart feature handling, points exactly on a boundary, sliver detection, and empty-area aggregation. An easy fixture of regular square tracts exercises only ordinary contiguity. An awkward fixture that deliberately includes an island, a multipart tract, boundary points, a sliver and an empty area exercises all six. The paths the easy fixture misses are precisely the ones where real boundary files break code. Put the awkward cases in the fixture on purpose square tracts awkward fixture ordinary contiguity tested tested island in weights never tested multipart feature never tested point exactly on a boundary never tested sliver between neighbours never tested area with no cases never tested add a new awkward case to the generator each time a real defect is found

Compliance & Audit Controls

The test suite is evidence. Record, for every deployment, the commit tested, the library and projection-data versions, the test results, and — for the regression layer — whether the stored output was updated and why. An agency asked to demonstrate that its published cluster alerts come from a validated method can point to that record.

Synthetic fixtures must be genuinely synthetic. A fixture made by jittering real case locations is a masked release, and it will be committed to a code repository with far wider access than the case file. Generate fixtures from distributions, and review any fixture derived from real data under the disclosure process.

Production Implementation Checklist

Frequently Asked Questions

How do I test a function whose output is a complicated geometry? Test its properties instead of its coordinates. A buffer contains its input; an intersection is contained in both operands; a union’s area is at least each operand’s; a reprojection round-trip returns the input. These assertions are easy to write and catch the realistic errors.

Is property-based testing worth the setup for spatial code? Yes, particularly for coordinate handling. Generated inputs find the antimeridian, the pole, the exactly-on-boundary point and the degenerate polygon far more reliably than hand-written cases, and the shrinking behaviour hands you a minimal failing example.

What should a regression test compare? Counts exactly, continuous values with a tight tolerance, and categorical outcomes — cluster membership, suppression flags — exactly. A change in any of the exact comparisons is a decision someone has to make deliberately.

How do I test code that only runs on protected data? Refactor it so the protected-data dependency is at the edge and the logic takes ordinary inputs, then test the logic on synthetic fixtures. The geocoding and exposure joins are the usual culprits, and they are also the stages most in need of tests.

When should the stored regression output be updated? When a change in output is intended and understood — a bug fix, a method improvement, a deliberate dependency upgrade — and never to make a failing test pass without that understanding. The update should be reviewed separately from the code change that caused it.