What to Test in a GIS Pipeline (and What You Cannot) #

Problem statement #

Every guide to testing says "write tests". Almost none of them say what to assert on when the output is a polygon.

The usual advice transfers badly. A web developer asserts response.status_code == 200 and moves on. The spatial equivalent — did this clip produce the right shape? — has no single right answer, because:

# these two polygons are "the same" by any reasonable standard
a = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
b = Polygon([(1, 1), (0, 1), (0, 0), (1, 0)])

a == b                    # False — different vertex order
a.equals(b)               # True  — same point set
a.wkt == b.wkt            # False
a.area == b.area          # True, but so is every other unit square

So people write no tests at all, or they write the one test that is easy to write — comparing output files byte for byte — and it breaks on the next GEOS upgrade for reasons nobody can explain.

The problem is not knowing how to test. It is knowing which properties of spatial output are worth asserting on, and which are noise that will make the suite fail for reasons that do not matter.

Quick answer #

Assert on properties that would be wrong if the code were wrong, and that stay the same when nothing meaningful changes:

def test_clip_keeps_crs_and_shrinks_area(parcels, boundary):
    out = clip_to_boundary(parcels, boundary)

    assert out.crs == parcels.crs                     # a silent CRS drop is the classic bug
    assert len(out) <= len(parcels)                   # clip can only remove
    assert out.is_valid.all()                         # clip can produce slivers
    assert out.geometry.area.sum() < parcels.geometry.area.sum()
    assert set(out.columns) == set(parcels.columns)   # attributes survived

Four cheap assertions. Between them they catch a dropped CRS, a broken predicate, an invalid result and a lost attribute table — which is most of what actually goes wrong.

The two kinds of failure #

Loud failures — a missing file, a bad column name — announce themselves the first time you run the script. You do not need a test suite to find them; you need one run.

Silent failures produce output. A spatial join with mismatched CRS returns zero rows, which looks like "no matches in this area". A buffer in degrees returns polygons, just enormous ones. A clip against the wrong boundary returns a perfectly valid subset of the wrong place. Tests exist for the second kind. If a test would only ever catch something an exception already catches, it is not earning its runtime.

Step-by-step solution #

1. Structure — does the output have the right shape? #

The cheapest and most valuable assertions. They run in microseconds and catch schema drift, silent column renames, and functions that quietly return the input unchanged.

def test_output_schema(result):
    assert isinstance(result, gpd.GeoDataFrame)
    assert result.geometry.name == "geometry"
    assert set(result.columns) >= {"parcel_id", "area_m2", "geometry"}
    assert result["parcel_id"].dtype == "int64"
    assert result["parcel_id"].is_unique

is_unique is the underrated one. A spatial join that duplicates rows is the single most common way a GIS result becomes quietly wrong — the map looks fine and the totals are double.

2. Geometry — is it the right kind of shape, and is it valid? #

def test_geometry_properties(result):
    assert result.crs is not None
    assert result.crs.is_projected                    # area/length are meaningless otherwise
    assert (result.geom_type == "Polygon").all()      # not MultiPolygon, not GeometryCollection
    assert result.is_valid.all()
    assert not result.geometry.is_empty.any()

is_empty deserves its own line. An empty geometry is not null, so dropna() misses it; it is valid, so is_valid passes it; and it writes to a file without complaint. It then contributes nothing to every downstream operation.

3. Values — are the numbers in the right range? #

Never assert an exact area. Assert a range, a ratio, or a relationship:

def test_buffer_area_is_plausible(points_27700):
    out = buffer_metres(points_27700, 50)

    expected = math.pi * 50 ** 2
    assert out.geometry.area.iloc[0] == pytest.approx(expected, rel=0.01)

    # the relationship matters more than the number
    wider = buffer_metres(points_27700, 100)
    assert (wider.geometry.area > out.geometry.area).all()

The second assertion survives a change in buffer resolution. The first would too, at rel=0.01 — but a stricter tolerance would break when GEOS changes how many segments it uses per quarter circle, which is exactly the kind of failure that teaches people to ignore the test suite.

4. Invariants — what must be true of any correct run? #

Invariants are the highest-value tests in spatial work because they hold for every input, so one test covers cases you never thought to write fixtures for.

