Spatial Data Testing & CI Pipelines

Problem Framing

A routing service that plans evacuation corridors ships a change to its buffer logic on a Friday afternoon during pre-landfall staging. The pull request passes review — the diff is three lines — and merges without a spatial test in the way. Over the weekend the change silently switches the buffer distance from metres to degrees because an upstream layer arrived in a geographic reference system instead of the projected one the function assumed. The unit that consumes the buffer now treats a 300-metre standoff as a 300-degree one, and by the time an analyst notices the evacuation zone swallowing three counties, the operational map has already been pushed to field tablets. No exception was raised; every function returned a valid geometry. This is the defining hazard of spatial software: wrong answers are still well-formed geometries, and a system that only tests “did it return a shape” will happily serve catastrophically wrong shapes. This page specifies the testing and continuous-integration discipline that makes that failure fail the build instead of the incident, implementing the Python Toolchains for Public Safety GIS contract for verifiable, reproducible spatial code under National Incident Management System (NIMS) and Federal Emergency Management Agency (FEMA) reporting obligations.

Prerequisites

This workflow assumes a senior engineer’s fluency with pytest and the Python geospatial stack, plus the following preconditions before the first test runs:

• Packages: pytest >= 7.4, geopandas >= 0.14, shapely >= 2.0, pyproj >= 3.6, and hypothesis >= 6.90 for property-based coverage. Shapely 2.x is required so that vectorized predicates and shapely.validation.make_valid behave consistently across the suite.
• A pinned toolchain: GDAL, PROJ, GEOS, and the PROJ data directory must be version-locked in the test environment. The reproducible container that provides this is owned by the Setting Up Dockerized GIS Environments workflow; the test suite assumes it is running inside that image so datum grids and topology resolve identically on every machine.
• A declared schema contract: the GeoDataFrames under test must map to an explicit column-and-dtype contract (required attributes, bounded domains, an expected geometry type, and a canonical EPSG code). Tests assert against that contract rather than against whatever a sample file happened to contain.
• Version-controlled fixtures: golden files and calibration parameters live under the same review and history discipline as code, per the Version Control for Spatial Workflows pattern, so that a change to expected output is a reviewable diff and not a silent overwrite.

Test Architecture

Safety-critical spatial code needs a layered test strategy because different classes of defect surface at different granularities. A malformed buffer is a unit-level topology bug; a datum shift that only appears near the antimeridian is a property-level bug; an ingest pipeline that reprojects twice is an integration bug. Organize coverage as a pyramid: a wide base of fast, dependency-free geometry and topology assertions; a layer of schema and contract tests that pin the data shape; property-based transform tests that generalize beyond hand-picked examples; golden-file regression that locks whole-output behaviour; and a thin top of end-to-end tests that exercise the real ingest-to-publish path. The entire pyramid runs inside a pinned GDAL and PROJ container on every push, and the continuous-integration gate fails closed — a single failing spatial assertion blocks the merge rather than emitting a warning that a reviewer can wave through.

Step-by-Step Implementation

Step 1 — Build deterministic GeoDataFrame fixtures

Tests that read a shared production extract are not tests; they are a slow, non-deterministic dependency on whatever that file contains today. Build small, self-describing fixtures from inline Well-Known Text so every test starts from a known spatial state, tagged with an explicit reference system. Keep them tiny — three or four features is enough to exercise a predicate — so the base of the pyramid stays fast.

import logging
from typing import Dict

import geopandas as gpd
import pytest
from shapely import wkt
from shapely.geometry.base import BaseGeometry

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

# Canonical projected CRS for the operational area (UTM zone 10N / WGS 84).
FIXTURE_CRS = "EPSG:32610"


