Spatial Epi Public Health GIS

Property-Based Tests for Spatial Transforms

This guide, part of Testing Geospatial Analysis Code, shows how to test coordinate and geometry transformations by generating hundreds of inputs and asserting properties, instead of hand-writing a few examples and their expected outputs.

Problem Context & Constraints

Hand-written tests for spatial code share a blind spot: the author chooses inputs they have thought about. The inputs that break spatial code are the ones nobody thought about — a longitude of exactly 180, a polygon whose two vertices coincide, a point that sits precisely on a tract boundary, a ring with a spike of zero width, a coordinate in the southern hemisphere in a pipeline written in the northern one.

Property-based testing inverts the approach. You describe the space of valid inputs and a property every output must satisfy, and the framework generates inputs — biased toward edges and extremes — until it finds one that violates the property. When it does, it shrinks that input to the simplest example that still fails, which is usually a far better bug report than the original.

Three constraints shape how to use it for spatial work. The input strategy must produce valid inputs, or the test spends its time rejecting garbage; generating random coordinate lists produces mostly invalid polygons. The property must be genuinely invariant — something that holds for every valid input, not merely for typical ones. And tolerances must be explicit, because floating-point geometry is approximate and an exact-equality property will fail on legitimate rounding.

Where Hand-Written and Generated Inputs Land Two panels showing the space of valid longitude and latitude. In the hand-written panel, the handful of test inputs cluster in the middle of the continental study area where the author works. In the generated panel, inputs are spread across the whole valid space and are deliberately dense along the edges: the antimeridian, the poles, the equator, the boundaries of the working projection's area of use. The failures in spatial code live at those edges, which is why generation finds them and hand-writing does not. The bugs live at the edges, which is exactly where people do not test Hand-written five points near home Generated dense at the edges and the equator The generator is biased toward extremes on purpose

Prerequisites

  • python 3.11
  • hypothesis 6.103.1, pytest 8.2.2
  • shapely 2.0.4, pyproj 3.6.1, geopandas 0.14.4, numpy 1.26.4
  • The functions under test, importable without touching protected data
  • A declared working coordinate reference system and its area of use

Step-by-Step Solution

Build strategies that generate valid coordinates and valid polygons, then assert invariants for each transformation the pipeline performs.

# Property-based tests for coordinate and geometry transformations.
# Pinned: hypothesis==6.103.1, pytest==8.2.2, shapely==2.0.4, pyproj==3.6.1, numpy==1.26.4
import math

import numpy as np
import pytest
from hypothesis import assume, given, settings, strategies as st
from pyproj import CRS, Transformer
from shapely.geometry import Point, Polygon
from shapely.validation import make_valid

WORK = CRS.from_epsg(5070)
AREA = WORK.area_of_use                      # west, south, east, north in degrees

# --- Strategies ------------------------------------------------------------

lon_in_area = st.floats(min_value=AREA.west, max_value=AREA.east, allow_nan=False)
lat_in_area = st.floats(min_value=AREA.south, max_value=AREA.north, allow_nan=False)


@st.composite
def star_polygon(draw, cx=1_000_000.0, cy=2_000_000.0):
    """Valid simple polygons: vertices at increasing angles, random radii.

    Sorting by angle guarantees a simple (non-self-intersecting) ring, so the
    test spends its budget on valid inputs instead of rejecting garbage.
    """
    n = draw(st.integers(min_value=3, max_value=40))
    angles = sorted(draw(st.lists(st.floats(0, 2 * math.pi, exclude_max=True),
                                  min_size=n, max_size=n, unique=True)))
    radii = draw(st.lists(st.floats(10.0, 50_000.0), min_size=n, max_size=n))
    pts = [(cx + r * math.cos(a), cy + r * math.sin(a)) for a, r in zip(angles, radii)]
    poly = Polygon(pts)
    assume(poly.is_valid and poly.area > 1.0)
    return poly


# --- Invariants --------------------------------------------------------------

@given(lon_in_area, lat_in_area)
@settings(max_examples=800, deadline=None)
def test_round_trip_within_a_centimetre(lon, lat):
    fwd = Transformer.from_crs(4326, WORK, always_xy=True)
    inv = Transformer.from_crs(WORK, 4326, always_xy=True)
    x, y = fwd.transform(lon, lat)
    lon2, lat2 = inv.transform(x, y)
    # About 1e-7 degrees is a centimetre; generous for any epidemiological use.
    assert math.isclose(lon, lon2, abs_tol=1e-7)
    assert math.isclose(lat, lat2, abs_tol=1e-7)


@given(star_polygon(), st.floats(1.0, 5_000.0))
def test_buffer_contains_and_grows(poly, d):
    b = poly.buffer(d)
    assert b.contains(poly.representative_point())
    assert b.area >= poly.area


@given(star_polygon(), star_polygon(cx=1_010_000.0))
def test_intersection_area_bounded(a, b):
    i = a.intersection(b)
    tol = 1e-6 * max(a.area, b.area)
    assert i.area <= a.area + tol and i.area <= b.area + tol