def test_dissolve_preserves_total_area(parcels):
    out = dissolve_by_ward(parcels)
    assert out.geometry.area.sum() == pytest.approx(parcels.geometry.area.sum(), rel=1e-6)

def test_clip_is_idempotent(parcels, boundary):
    once = clip_to_boundary(parcels, boundary)
    twice = clip_to_boundary(once, boundary)
    assert len(once) == len(twice)

def test_reproject_round_trip(parcels):
    there_and_back = parcels.to_crs(4326).to_crs(parcels.crs)
    assert there_and_back.geometry.geom_equals_exact(parcels.geometry, tolerance=0.001).all()

Dissolve must conserve area. Clip must be idempotent. Reprojection must round-trip. None of these needs a carefully constructed expected result — they compare the code against itself.

5. Golden results — the last resort #

A golden test stores a known-good output and compares against it. It catches anything the rungs above miss, and it fails for every irrelevant reason too.

def test_matches_golden(tmp_path):
    out = run_pipeline(FIXTURE_DIR, tmp_path)
    golden = gpd.read_file("tests/golden/wards_summary.gpkg")

    assert_geodataframe_equal(
        out.sort_values("ward_id").reset_index(drop=True),
        golden.sort_values("ward_id").reset_index(drop=True),
        check_less_precise=True,
    )

Use one, on one representative dataset, and expect to regenerate it deliberately when the pipeline changes. A suite of twenty golden tests is a suite nobody maintains.

Code examples #

Example 1: a reusable assertion helper #

Most spatial tests repeat the same four checks. Extract them.

# tests/asserts.py
import geopandas as gpd

def assert_clean_layer(gdf, *, crs=None, geom_type=None, min_rows=1):
    """The assertions every layer in this pipeline must satisfy."""
    assert isinstance(gdf, gpd.GeoDataFrame), f"expected GeoDataFrame, got {type(gdf)}"
    assert len(gdf) >= min_rows, f"only {len(gdf)} rows"
    assert gdf.crs is not None, "layer has no CRS"
    if crs is not None:
        assert gdf.crs == crs, f"expected {crs}, got {gdf.crs}"
    if geom_type is not None:
        bad = gdf[gdf.geom_type != geom_type]
        assert bad.empty, f"{len(bad)} rows are not {geom_type}"
    assert gdf.is_valid.all(), f"{(~gdf.is_valid).sum()} invalid geometries"
    assert not gdf.geometry.is_empty.any(), "layer contains empty geometries"
    assert not gdf.geometry.isna().any(), "layer contains null geometries"

Now each test is one line of intent plus one line of specifics:

def test_clip_output(parcels, boundary):
    out = clip_to_boundary(parcels, boundary)
    assert_clean_layer(out, crs=parcels.crs, geom_type="Polygon")
    assert len(out) <= len(parcels)

Example 2: asserting on the change, not the result #

Where the absolute answer is hard to pin down, the delta usually is not.

def test_cleaning_removes_only_what_it_should(dirty_parcels):
    before = {
        "rows": len(dirty_parcels),
        "invalid": (~dirty_parcels.is_valid).sum(),
        "null_geom": dirty_parcels.geometry.isna().sum(),
    }
    out = clean(dirty_parcels)

    assert (~out.is_valid).sum() == 0
    assert out.geometry.isna().sum() == 0
    # everything removed should be accounted for by the two counts above
    assert len(out) == before["rows"] - before["null_geom"]

That last line is the real test: it asserts that cleaning removed the null geometries and repaired the invalid ones rather than dropping them. A clean() that silently discarded 200 invalid parcels would pass a naive "output is valid" test and fail this one.

Example 3: what a test for a silent CRS bug looks like #

def test_area_calculation_refuses_geographic_crs():
    gdf = gpd.GeoDataFrame(
        geometry=[Polygon([(-3.2, 55.9), (-3.1, 55.9), (-3.1, 56.0), (-3.2, 56.0)])],
        crs="EPSG:4326",
    )
    with pytest.raises(ValueError, match="projected"):
        add_area_column(gdf)

def test_area_calculation_accepts_projected_crs(parcels_27700):
    out = add_area_column(parcels_27700)
    assert out["area_m2"].min() > 0
    assert out["area_m2"].max() < 1e9        # nothing in this dataset is 1,000 km²