@pytest.fixture
def incident_zones() -> gpd.GeoDataFrame:
    """Three known incident polygons in a metric CRS for topology tests."""
    features: Dict[str, str] = {
        "staging": "POLYGON ((500000 4200000, 500300 4200000, 500300 4200300, 500000 4200300, 500000 4200000))",
        "hot_zone": "POLYGON ((500250 4200100, 500600 4200100, 500600 4200450, 500250 4200450, 500250 4200100))",
        "detached": "POLYGON ((501000 4201000, 501200 4201000, 501200 4201200, 501000 4201200, 501000 4201000))",
    }
    try:
        geoms: list[BaseGeometry] = [wkt.loads(w) for w in features.values()]
    except Exception as exc:  # shapely raises GEOSException subclasses on bad WKT
        logger.error("Fixture WKT failed to parse: %s", exc)
        raise
    gdf = gpd.GeoDataFrame(
        {"zone_id": list(features.keys()), "priority": [1, 1, 3]},
        geometry=geoms,
        crs=FIXTURE_CRS,
    )
    logger.info("Built incident_zones fixture with %d features in %s", len(gdf), gdf.crs)
    return gdf

Because the fixture declares its CRS and its attributes inline, any test that consumes it is testing against a contract, not against a sample. When a bug report arrives, the fixture that reproduces it becomes a permanent regression case with no external file to lose.

Step 2 — Assert geometry and topology invariants

The core defense against well-formed-but-wrong output is a set of assertion helpers that check the properties a spatial operation must preserve: validity, the expected reference system, the expected geometry type, and topological relationships between features. Encode these as reusable functions so every test reads as a statement of intent rather than a wall of .iloc indexing.

import logging

import geopandas as gpd
from shapely.geometry.base import BaseGeometry
from shapely.validation import explain_validity

logger = logging.getLogger(__name__)


class SpatialAssertionError(AssertionError):
    """Raised when a geometry or topology invariant is violated."""


def assert_all_valid(gdf: gpd.GeoDataFrame) -> None:
    """Every geometry must be non-empty and OGC-valid."""
    invalid = gdf[~gdf.geometry.is_valid | gdf.geometry.is_empty]
    if not invalid.empty:
        reasons = {int(i): explain_validity(g) for i, g in invalid.geometry.items()}
        logger.error("Invalid geometries at rows %s", reasons)
        raise SpatialAssertionError(f"{len(invalid)} invalid/empty geometries: {reasons}")


def assert_crs(gdf: gpd.GeoDataFrame, expected_epsg: int) -> None:
    """The layer must carry the exact expected EPSG code, never None."""
    if gdf.crs is None:
        raise SpatialAssertionError("Layer has no CRS; refusing to trust its coordinates")
    actual = gdf.crs.to_epsg()
    if actual != expected_epsg:
        raise SpatialAssertionError(f"CRS mismatch: expected EPSG:{expected_epsg}, got EPSG:{actual}")


def assert_disjoint(a: BaseGeometry, b: BaseGeometry) -> None:
    """Two features that must never overlap (e.g. mutually exclusive zones)."""
    if a.intersects(b) and a.intersection(b).area > 1e-9:
        raise SpatialAssertionError("Geometries overlap but were required to be disjoint")


def test_buffer_preserves_validity_and_crs(incident_zones: gpd.GeoDataFrame) -> None:
    buffered = incident_zones.copy()
    buffered["geometry"] = buffered.geometry.buffer(50.0)  # 50 metres, metric CRS
    assert_all_valid(buffered)
    assert_crs(buffered, 32610)
    # A 50 m buffer must strictly enlarge every polygon.
    assert (buffered.geometry.area > incident_zones.geometry.area).all()

The final assertion in the test catches exactly the metres-versus-degrees hazard from the problem scenario: if the buffer were applied in a geographic CRS, the area comparison and the assert_crs guard would fail together, and the build would go red before the change reached review.

Step 3 — Lock behaviour with golden-file regression

Unit assertions catch violations of stated invariants, but they cannot catch a subtle drift in a whole processed layer — a reclassification that flips two features, a dissolve that merges one polygon too many. Golden-file regression captures a canonical expected output once, reviews it as a diff, and then compares every subsequent run against it with a geometry-aware tolerance so floating-point noise does not cause false failures.

