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:

from shapely.geometry import Polygon

# 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)])

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

a == b       False
a.equals(b)  True
a.wkt == b.wkt False
a.area == b.area True

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.

Every test on this page is written to a file and run with pytest against real data: Oakham's Land Registry parcels and Rutland's ONS wards. 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:

%%writefile tests/test_quick.py
from pipeline import clip_to_boundary


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

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings tests/test_quick.py

.                                                                        [100%]
1 passed in 0.33s

Five 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. Here is the first of those, on the real wards:

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    joined = gpd.sjoin(parcels, wards.to_crs(4326), predicate="within")
print(f"{len(joined)} rows from a join of {len(parcels):,} parcels to wards in the wrong CRS; "
      f"{len(caught)} warning, no exception")

0 rows from a join of 5,006 parcels to wards in the wrong CRS; 1 warning, no exception

Structure, geometry, values, invariants, golden results

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.

%%writefile tests/test_structure.py
import geopandas as gpd
import pytest

from pipeline import area_by_ward


@pytest.fixture
def result(parcels, wards):
    return area_by_ward(parcels, wards)


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

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings -rf tests/test_structure.py 2>&1 | grep -E "^FAILED|passed|failed" | cut -c1-110

FAILED tests/test_structure.py::test_output_schema - assert False
1 failed in 0.35s

is_unique is the underrated one, and here it fails for real: area_by_ward() joins parcels to wards with intersects, and every parcel that straddles a ward boundary comes back once per ward. A spatial join that duplicates rows is the most common way a GIS result becomes quietly wrong — the map looks fine and the totals are inflated. The dtype line matters too: the INSPIRE IDs arrive as int32, and a test written from memory as int64 would fail for a reason that has nothing to do with the code.

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

%%writefile tests/test_geometry.py
from pipeline import clip_to_boundary


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


def test_clip_makes_multipolygons(parcels, boundary):
    result = clip_to_boundary(parcels, boundary)
    print("\n", result.geom_type.value_counts().to_dict())

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings -s tests/test_geometry.py

.
 {'Polygon': 1631}
.
2 passed in 0.32s

Every row came back a Polygon, and that is itself worth a second look. A plain gpd.clip of these parcels to Oakham North East returns 1,609 rows, three of them MultiPolygon and five GeometryCollection; with keep_geom_type=True it returns 1,631 all-Polygon rows, because it splits those multi-part results into separate rows that share an INSPIREID. Assert the set of types the next step can handle, and the uniqueness of the key. 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:

%%writefile tests/test_values.py
import math

import pytest

from pipeline import buffer_metres


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()


def test_exact_area_is_brittle(points_27700):
    assert buffer_metres(points_27700, 50).geometry.area.iloc[0] == pytest.approx(math.pi * 50 ** 2, rel=1e-4)

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings tests/test_values.py 2>&1 | grep -E "^E +assert|passed|failed" | head -3

E       assert np.float64(7841.371226366008) == 7853.981633974483 ± 0.785398
E         comparison failed
1 failed, 1 passed in 0.18s

The relationship survives a change in buffer resolution, and so does the area at rel=0.01. At rel=1e-4 the same area fails — 7,841.4 m² against πr² = 7,854.0 m² — because a default buffer is a 64-sided polygon, 0.16% smaller than the circle. A tolerance tighter than the approximation tests GEOS's segment count rather than your code, 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.

%%writefile tests/test_invariants.py
import pytest

from pipeline import clip_to_boundary, dissolve_all


def test_dissolve_preserves_total_area_of_a_coverage(wards):
    out = dissolve_all(wards)
    assert out.geometry.area.sum() == pytest.approx(wards.geometry.area.sum(), rel=1e-6)


