Property-based testing for geometry normalization

Your normalisation function has fourteen example-based tests and they all pass. It also has a bug: on a polygon whose lowest-leftmost coordinate appears twice, the start-vertex selection is unstable, so the same geometry canonicalises two different ways depending on ring order — which means the change detector reports a rezone every time the county republishes. Nobody wrote that example because nobody thought of it. A generator would have found it in about forty tries. This guide writes generators that produce the shapes municipal data actually contains, chooses invariants that genuinely hold, and reads the counterexamples they produce. It is the practical form of layer two in testing spatial data pipelines.

Diagnosis: which properties are actually invariant?

The temptation is to assert everything, and over-strong properties fail on legitimate input, get weakened until they assert nothing, and are then deleted. Sort candidate properties into three groups first.

Six candidate properties, sorted by whether they actually hold A six-row by three-column matrix of candidate properties for a geometry normalisation function. Rows are idempotence, area preservation, topological equality, digest equality under ring reversal, vertex count preservation, and bounding box preservation. Columns are whether the property is invariant, invariant only within a tolerance, or not invariant at all, and what the correct form of the assertion is. Vertex count preservation is not invariant because rounding legitimately merges coincident vertices, and the correct assertion is that the count never rises. Six candidate properties, sorted by whether they actually hold invariant? why correct assertion canonical(canonical(g)) == canonical(g) yes, always a normal form is a fixed point assert equality area is preserved within a tolerance grid rounding moves vertices tolerance = perimeter × grid topologically equal to the input within precision rounding is lossy equals within precision digest survives ring reversal yes, always direction carries no meaning assert digest equality vertex count is preserved no rounding merges coincident vertices count never RISES bounding box is preserved no moves by up to half a grid cell within half a cell invariant within a tolerance not invariant
The last two rows look invariant and are not. Rounding to a grid legitimately merges coincident vertices and can move the bounding box by half a cell, so asserting preservation produces failures on correct input — which get "fixed" by weakening the property until it asserts nothing.

Genuinely invariant — must hold for every valid input:

  • canonical(canonical(g)) == canonical(g) — idempotence.
  • canonical(g).area == g.area — normalisation reorders and rounds; it must not resize.
  • canonical(g).equals(g) topologically, within the declared precision.
  • digest(canonical(g)) == digest(canonical(reverse_rings(g))) — ring direction carries no meaning.

Invariant only within a stated tolerance — true given a declared precision:

  • vertex count is preserved, unless rounding collapses two vertices onto one grid cell.
  • round-tripping through WKB is exact, given the same precision model.

Not invariant, though it looks like it — the trap:

  • “the vertex count never changes” is false, because rounding legitimately merges coincident vertices.
  • “the bounding box is unchanged” is false at the precision boundary, by up to half a grid cell.

Getting a property into the wrong group costs an afternoon of arguing with a correct counterexample.

Step-by-step implementation

1. Write generators that produce municipal pathologies

A generator of random polygons mostly produces shapes county data never contains. A generator built from the known pathologies produces failures you will actually see in production.

Bugs found by a generic generator against a municipal one Stacked bar chart comparing two generators over 500 examples each, by the class of bug found. A generic random-polygon generator finds 1 real bug and spends most of its examples on shapes county data never contains, producing 7 failures caused by the generator itself. A generator built from known municipal pathologies finds 6 real bugs, including the duplicate-leading-vertex case, and produces no generator-caused failures. Bugs found by a generic generator against a municipal one 0 100 200 300 400 500 430 generic random polygons 60 municipal pathology generator Of 500 generated examples failures caused by the generator examples on shapes that never occur real bugs found "Real bugs" are ones that reproduced on data from an actual county export.
The generic generator spends its budget on shapes that do not occur and then produces failures that are the generator's fault, which trains everyone to ignore it. Building the generator from known pathologies — duplicate vertices, clockwise rings, slivers, holes — finds six times as many real bugs in the same 500 examples.
from hypothesis import strategies as st
from shapely.geometry import Polygon

# Realistic municipal coordinates: a projected CRS in metres, plausible easting and
# northing. Random floats over the whole double range test the float library, not you.
easting = st.floats(min_value=400_000, max_value=700_000, allow_nan=False,
                    allow_infinity=False, width=64)
northing = st.floats(min_value=3_000_000, max_value=5_000_000, allow_nan=False,
                     allow_infinity=False, width=64)