import logging
from pathlib import Path

import geopandas as gpd
from shapely import equals_exact

logger = logging.getLogger(__name__)

GOLDEN_DIR = Path(__file__).parent / "golden"


def load_or_write_golden(gdf: gpd.GeoDataFrame, name: str, update: bool) -> gpd.GeoDataFrame:
    """Return the stored golden layer, or (re)write it when explicitly updating."""
    path = GOLDEN_DIR / f"{name}.geojson"
    if update or not path.exists():
        GOLDEN_DIR.mkdir(parents=True, exist_ok=True)
        gdf.to_file(path, driver="GeoJSON")
        logger.warning("Wrote golden file %s (%d features)", path, len(gdf))
        return gdf
    try:
        return gpd.read_file(path)
    except Exception as exc:  # fiona/pyogrio raise on unreadable golden
        logger.error("Failed to read golden file %s: %s", path, exc)
        raise


def assert_matches_golden(actual: gpd.GeoDataFrame, name: str, tolerance: float = 1e-6) -> None:
    """Diff a result against its golden layer, tolerant to sub-tolerance coordinate noise."""
    golden = load_or_write_golden(actual, name, update=False)
    if len(actual) != len(golden):
        raise AssertionError(f"Feature count drift: golden={len(golden)}, actual={len(actual)}")
    # Align on a stable key so row order cannot cause spurious diffs.
    a = actual.sort_values("zone_id").reset_index(drop=True)
    g = golden.sort_values("zone_id").reset_index(drop=True)
    mismatches = [
        row.zone_id
        for row, gold in zip(a.itertuples(), g.geometry)
        if not equals_exact(row.geometry, gold, tolerance=tolerance)
    ]
    if mismatches:
        raise AssertionError(f"Geometry drift beyond tolerance for zones: {mismatches}")
    logger.info("Golden match for %s within tolerance %g", name, tolerance)

Updating a golden file is a deliberate act gated behind an update flag (wired to a --update-golden CLI option or an environment variable), never an automatic overwrite. That way a changed golden layer appears in the pull request as a reviewable geometry diff, and an unexplained change to the expected output is caught by a human before it is caught by an incident.

Step 4 — Generalize transforms with property-based tests

Hand-picked examples test the coordinates you thought of; coordinate-transform bugs live in the coordinates you did not — near the poles, across the antimeridian, at the axis-order boundary between geographic and projected systems. Property-based testing with hypothesis generates inputs across the whole valid domain and asserts invariants that must hold for all of them, such as reproject-there-and-back round-tripping within tolerance. The deep treatment of this technique lives in Writing Property-Based Tests for Coordinate Transforms; the pattern below is the load-bearing core.

import logging

from hypothesis import given, settings
from hypothesis import strategies as st
from pyproj import Transformer
from pyproj.exceptions import ProjError

logger = logging.getLogger(__name__)

# always_xy=True forces (lon, lat) ordering on both sides, eliminating axis-order ambiguity.
_FWD = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
_INV = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True)

# Web Mercator is only defined to about +/- 85.06 degrees latitude.
_lon = st.floats(min_value=-179.9, max_value=179.9, allow_nan=False, allow_infinity=False)
_lat = st.floats(min_value=-85.0, max_value=85.0, allow_nan=False, allow_infinity=False)


@settings(max_examples=400, deadline=None)
@given(lon=_lon, lat=_lat)
def test_reproject_round_trip_is_stable(lon: float, lat: float) -> None:
    """Round-tripping WGS 84 -> Web Mercator -> WGS 84 must return the origin within tolerance."""
    try:
        x, y = _FWD.transform(lon, lat)
        lon2, lat2 = _INV.transform(x, y)
    except ProjError as exc:
        logger.error("Transform failed for (%s, %s): %s", lon, lat, exc)
        raise
    # 1e-6 degrees is roughly 0.1 m at the equator: well below operational tolerance.
    assert abs(lon2 - lon) < 1e-6, f"lon drift {abs(lon2 - lon)} at ({lon}, {lat})"
    assert abs(lat2 - lat) < 1e-6, f"lat drift {abs(lat2 - lat)} at ({lon}, {lat})"