@given(star_polygon())
def test_area_preserved_by_equal_area_projection(poly):
    """Re-projecting through geographic coordinates and back must not change area."""
    to_geo = Transformer.from_crs(WORK, 4326, always_xy=True)
    back = Transformer.from_crs(4326, WORK, always_xy=True)
    xs, ys = poly.exterior.coords.xy
    gx, gy = to_geo.transform(np.asarray(xs), np.asarray(ys))
    rx, ry = back.transform(gx, gy)
    again = Polygon(zip(rx, ry))
    assert math.isclose(again.area, poly.area, rel_tol=1e-6)


@given(star_polygon())
def test_make_valid_is_idempotent(poly):
    once = make_valid(poly)
    twice = make_valid(once)
    assert once.equals(twice)

The polygon strategy is the part worth reusing. Generating vertices at sorted angles around a centre guarantees a simple ring, so every generated input is valid and the test budget is spent checking behaviour rather than discarding malformed shapes. Varying the radius per vertex produces concave, spiky and near-degenerate polygons, which is where overlay and buffer code breaks.

The area-preservation test is specific to epidemiology, where equal-area projections are chosen precisely so that rates are comparable across a map. It asserts that sending a polygon out to geographic coordinates and back does not change its area, which catches a wrong axis order, a wrong datum, or a transform silently falling back to a less accurate method.

Reading a Shrunk Failure A sequence of three polygons from a failing property test. The first is the original generated input, a thirty-one vertex star with irregular radii, on which the buffer containment property failed. The framework removes vertices and simplifies radii while the property still fails, producing an intermediate twelve-vertex shape and finally a four-vertex polygon with one extremely short edge. The minimal example makes the cause obvious: a near-zero-length edge produced a degenerate buffer, which the thirty-one-vertex original concealed. The framework hands you the smallest input that still fails generated: 31 vertices shrinking: 12 vertices minimal: 4 vertices one near-zero edge The 31-vertex original hid the cause; the 4-vertex minimum states it Save the minimal example as a permanent regression case

Validation & Edge Cases

Strategies that reject most inputs. If assume discards most generated values, the framework warns and the test is weak. Build the validity into the generator, as the angular strategy does, rather than filtering after the fact.

Properties that are almost invariant. “The centroid lies inside the polygon” is false for concave shapes and will fail correctly. Use a representative point when inside-ness is what you mean.

Tolerances that are too tight. Floating-point geometry accumulates error. A relative tolerance scaled to the operand size is usually right; an absolute tolerance in metres fails on large polygons and passes wrongly on tiny ones.

Inputs outside the working projection’s area of use. A projected system is only valid over its area of use, and generating coordinates outside it tests behaviour nobody should rely on. Constrain the strategy to the area of use and test out-of-area handling separately, as an explicit rejection.

Invariants for the analysis layer, not only geometry. The same approach applies above the geometry library. A spatial weights builder should produce a symmetric matrix for contiguity on any generated tessellation; a rate calculation should return the same answer when every numerator and denominator is multiplied by the same constant only if it is a ratio, and a different answer if it applies a small-number adjustment; an aggregation to a coarser level should never produce a total different from the sum of its parts. Each of these is a one-line property that checks a piece of epidemiological logic against every input the generator can think of.

Slow properties that time out. Overlay-heavy properties on complex generated polygons can take seconds per example. Disable the per-example deadline for those tests, reduce the example count, and move them to a nightly suite rather than letting a timeout masquerade as a failure on every commit.

Nondeterminism in the code under test. A property that fails intermittently is often exposing an unseeded random component or an unstable sort inside the function. That is a real defect worth fixing, not a flaky test to be retried until it passes.

Found failures that are not regressions yet. When the framework finds a failure, save the shrunk example as an explicit test case. Otherwise the database of found examples is the only record, and it lives on one machine.

Falsifying example: test_buffer_contains_and_grows(
    poly=,
    d=1.0,
)

The tolerance question in that list is concrete: a fixed absolute tolerance is simultaneously too tight for large polygons and too loose for small ones:

Scale the Tolerance to the Operand Accumulated floating-point error in a computed area plotted against polygon area, rising in proportion to size. A fixed absolute tolerance is a horizontal line: above it on the right, large polygons fail tests despite being correct, and far below it on the left, tiny polygons pass tests even when their computed area is wrong by more than itself. A relative tolerance scales with area and sits just above the error curve across the whole range. A fixed tolerance fails big shapes and forgives small ones floating-point error fixed absolute tolerance relative tolerance small: wrongly passes large: wrongly fails polygon area area error

Compliance Notes

Property-based tests are code-quality controls rather than disclosure controls, but they carry one compliance obligation: they must never generate or embed real addresses or case locations. Generated coordinates within the study area are synthetic by construction; a test that loads a sample of real points to “make it realistic” has put a fragment of the case file into the repository and into continuous integration logs.

Record the strategies, the number of examples per test, the tolerances and the pinned library versions in the test configuration, so a reviewer assessing the pipeline’s validation can see exactly what was checked and how strictly. Store every shrunk failing example as a named regression test, which turns a transient discovery into a permanent part of the evidence that the method was validated.