The upper bound is a sanity rail, not a precise expectation. It fires the moment somebody computes area in degrees and gets a number seven orders of magnitude too small — or too large, depending on which way the mistake went.

Explanation #

The reason spatial testing feels harder than it is comes down to one thing: geometry has many representations of the same truth. The same square can be wound clockwise or anticlockwise, start at any vertex, carry a redundant collinear point, or be stored as a Polygon or a one-part MultiPolygon. All of those are the same shape. None of them are the same bytes.

Every brittle spatial test is a test that accidentally asserted on the representation instead of the truth. Byte comparison, WKT comparison, exact coordinate comparison and == on geometries all do this. They pass today and fail when GEOS changes its vertex ordering, when GDAL changes its default precision, or when someone runs the pipeline on a machine with a different PROJ version.

Assertions on derived properties — area, length, count, validity, type, CRS, bounds — are stable under all of those changes and unstable under the changes you care about. That is the whole trick.

There is a second, quieter reason to prefer property assertions: they document intent. assert out.crs.is_projected tells the next reader that this function's output is meant to be measurable. A byte comparison tells them nothing at all.

Edge cases or notes #

• Empty results are a legitimate outcome, so assert deliberately. assert len(out) > 0 is right when the fixture guarantees matches, and wrong when the function is supposed to return nothing for a disjoint input. Write both tests.
• geom_equals_exact needs a tolerance and respects vertex order; equals ignores order but is slower. For test assertions, equals is usually what you mean.
• Areas differ between GEOS versions at the tenth decimal place. Any tolerance tighter than rel=1e-9 will eventually fail on somebody else's machine.
• Test the writer, not just the transform. Shapefile truncates column names to 10 characters and cannot store a datetime; a transform that passes every test can still produce a file that loses data on write. See how to fix truncated shapefile column names.
• Do not test GeoPandas. test_that_buffer_returns_polygons tests a dependency, not your code. Test the function that calls it, including the argument you might get wrong.
• Randomised fixtures need a fixed seed. A test that fails one run in fifty is worse than no test.

Internal links #

• How to test a GIS pipeline with pytest — the mechanics: fixtures, structure, running the suite
• Tests pass locally but fail in CI — when the assertions are right but the environment is not
• How to run a Python GIS pipeline in CI with GitHub Actions — where these tests actually run
• How to validate pipeline inputs and outputs automatically — the same assertions, applied at runtime rather than in tests
• Spatial data quality: the six dimensions that matter — what "correct" means before you assert on it
• What makes a geometry valid? — what is_valid is actually checking
• Coordinate precision and floating point in GIS explained — why exact coordinate comparison fails

FAQ #

What is the single most valuable assertion in a GIS test suite? #

assert out.crs == expected_crs. A dropped or wrong CRS is the most common silent failure in Python GIS, it produces plausible output every time, and one line catches it.

Should I compare geometries with ==? #

No. == on two geometry objects compares them element-wise as a Series, and on two individual shapes compares representation rather than shape. Use .equals() for "same shape", or .geom_equals_exact(other, tolerance=...) when vertex order also matters.

How do I test something that depends on a real dataset? #

Cut a tiny extract — five to twenty features — commit it, and test against that. If the logic genuinely needs scale, that is a performance test, not a correctness test, and it belongs behind a marker so it does not run on every commit.

Is it worth testing that a file was written? #

Testing that the file exists is nearly worthless; testing that reading it back gives you the same layer is worth a lot. The round trip is where shapefile column truncation, dropped CRS and coerced dtypes show up.

How many tests does a pipeline actually need? #

One per silent failure mode you can name. That is usually five to fifteen, and it is far more valuable than a hundred tests that assert good input produces some output.

What about testing map output? #

Compare the data that feeds the map, not the image. Image comparison fails on font rendering, matplotlib versions and antialiasing — all noise. If you must, assert on figure structure: number of artists, axis limits, colormap name.

Do property-based tests make sense here? #

Yes, for invariants. Generating random valid polygons with hypothesis and asserting that dissolve conserves area, or that buffer then negative-buffer approximately restores the original, finds edge cases that hand-written fixtures never reach.