def test_dissolve_preserves_total_area_of_parcels(parcels):
    out = dissolve_all(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()

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings -rf tests/test_invariants.py 2>&1 | grep -E "^FAILED|passed|failed" | cut -c1-110

E         comparison failed
FAILED tests/test_invariants.py::test_dissolve_preserves_total_area_of_parcels
FAILED tests/test_invariants.py::test_clip_is_idempotent - assert 1631 == 1628
2 failed, 2 passed in 0.89s

Two of the four fail, and both failures are findings. Dissolve must conserve area — of a coverage: Rutland's wards pass, and Oakham's parcels do not, because registered parcels overlap, so the dissolved area is smaller than the sum of the parts. And the clip is not idempotent: the first pass returns 1,631 rows, the second 1,628, because keep_geom_type=True had split multi-part results into rows with duplicate IDs and the second pass handles those rows differently again. An invariant states an assumption about the data and the code, and a real fixture is what tests it. The reprojection round trip holds to a millimetre here — like the others, it compares the code against itself, with no carefully constructed expected result.

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.

%%writefile tests/test_golden.py
from pathlib import Path

import geopandas as gpd
import pytest
from geopandas.testing import assert_geodataframe_equal

from pipeline import clip_to_boundary

GOLDEN = Path("tests/golden/oakham_ne.gpkg")


def test_matches_golden(parcels, boundary):
    out = clip_to_boundary(parcels, boundary).sort_values("INSPIREID").reset_index(drop=True)
    if not GOLDEN.exists():
        GOLDEN.parent.mkdir(parents=True, exist_ok=True)
        out.to_file(GOLDEN)
        pytest.skip("golden file written — review and commit it")
    golden = gpd.read_file(GOLDEN).sort_values("INSPIREID").reset_index(drop=True)
    assert_geodataframe_equal(out, golden, check_less_precise=True, check_geom_type=False)

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings -rs tests/test_golden.py | tail -1
!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings tests/test_golden.py | tail -1

1 skipped in 0.26s
1 passed in 0.32s

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.

Assertions on Oakham parcels and a silent CRS bug

The data used on this page

HM Land Registry's INSPIRE Index Polygons for Rutland (© Crown copyright and database rights; OGL), cut to Oakham, and Rutland's 15 wards from the ONS December 2024 boundaries (OGL). The setup blocks write the module under test and the pytest fixtures; every test file on this page imports from them, so they run before anything else.

import sys
import warnings
from pathlib import Path

import geopandas as gpd
import requests

for folder in ("data", "tests"):
    Path(folder).mkdir(exist_ok=True)
if not Path("data/Rutland.zip").exists():
    Path("data/Rutland.zip").write_bytes(requests.get(
        "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip", timeout=300).content)
if not Path("data/wards.geojson").exists():
    Path("data/wards.geojson").write_bytes(requests.get(
        "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
        "Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
        params={"where": "LAD24CD='E06000017'", "outFields": "WD24NM", "outSR": 27700, "f": "geojson"},
        timeout=120).content)
wards = gpd.read_file("data/wards.geojson")
oakham = wards[wards["WD24NM"].str.startswith("Oakham")].union_all()
parcels = gpd.read_file("/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml",
                        columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO")
parcels = parcels[parcels.intersects(oakham)].reset_index(drop=True)
parcels.to_file("data/oakham_parcels.gpkg")
print(f"{len(parcels):,} Oakham parcels, {len(wards)} wards")

5,006 Oakham parcels, 15 wards

%%writefile pipeline.py
"""The functions this page tests."""
import geopandas as gpd


def clip_to_boundary(gdf, boundary):
    return gpd.clip(gdf, boundary.to_crs(gdf.crs), keep_geom_type=True)


def area_by_ward(parcels, wards):
    joined = gpd.sjoin(parcels, wards[["WD24NM", "geometry"]], predicate="intersects").drop(columns="index_right")
    return joined.assign(area_m2=joined.area)


def buffer_metres(gdf, distance):
    if not gdf.crs.is_projected:
        raise ValueError("buffer_metres needs a projected CRS")
    return gdf.assign(geometry=gdf.buffer(distance))


def dissolve_all(gdf):
    return gdf[["geometry"]].dissolve()


def clean(gdf):
    out = gdf[gdf.geometry.notna()].copy()
    out["geometry"] = out.geometry.make_valid()
    return out


def add_area_column(gdf):
    if not gdf.crs.is_projected:
        raise ValueError("add_area_column needs a projected CRS")
    return gdf.assign(area_m2=gdf.area)

%%writefile tests/conftest.py
import geopandas as gpd
import pytest
import shapely
from shapely.geometry import Point


@pytest.fixture(scope="session")
def parcels():
    return gpd.read_file("data/oakham_parcels.gpkg")


@pytest.fixture(scope="session")
def wards():
    return gpd.read_file("data/wards.geojson")


@pytest.fixture(scope="session")
def boundary(wards):
    return wards[wards["WD24NM"] == "Oakham North East"]


@pytest.fixture
def points_27700():
    return gpd.GeoDataFrame(geometry=[Point(486000, 308800)], crs="EPSG:27700")


@pytest.fixture
def dirty_parcels(parcels):
    """Ten real parcels: two with the geometry removed, two made into bow-ties."""
    dirty = parcels.head(10).copy()
    dirty.loc[[0, 1], "geometry"] = None
    for i in (2, 3):
        ring = list(dirty.geometry[i].exterior.coords)
        ring[1], ring[2] = ring[2], ring[1]
        dirty.loc[i, "geometry"] = shapely.Polygon(ring)
    return dirty

Example 1: a reusable assertion helper

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

%%writefile tests/asserts.py
import geopandas as gpd


def assert_clean_layer(gdf, *, crs=None, geom_types=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_types is not None:
        bad = gdf[~gdf.geom_type.isin(geom_types)]
        assert bad.empty, f"{len(bad)} rows are not {geom_types}"
    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:

%%writefile tests/test_helper.py
from asserts import assert_clean_layer
from pipeline import clip_to_boundary


def test_clip_output(parcels, boundary):
    out = clip_to_boundary(parcels, boundary)
    assert_clean_layer(out, crs=parcels.crs, geom_types={"Polygon", "MultiPolygon"})
    assert len(out) <= len(parcels)


def test_helper_reports_the_problem(parcels, boundary):
    out = clip_to_boundary(parcels, boundary)
    assert_clean_layer(out, geom_types={"MultiPolygon"})

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings -rf tests/test_helper.py 2>&1 | grep -E "^E +AssertionError|passed|failed" | sed 's/^E *//'

AssertionError: 1631 rows are not {'MultiPolygon'}
1 failed, 1 passed in 0.27s

The second test is written to fail, to show what the helper says when a layer does not meet it: one line naming the check and the count, instead of a bare assert False.

Example 2: asserting on the change, not the result

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

%%writefile tests/test_delta.py
from pipeline import clean


def test_cleaning_removes_only_what_it_should(dirty_parcels):
    before = {"rows": len(dirty_parcels),
              "invalid": int((~dirty_parcels.is_valid).sum()),
              "null_geom": int(dirty_parcels.geometry.isna().sum())}
    out = clean(dirty_parcels)
    print("\n", before, "->", len(out), "rows")
    assert (~out.is_valid).sum() == 0
    assert out.geometry.isna().sum() == 0
    # everything removed should be accounted for by the null count
    assert len(out) == before["rows"] - before["null_geom"]

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings -s tests/test_delta.py | tail -3

 {'rows': 10, 'invalid': 2, 'null_geom': 2} -> 8 rows
.
1 passed in 0.13s

The fixture is ten real Oakham parcels with deliberate damage — two geometries removed and two turned into bow-ties by swapping a pair of vertices, because the INSPIRE polygons themselves are all valid. That last assertion is the real test: it checks that cleaning removed the null geometries and repaired the invalid ones rather than dropping them. A clean() that silently discarded 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

%%writefile tests/test_crs_bug.py
import geopandas as gpd
import pytest
from shapely.geometry import Polygon

from pipeline import add_area_column


def test_area_calculation_refuses_geographic_crs():
    gdf = gpd.GeoDataFrame(geometry=[Polygon([(-0.74, 52.66), (-0.72, 52.66), (-0.72, 52.68), (-0.74, 52.68)])],
                           crs="EPSG:4326")
    with pytest.raises(ValueError, match="projected"):
        add_area_column(gdf)


def test_area_calculation_accepts_projected_crs(parcels):
    out = add_area_column(parcels)
    assert out["area_m2"].min() > 0
    assert out["area_m2"].max() < 1e9        # nothing in Oakham is 1,000 km²
    print("\n largest parcel:", round(out["area_m2"].max()), "m²")

!{sys.executable} -m pytest -q -p no:cacheprovider --color=no --disable-warnings -s tests/test_crs_bug.py | tail -2

.
2 passed in 0.15s

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.

Geometry has many representations of one truth

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.

Empty results, exact equality and platform drift

• 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.
• Floating-point results can differ between library versions and platforms. Keep tolerances at the level your question needs, not at the level today's run happens to reach.
• 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.