@st.composite
def parcel_quads(draw):
    """A convex quadrilateral with a skew, which is what most platted lots are."""
    x0, y0 = draw(easting), draw(northing)
    w = draw(st.floats(min_value=0.5, max_value=400.0))
    h = draw(st.floats(min_value=0.5, max_value=400.0))
    skew = draw(st.floats(min_value=-60.0, max_value=60.0))
    clockwise = draw(st.booleans())
    ring = [(x0, y0), (x0 + w, y0 + skew), (x0 + w, y0 + h + skew), (x0, y0 + h)]
    return Polygon(ring[::-1] if clockwise else ring)


@st.composite
def parcels_with_duplicate_vertices(draw):
    """County exports contain repeated consecutive vertices in quantity. This is the
    generator that finds unstable start-vertex selection."""
    base = draw(parcel_quads())
    coords = list(base.exterior.coords)[:-1]
    idx = draw(st.integers(min_value=0, max_value=len(coords) - 1))
    coords.insert(idx, coords[idx])                  # exact duplicate
    return Polygon(coords + [coords[0]])


@st.composite
def parcels_with_holes(draw):
    outer = draw(parcel_quads())
    minx, miny, maxx, maxy = outer.bounds
    if (maxx - minx) < 4 or (maxy - miny) < 4:
        return outer                                  # too small for a hole
    pad = draw(st.floats(min_value=1.0, max_value=min(maxx - minx, maxy - miny) / 3))
    hole = Polygon([(minx + pad, miny + pad), (maxx - pad, miny + pad),
                    (maxx - pad, maxy - pad), (minx + pad, maxy - pad)])
    return Polygon(outer.exterior.coords, [hole.exterior.coords])


municipal_parcels = st.one_of(parcel_quads(), parcels_with_duplicate_vertices(),
                              parcels_with_holes())

2. Assert the invariants, with the tolerance stated in the assertion 

import math

from hypothesis import assume, given, settings

PRECISION = 3                     # millimetre grid
GRID = 10 ** -PRECISION
AREA_REL_TOL = 1e-9


@given(municipal_parcels)
@settings(max_examples=500, deadline=None)
def test_canonical_is_idempotent(poly):
    assume(poly.is_valid)                     # the contract is for valid input
    once = canonical(poly)
    assert canonical(once).equals_exact(once, 0), "canonical form is not stable"


@given(municipal_parcels)
@settings(max_examples=500, deadline=None)
def test_canonical_preserves_area(poly):
    assume(poly.is_valid)
    # Rounding to the grid can move each vertex by up to half a cell, so the area
    # tolerance is derived from the perimeter rather than picked.
    slack = poly.length * GRID
    assert math.isclose(canonical(poly).area, poly.area,
                        rel_tol=AREA_REL_TOL, abs_tol=slack + 1e-9)


@given(municipal_parcels)
@settings(max_examples=500, deadline=None)
def test_ring_direction_carries_no_meaning(poly):
    assume(poly.is_valid)
    reversed_rings = Polygon(list(poly.exterior.coords)[::-1],
                             [list(h.coords)[::-1] for h in poly.interiors])
    assert digest(canonical(poly)) == digest(canonical(reversed_rings))


@given(municipal_parcels)
@settings(max_examples=300, deadline=None)
def test_vertex_count_only_ever_falls(poly):
    """The CORRECT version of "vertex count is preserved": rounding may merge
    coincident vertices, so the count can fall and must never rise."""
    assume(poly.is_valid)
    before = len(poly.exterior.coords) - 1
    after = len(canonical(poly).exterior.coords) - 1
    assert after <= before

The derived slack in the area test is the pattern worth copying. A hard-coded abs_tol=0.01 either fails on a large perimeter or hides a real error on a small one; deriving it from the perimeter and the grid makes the tolerance mean something.

3. Read the counterexample rather than weakening the property

Hypothesis shrinks a failure to something minimal, and the shrunk case usually names the bug directly:

Falsifying example: test_canonical_is_idempotent(
    poly=<POLYGON ((400000 3000000, 400000 3000000, 400000.5 3000000, ...))>
)

Two identical leading vertices. The start-vertex selection uses min() over the coordinate list, min() returns the first of equal minima, and the index of that minimum differs between the original and the once-canonicalised ring. The fix is a total order that breaks ties deterministically:

def _start_index(coords):
    """A total order over coordinates. Ties on (x, y) are broken by index so the
    choice is deterministic even when a coordinate appears twice — which is what
    the duplicate-vertex generator found."""
    return min(range(len(coords)), key=lambda i: (coords[i][0], coords[i][1], i))