When this test fails, hypothesis shrinks the counterexample to the smallest coordinate that breaks the invariant and prints it, so a datum or axis-order regression arrives as a concrete (lon, lat) pair rather than a vague suspicion. Schema and contract tests belong at the same tier: assert the required columns, dtypes, and value domains of every layer the code produces so a renamed attribute or a widened priority range fails immediately.

Step 5 — Gate merges in CI on a pinned toolchain

A green suite on a developer laptop proves nothing if the CI runner resolves datum shifts differently. Run the whole pyramid inside the pinned GDAL and PROJ container on every push and pull request, and configure the job so that any failure blocks the merge. Fail closed: a transform error, an invalid golden read, or a single red assertion must stop the pipeline, never degrade to a warning.

# .github/workflows/spatial-tests.yml
name: spatial-tests
on: [push, pull_request]

jobs:
  test:
    runs-on: ubuntu-latest
    # Pinned image owns exact GDAL, PROJ, GEOS, and PROJ-data versions.
    container: ghcr.io/incidentgis/gis-test:gdal-3.9.2-proj-9.4.1
    steps:
      - uses: actions/checkout@v4
      - name: Verify pinned geospatial stack
        run: |
          python -c "import pyproj, shapely; print('PROJ', pyproj.proj_version_str, 'GEOS', shapely.geos_version_string)"
      - name: Run spatial test suite
        run: pytest -q --maxfail=1 --hypothesis-seed=0

Pinning the --hypothesis-seed makes property-based runs reproducible in CI so a failure can be replayed exactly, while the container tag freezes the datum grids. Building and versioning that image is the responsibility of the Dockerized GIS environment workflow; this job merely consumes a known-good tag and refuses to merge when the suite is red.

Configuration Reference

Tune these knobs per repository; a fast-moving analysis library and a stable ingest service will land on different tolerances and example counts.

Verification & Smoke Test

Before promoting a change to the test harness itself, confirm the assertion helpers actually reject the failures they claim to catch. A test suite that never fails is indistinguishable from no suite at all, so the smoke test deliberately feeds bad input and asserts that the guard fires.

import logging

import geopandas as gpd
from shapely import wkt

logger = logging.getLogger(__name__)


def smoke_test() -> None:
    # 1. assert_crs must reject a layer with no CRS.
    naked = gpd.GeoDataFrame({"zone_id": ["x"]}, geometry=[wkt.loads("POINT (0 0)")], crs=None)
    try:
        assert_crs(naked, 32610)
        raise AssertionError("expected SpatialAssertionError for missing CRS")
    except SpatialAssertionError:
        pass

    # 2. assert_all_valid must reject a self-intersecting bowtie polygon.
    bowtie = gpd.GeoDataFrame(
        {"zone_id": ["b"]},
        geometry=[wkt.loads("POLYGON ((0 0, 1 1, 1 0, 0 1, 0 0))")],
        crs="EPSG:32610",
    )
    try:
        assert_all_valid(bowtie)
        raise AssertionError("expected SpatialAssertionError for invalid geometry")
    except SpatialAssertionError:
        pass

    logger.info("harness smoke test passed: guards reject bad input")


smoke_test()

A CLI equivalent confirms the pinned stack is present and the suite collects without import errors before the full run executes in continuous integration:

python -c "import geopandas, shapely, pyproj, hypothesis; print('stack ok')"
pytest --collect-only -q     # non-zero exit if any test module fails to import

Integration With Adjacent Workflows