4. Keep the failing case as a permanent example 

@pytest.mark.parametrize("poly", [
    # From a Hypothesis counterexample, 2026-08-11: duplicate leading vertex made
    # start-vertex selection unstable, so the change detector saw a phantom rezone
    # on every republication.
    Polygon([(400000, 3000000), (400000, 3000000), (400000.5, 3000000),
             (400000.5, 3000000.5), (400000, 3000000.5)]),
])
def test_regression_duplicate_leading_vertex(poly):
    once = canonical(poly)
    assert canonical(once).equals_exact(once, 0)

Property tests find bugs; example tests pin them. Promoting each counterexample into a named regression test with a comment saying what it broke is what keeps the fix from being undone by a later refactor — and it runs in microseconds rather than five hundred generated cases.

Verification & testing

The meta-question for a property suite is whether it would actually catch a regression, and that is testable by mutation:

def test_the_suite_catches_an_unstable_start_vertex(monkeypatch):
    """Deliberately break canonical() and assert the property test notices. A
    property suite that survives an injected bug is decoration."""
    def unstable_start(coords):
        return min(range(len(coords)), key=lambda i: (coords[i][0], coords[i][1]))

    monkeypatch.setattr("normalise._start_index", unstable_start)
    with pytest.raises(AssertionError):
        test_canonical_is_idempotent()


def test_generators_produce_valid_polygons():
    """A generator emitting invalid geometry makes every property fail for the
    generator's reasons rather than the code's."""
    for _ in range(200):
        poly = municipal_parcels.example()
        assert poly.is_valid or poly.is_empty

Failure recovery

A property fails on legitimate input. The property is wrong, not the code — that is the common case and it is a useful finding. Move it from “invariant” to “invariant within a tolerance,” derive the tolerance from something physical, and record why in a comment. Do not delete it.

Reading a counterexample instead of weakening the property A four-row by three-column matrix of counterexample situations. Rows are a shrunk example with a duplicate leading vertex, a property failing on legitimately valid input, a counterexample nobody can interpret, and a suite that has become too slow. Columns are what it means, the wrong response, and the right response. The wrong response to a property failing on valid input is to delete the property; the right one is to move it into the tolerance category with a derived bound. Reading a counterexample instead of weakening the property what it means the tempting wrong move the right move Duplicate leading vertex unstable tie-break add an example test only total order, then promote the example Fails on legitimately valid input the property is too strong delete the property derive a tolerance, keep it Counterexample is uninterpretable several pathologies at once ignore it narrow the generator, re-run Suite has become too slow expensive GEOS calls cut examples everywhere split fast and nightly suites right the situation wrong
The middle column is the failure mode of property testing as a practice: every one of these has an easy wrong response that makes the suite quieter. Deleting a property that failed on valid input removes the only test that would have caught the next real bug in that area.

A counterexample nobody can interpret. Shrink further by narrowing the generator: constrain to one pathology at a time and re-run. A failure that only appears when a duplicate vertex and a hole and a clockwise ring coincide is three bugs wearing a coat.

The suite becomes too slow. Cut example counts on the expensive properties and raise them on the cheap ones; area and idempotence are fast, and anything calling into GEOS overlay is not. Move heavy generated tests into a nightly suite and keep the promoted regression examples in the fast one.

Frequently asked questions

Why not just write more example tests?

Because the bugs that survive review are the cases nobody thought of, and examples are limited to what you thought of. The duplicate-leading-vertex bug is a good illustration: it is obvious in hindsight, appears in real county exports, and nobody writes it as an example. A generator finds it in tens of tries.

Should generators produce invalid geometry?

Not in the same strategy as valid input. A generator emitting bow-ties makes every property fail for the generator’s reasons, which trains you to ignore failures. Test invalid input explicitly, as examples, asserting that it is rejected — that is a different contract and deserves its own tests.

How do you pick a tolerance for a property?

Derive it from something physical rather than tuning it until the suite is green. For area under grid rounding, perimeter times grid size bounds the error; for a reprojection round trip, the transformation’s own claimed accuracy does. A tolerance with a derivation in a comment survives review; a tuned constant creeps upward until the property asserts nothing.

What do you do with a counterexample after fixing the bug?

Promote it into a named regression test with a comment recording what it broke and when it was found. The property test keeps hunting for new cases; the example pins this one cheaply and permanently, so a later refactor that reintroduces the bug fails in microseconds rather than after five hundred generated examples.