This testing discipline is the gate every other Python workflow on the site passes through before its output is trusted. The golden files and calibration parameters it depends on are managed under the spatial version-control contract, so a re-baselined expected layer is a reviewable diff with full history rather than a silent overwrite. The pinned GDAL and PROJ container the CI job runs inside is produced by the reproducible Docker image workflow, which owns the datum-grid reproducibility that makes golden-file comparison deterministic across machines. Upstream, the schema contracts these tests assert against are the same ones enforced at ingestion time, so a contract test failing in CI is an early warning that a producer has drifted from the agreed data shape long before that drift reaches an operational dashboard. Every merge blocked by a red spatial suite is a wrong geometry that never became an Open Geospatial Consortium (OGC) feature service response an incident commander would have acted on.

Troubleshooting

Symptom: golden-file tests pass locally but fail in CI with tiny coordinate differences. The local and CI environments resolve datum shifts with different PROJ data grids or GEOS versions. Run the suite inside the pinned toolchain container everywhere, print pyproj.proj_version_str and shapely.geos_version_string at the top of the job, and re-baseline goldens only from inside that image so the expected output matches what CI computes.

Symptom: a property-based transform test fails intermittently on different runs. The hypothesis seed is unset, so each run explores different coordinates and occasionally hits a boundary the code mishandles. Pin --hypothesis-seed in CI to make failures reproducible, then read the shrunk counterexample hypothesis prints — it is the minimal (lon, lat) that breaks the invariant — and add it as an explicit regression example so the bug stays caught.

Symptom: geometry-equality assertions flake even when the logic is unchanged. Exact coordinate equality is comparing floating-point values perturbed by reprojection and serialization. Replace == with shapely.equals_exact at a declared tolerance or a symmetric-difference area threshold, and align both sides on a stable key before comparing so row-order differences cannot masquerade as geometry drift.

Symptom: the suite is too slow to run on every push, so developers skip it. The pyramid is inverted — too many slow end-to-end tests and too few fast unit assertions. Move coverage down into geometry, topology, and property-based tests that need no external services, mark the handful of genuine integration tests with a slow marker, and keep the default pytest invocation to the fast tiers so the merge gate stays cheap enough to always run.

Symptom: a green suite still let a wrong buffer distance through. The tests assert that output is valid but never assert the magnitude of the result. Add value-level assertions — buffered area strictly greater than input area, distances within an expected band, feature counts unchanged by a no-op — because for spatial code, structurally valid and operationally correct are different claims that both need a test.

Frequently Asked Questions

Why can’t I just assert exact coordinate equality in spatial tests? Floating-point reprojection, GEOS version differences, and serialization round-trips perturb coordinates in the last few decimal places, so exact equality produces flaky failures that erode trust in the suite. Compare geometries with a tolerance-aware predicate such as equals_exact with a decimal precision, or a symmetric-difference area threshold, so that meaningful changes fail the build while sub-millimetre noise does not.

Do I need to pin GDAL and PROJ versions in CI? Yes. Different PROJ data grids and GEOS releases can shift transformed coordinates by up to a metre and change which geometries are considered valid, so an unpinned runner makes golden-file tests pass locally and fail in CI for reasons unrelated to the code change. Pin the GDAL, PROJ, GEOS, and PROJ-data versions in the test container image so every run resolves datum shifts and topology identically.

How much spatial testing belongs in unit tests versus integration tests? Most coverage should be fast unit-level geometry and topology assertions and property-based transform checks that need no external services; these form the wide base of the test pyramid. Golden-file regression sits in the middle, and a small number of slower end-to-end tests exercise the real ingest-to-publish path. Keeping the base wide keeps the suite fast enough to run on every push, which is what makes it a genuine merge gate.

Related

• Writing Property-Based Tests for Coordinate Transforms
• Setting Up Dockerized GIS Environments
• Version Control for Spatial Workflows

Up: Python Toolchains for Public Safety GIS