How to Test a GIS Pipeline with pytest

Problem statement

The pipeline runs every night: it takes Rutland's school sites from OpenStreetMap, buffers each by 50 metres, and lists the Land Registry parcels inside those zones for a planning check. Then someone "simplifies" one line:

- zones = buffer_metres(schools.to_crs(27700), 50)
+ zones = schools.buffer(50)            # schools are in EPSG:4326: now 50 degrees

This is what the two versions produce on the real data:

import warnings

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    wrong_zones = schools.buffer(50)
right_zones = schools.to_crs(27700).buffer(50)

right = parcels.intersects(right_zones.union_all()).sum()
wrong = parcels.to_crs(4326).intersects(wrong_zones.union_all()).sum()
print(f"parcels near a school: {right:,} with metres, {wrong:,} with degrees (of {len(parcels):,})")
print("one 'zone' spans latitudes", [round(v, 1) for v in wrong_zones.iloc[0].bounds[1::2]])
print("warnings raised:", [str(w.message)[:60] + "…" for w in caught])

parcels near a school: 1,398 with metres, 21,253 with degrees (of 21,253)
one 'zone' spans latitudes [2.7, 102.6]
warnings raised: ["Geometry is in a geographic CRS. Results from 'buffer' are l…"]

No exception. One UserWarning, in a log nobody reads at night. Every downstream step runs happily on "zones" 100 degrees wide that reach latitude 102, and the report says every parcel in the county is near a school. GIS code is unusually good at failing silently — a CRS mistake produces valid-looking geometry in the wrong place, a bad join produces fewer rows rather than an error, and a broken clip produces an empty layer that writes perfectly.

Tests are the only thing that notices. But testing spatial code has its own problems:

• fixtures need geometry, and real datasets are too big to commit
• floating-point coordinates rarely compare exactly
• the interesting properties are geometric (area, topology, CRS), not textual
• I/O and databases are slow, so a naive test suite takes twenty minutes
• the failure you most want to catch — plausible but wrong output — has no exception to assert on

Quick answer

Test the transformation functions, not the script, and assert on geometric properties:

1. structure the pipeline as small pure functions that take and return GeoDataFrames
2. build tiny in-memory fixtures with shapely — a few features, not a county
3. assert on CRS, feature count, geometry type, validity and area, not on exact coordinates
4. use pytest.approx or a tolerance for coordinate comparisons
5. add one end-to-end test on a small real fixture file that exercises real I/O

%%writefile test_quick.py
import geopandas as gpd
import pytest
from shapely.geometry import Point


def buffer_metres(gdf, distance_m):
    if gdf.crs is None:
        raise ValueError("input has no CRS")
    if gdf.crs.is_geographic:
        raise ValueError(f"buffer_metres needs a projected CRS, got {gdf.crs.name}")
    return gdf.assign(geometry=gdf.geometry.buffer(distance_m))


@pytest.fixture
def points_27700():
    return gpd.GeoDataFrame({"id": [1, 2]}, geometry=[Point(486000, 308800), Point(487500, 309400)],
                            crs="EPSG:27700")


def test_buffer_metres_uses_metres(points_27700):
    out = buffer_metres(points_27700, 50)
    assert out.crs == points_27700.crs
    assert (out.geom_type == "Polygon").all()
    assert out.geometry.area.iloc[0] == pytest.approx(3.14159 * 50**2, rel=0.01)


def test_buffer_metres_rejects_geographic_crs():
    gdf = gpd.GeoDataFrame(geometry=[Point(-0.73, 52.67)], crs="EPSG:4326")
    with pytest.raises(ValueError, match="projected CRS"):
        buffer_metres(gdf, 50)

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

..                                                                       [100%]
2 passed in 0.20s

The second test is the one that would have caught the incident: it asserts that the function refuses to buffer in degrees at all, turning a warning nobody reads into an error nobody can miss. The first passes with a 1% tolerance because a buffer is a polygon with 32 segments, not a circle.

What to test at each level

Make it testable, assert properties, keep slow tests apart

Make the pipeline testable first

A script that reads, transforms and writes in one function can only be tested end to end. Split it into pure transformations and a thin shell:

%%writefile src/transforms.py
"""Pure: GeoDataFrame in, GeoDataFrame out, no I/O."""
import geopandas as gpd


def buffer_metres(gdf: gpd.GeoDataFrame, distance_m: float) -> gpd.GeoDataFrame:
    if gdf.crs is None:
        raise ValueError("input has no CRS")
    if gdf.crs.is_geographic:
        raise ValueError(f"buffer_metres needs a projected CRS, got {gdf.crs.name}")
    return gdf.assign(geometry=gdf.geometry.buffer(distance_m))


def drop_invalid(gdf: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, int]:
    bad = gdf.geometry.isna() | gdf.geometry.is_empty | ~gdf.geometry.is_valid
    return gdf.loc[~bad].copy(), int(bad.sum())


def select_near(features: gpd.GeoDataFrame, zones: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Features that intersect any zone, each once, with their own columns only."""
    if features.crs != zones.crs:
        raise ValueError(f"CRS mismatch: {features.crs} vs {zones.crs}")
    hit = features.sindex.query(zones.geometry, predicate="intersects")[1]
    return features.iloc[sorted(set(hit))].copy()

%%writefile src/pipeline.py
"""The impure shell: reads, calls the transforms, writes."""
import geopandas as gpd

from src.transforms import buffer_metres, drop_invalid, select_near


def run(config):
    schools = gpd.read_file(config["schools"])
    parcels = gpd.read_file(config["parcels"])
    schools, dropped = drop_invalid(schools)
    zones = buffer_metres(schools.to_crs(config["crs"]), config["buffer_m"])
    near = select_near(parcels.to_crs(config["crs"]), zones)
    near.to_file(config["output"], driver="GPKG")
    return {"schools": len(schools), "dropped": dropped, "parcels": len(near)}

Pure functions are where the logic lives and where the tests belong. The shell gets one test, at the end.

Build small fixtures in code

Committing the county's 21,253 parcels to test a selection is how test suites become unusable. Construct exactly what the test needs:

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

BNG = "EPSG:27700"


@pytest.fixture
def square():
    """A 100 m × 100 m square in EPSG:27700."""
    return Polygon([(0, 0), (100, 0), (100, 100), (0, 100)])


@pytest.fixture
def parcels(square):
    return gpd.GeoDataFrame(
        {"id": [1, 2, 3], "class": ["residential", "commercial", "residential"]},
        geometry=[square,
                  Polygon([(200, 0), (300, 0), (300, 100), (200, 100)]),
                  Polygon([(0, 200), (50, 200), (50, 250), (0, 250)])],
        crs=BNG)


@pytest.fixture
def bowtie():
    """A self-intersecting polygon — invalid by construction."""
    return Polygon([(0, 0), (10, 10), (10, 0), (0, 10)])


@pytest.fixture
def grid_polygons():
    """A 3×3 grid of 100 m squares — predictable areas and adjacency."""
    cells = [(f"cell_{i}{j}", Polygon([(i * 100, j * 100), (i * 100 + 100, j * 100),
                                       (i * 100 + 100, j * 100 + 100), (i * 100, j * 100 + 100)]))
             for i in range(3) for j in range(3)]
    return gpd.GeoDataFrame({"id": [c[0] for c in cells]}, geometry=[c[1] for c in cells], crs=BNG)


@pytest.fixture
def random_points():
    rng = np.random.default_rng(42)               # seeded: reproducible failures
    xs, ys = rng.uniform(0, 300, 50), rng.uniform(0, 300, 50)
    return gpd.GeoDataFrame({"value": rng.normal(100, 15, 50)},
                            geometry=[Point(x, y) for x, y in zip(xs, ys)], crs=BNG)


@pytest.fixture
def messy_frame(grid_polygons, bowtie):
    """Everything a cleaning function should cope with: nine good cells, a bowtie, a None, an empty."""
    rows = list(grid_polygons.geometry) + [bowtie, None, Polygon()]
    return gpd.GeoDataFrame({"id": range(len(rows))}, geometry=rows, crs=BNG)


@pytest.fixture
def oakham_config(tmp_path):
    """The small real fixture: Oakham's OSM school sites and Land Registry parcels."""
    return {"schools": "tests/data/oakham_schools.gpkg", "parcels": "tests/data/oakham_parcels.gpkg",
            "output": str(tmp_path / "near_schools.gpkg"), "crs": BNG, "buffer_m": 50}

A square of side 100 has an area of exactly 10,000 — which makes assertions readable and failures obvious. A seeded random generator matters too: an unseeded fixture produces a test that fails once a fortnight and cannot be reproduced. The one real fixture, oakham_config, points at two small files cut from the county data in the setup block below — the real data belongs in the end-to-end test, not in every unit test.

Assert on properties, not on coordinates

%%writefile tests/test_transforms.py
import geopandas as gpd
import pytest
from shapely import equals_exact
from shapely.geometry import Point, Polygon

from src.transforms import buffer_metres, drop_invalid, select_near


def test_select_near_keeps_only_overlapping_parcels(parcels, square):
    zones = gpd.GeoDataFrame(geometry=[square.buffer(10)], crs=parcels.crs)
    near = select_near(parcels, zones)
    assert len(near) == 1                                               # count
    assert near.crs == parcels.crs                                      # CRS preserved
    assert near.geometry.is_valid.all()                                 # validity
    assert (near.geom_type == "Polygon").all()                          # type
    assert near.geometry.area.sum() <= parcels.geometry.area.sum()      # invariant
    assert list(near.columns) == list(parcels.columns)                  # schema


def test_select_near_counts_each_parcel_once(grid_polygons):
    zones = gpd.GeoDataFrame(geometry=[Point(100, 100).buffer(30), Point(110, 110).buffer(30)],
                             crs=grid_polygons.crs)
    assert len(select_near(grid_polygons, zones)) == 4


def test_select_near_refuses_mixed_crs(parcels, square):
    zones = gpd.GeoDataFrame(geometry=[square], crs=parcels.crs).to_crs(4326)
    with pytest.raises(ValueError, match="CRS mismatch"):
        select_near(parcels, zones)


def test_centroid_position(parcels):
    centroid = parcels.geometry.iloc[0].centroid
    assert centroid.x == pytest.approx(50.0, abs=1e-6)
    assert centroid.y == pytest.approx(50.0, abs=1e-6)


def test_reproject_round_trip_in_rutland(parcels):
    here = parcels.translate(486000, 308800)                  # the fixture moved to Oakham
    back = here.to_crs(4326).to_crs(here.crs)
    moved = [a.hausdorff_distance(b) for a, b in zip(here, back)]
    assert not equals_exact(here.iloc[0], back.iloc[0], tolerance=0)     # not bit for bit
    assert 0 < max(moved) < 0.01                                          # under a centimetre


def test_reproject_round_trip_off_the_map(parcels):
    back = parcels.to_crs(4326).to_crs(parcels.crs)            # the fixture at the grid origin
    assert parcels.geometry.iloc[0].hausdorff_distance(back.geometry.iloc[0]) > 80


def test_same_shape_different_start_vertex(square):
    rotated = Polygon([(100, 0), (100, 100), (0, 100), (0, 0)])
    assert square != rotated                                  # == is structural
    assert square.equals(rotated)                             # same point set
    assert square.symmetric_difference(rotated).area == pytest.approx(0, abs=1e-9)


def test_missing_crs_raises():
    with pytest.raises(ValueError, match="no CRS"):
        buffer_metres(gpd.GeoDataFrame(geometry=[Point(0, 0)]), 10)


def test_empty_input_returns_empty_not_error(parcels):
    out = buffer_metres(parcels.iloc[0:0], 10)
    assert out.empty and out.crs == parcels.crs


def test_invalid_geometry_is_counted_not_silently_dropped(messy_frame):
    clean, dropped = drop_invalid(messy_frame)
    assert dropped == 3 and len(clean) == 9

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

..........                                                               [100%]
10 passed in 0.14s

The first test makes six assertions, none of which mention a coordinate. These are the properties that break when someone reprojects wrongly, loses the CRS, or drops a column — and they keep passing when a GEOS upgrade moves a vertex in the fifteenth decimal place. The second guards the one thing a spatial-index query gets wrong most often: two overlapping zones hit the same cells, and without the set() in select_near() each would be listed twice.

The comparison tests pin down different notions of "equal". A reprojection round trip through WGS 84 does not return the same coordinates: in Oakham each vertex comes back about a millimetre from where it started. At the grid's false origin, where the 100 m square fixture sits, it comes back 86 m away. The point converts to 7.56° W, 49.77° N, off the Isles of Scilly, and on the way back PROJ finds that position outside the area of use of the OSGB36 datum shift it used going out, so it falls back to a "ballpark" transformation with no datum shift at all. A fixture with round coordinates is ideal for areas and counts, and misleading for anything that reprojects. A square written from a different starting vertex is not == — in Shapely 2 that is exact structural equality — but it is equals(), and the area of the symmetric difference is zero. Comparing that area against a tolerance is the robust way to ask "is this the same shape?".

The empty-input test earns its place: a pipeline that crashes on a day with no data is a pipeline that pages someone at 03:00 on a bank holiday.

Add one end-to-end test with real I/O

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

from src.pipeline import run


def test_pipeline_end_to_end(oakham_config):
    result = run(oakham_config)
    written = gpd.read_file(oakham_config["output"])
    assert len(written) == result["parcels"] > 0
    assert result["schools"] == 9 and result["dropped"] == 0
    assert written.crs.to_epsg() == 27700
    assert written.geometry.is_valid.all()
    assert written["INSPIREID"].is_unique


def test_degrees_mistake_is_an_error_not_a_report(oakham_config):
    with pytest.raises(ValueError, match="projected CRS"):
        run({**oakham_config, "crs": "EPSG:4326"})

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

..                                                                       [100%]
2 passed in 0.17s

One test that exercises reading, transforming, writing and reading back catches driver problems, schema surprises and permission errors that pure-function tests cannot. The second is the incident again, end to end: asked to work in degrees, the pipeline now stops instead of reporting every parcel in Oakham.

Keep the slow tests separate

%%writefile pyproject.toml
[tool.pytest.ini_options]
testpaths = ["tests"]
markers = [
    "slow: takes more than a second",
    "db: needs a live PostGIS connection",
    "integration: touches the network or real data",
]
addopts = "-q --strict-markers"

%%writefile tests/test_postgis.py
import geopandas as gpd
import pytest
from sqlalchemy import create_engine, text


@pytest.fixture(scope="session")
def postgis_engine():
    engine = create_engine("postgresql+psycopg://user:guide@localhost:55433/gis")
    try:
        with engine.connect() as conn:
            conn.execute(text("SELECT PostGIS_Version()"))
    except Exception as exc:
        pytest.skip(f"no PostGIS: {type(exc).__name__}")
    return engine


@pytest.mark.db
def test_postgis_round_trip(postgis_engine, parcels):
    parcels.to_postgis("test_parcels", postgis_engine, if_exists="replace")
    back = gpd.read_postgis("SELECT * FROM test_parcels", postgis_engine, geom_col="geometry")
    assert len(back) == len(parcels)
    assert back.crs.to_epsg() == 27700
    assert back.geometry.area.sum() == pytest.approx(parcels.geometry.area.sum())

!{sys.executable} -m pytest -p no:cacheprovider --color=no --disable-warnings -m "not db and not slow" | tail -1
!{sys.executable} -m pytest -p no:cacheprovider --color=no --disable-warnings | tail -1

12 passed, 1 deselected in 0.28s
13 passed in 0.41s

pytest -m "not db and not slow" is the fast suite you run constantly; plain pytest runs everything, in CI. The database fixture skips rather than fails when PostGIS is not there, so a laptop without Docker still gets a green fast suite — and --strict-markers turns a misspelt @pytest.mark.dbb into an error instead of a test that silently always runs.

Tests on Rutland's schools and parcels

The data used on this page

Rutland's school sites from Geofabrik's OpenStreetMap extract of 1 January 2026 (© OpenStreetMap contributors, ODbL), HM Land Registry's INSPIRE parcels for Rutland (© Crown copyright and database rights; OGL), and Oakham's three wards from the ONS December 2024 boundaries (OGL). The setup cuts the two small real fixture files the end-to-end test reads.

import sys
from pathlib import Path

import geopandas as gpd
import requests

Path("data").mkdir(exist_ok=True)
for folder in ("src", "tests", "tests/data"):
    Path(folder).mkdir(exist_ok=True)
Path("src/__init__.py").touch()
downloads = {
    "data/rutland.shp.zip": "https://download.geofabrik.de/europe/united-kingdom/england/rutland-260101-free.shp.zip",
    "data/Rutland.zip": "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip",
}
for target, url in downloads.items():
    if not Path(target).exists():
        Path(target).write_bytes(requests.get(url, 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)

schools = gpd.read_file("/vsizip/data/rutland.shp.zip/gis_osm_pois_a_free_1.shp",
                        where="fclass = 'school'", columns=["name", "fclass"])
parcels = gpd.read_file("/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml",
                        columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO")
wards = gpd.read_file("data/wards.geojson")

oakham = wards[wards["WD24NM"].str.startswith("Oakham")].union_all()
schools[schools.to_crs(27700).intersects(oakham)].to_file("tests/data/oakham_schools.gpkg")
parcels[parcels.intersects(oakham)].to_file("tests/data/oakham_parcels.gpkg")
print(f"{len(schools)} school sites, {len(parcels):,} parcels;",
      {p.name: f"{p.stat().st_size // 1024} KB" for p in Path("tests/data").iterdir()})

26 school sites, 21,253 parcels; {'oakham_schools.gpkg': '104 KB', 'oakham_parcels.gpkg': '2908 KB'}

Example 1: property-based tests for geometric invariants

%%writefile tests/test_properties.py
import geopandas as gpd
import pytest
from hypothesis import given, settings, strategies as st
from shapely.geometry import Point

from src.transforms import buffer_metres


@settings(max_examples=200, deadline=None)
@given(x=st.floats(min_value=0, max_value=700_000, allow_nan=False),
       y=st.floats(min_value=0, max_value=1_300_000, allow_nan=False),
       distance=st.floats(min_value=0.1, max_value=1000))
def test_buffer_area_grows_with_distance(x, y, distance):
    gdf = gpd.GeoDataFrame(geometry=[Point(x, y)], crs="EPSG:27700")
    small = buffer_metres(gdf, distance).geometry.area.iloc[0]
    large = buffer_metres(gdf, distance * 2).geometry.area.iloc[0]
    assert large > small
    assert large == pytest.approx(small * 4, rel=0.001)      # area scales with r²


@settings(max_examples=200, deadline=None)
@given(x=st.floats(min_value=470_000, max_value=500_000), y=st.floats(min_value=290_000, max_value=320_000))
def test_reprojection_round_trip_in_rutland(x, y):
    gdf = gpd.GeoDataFrame(geometry=[Point(x, y)], crs="EPSG:27700")
    back = gdf.to_crs(4326).to_crs(27700).geometry.iloc[0]
    assert back.distance(Point(x, y)) < 0.01

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

..                                                                       [100%]
2 passed in 0.83s

Property-based testing suits geometry unusually well, because so many operations have laws — buffering scales area quadratically, a reprojection round trip is the identity to within a tolerance, a dissolve never increases total area. Hypothesis generates 200 inputs for each law and, when one fails, shrinks it to the simplest input that still fails. The round-trip test first drew its points from the whole British National Grid, with a micrometre tolerance, and failed at once; Hypothesis shrank the failure to x=0.0, y=0.0 and a miss of 86.17 m. The law only holds where the transformation does, so the domain is now Rutland's square and the tolerance a centimetre.

Example 2: assert on the run summary, not just the files

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

from src.pipeline import run


def test_run_summary_is_honest(oakham_config):
    result = run(oakham_config)
    assert set(result) >= {"schools", "dropped", "parcels"}
    written = gpd.read_file(oakham_config["output"])
    assert len(written) == result["parcels"], "summary disagrees with the output file"

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

.                                                                        [100%]
1 passed in 0.14s

The final assertion catches a whole class of reporting bug: a summary that says 687 while the file holds 686.

Example 3: a regression test against a golden summary

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

import geopandas as gpd
import pytest

from src.pipeline import run

GOLDEN = Path(__file__).parent / "golden" / "oakham_summary.json"


def summarise(gdf) -> dict:
    return {
        "features": len(gdf),
        "crs": gdf.crs.to_string(),
        "geom_types": sorted(gdf.geom_type.unique().tolist()),
        "total_area_m2": round(float(gdf.geometry.area.sum()), 1),
        "bounds": [round(float(v), 1) for v in gdf.total_bounds],
        "columns": sorted(c for c in gdf.columns if c != "geometry"),
    }


def test_output_matches_golden(oakham_config):
    run(oakham_config)
    actual = summarise(gpd.read_file(oakham_config["output"]))
    if not GOLDEN.exists():                       # first run records the baseline
        GOLDEN.parent.mkdir(parents=True, exist_ok=True)
        GOLDEN.write_text(json.dumps(actual, indent=2))
        pytest.skip("golden file created — review and commit it")
    assert actual == json.loads(GOLDEN.read_text())

!{sys.executable} -m pytest -p no:cacheprovider --color=no --disable-warnings -rs tests/test_golden.py | tail -2
!{sys.executable} -m pytest -p no:cacheprovider --color=no --disable-warnings tests/test_golden.py | tail -1
print(Path("tests/golden/oakham_summary.json").read_text())

SKIPPED [1] tests/test_golden.py:29: golden file created — review and commit it
1 skipped in 0.17s
1 passed in 0.15s
{
  "features": 687,
  "crs": "EPSG:27700",
  "geom_types": [
    "Polygon"
  ],
  "total_area_m2": 6832596.5,
  "bounds": [
    484383.4,
    306461.8,
    487165.6,
    310279.2
  ],
  "columns": [
    "INSPIREID"
  ]
}

A golden summary — counts, areas, bounds, schema — is small enough to commit and specific enough to catch an accidental change in behaviour. Storing the summary rather than the output file keeps the repository light and the diff readable: if a change to the buffer distance or the parcel data moves the numbers, the review shows which ones.

Testing logic and geometry at once

Testing GIS code is testing two things at once: ordinary program logic, and geometric correctness. Ordinary tools handle the first. The second needs assertions phrased in the language of the domain, because the failures that matter most produce output that is structurally perfect and geographically wrong — like a list of every parcel in the county, all "near a school".

The precondition for all of it is structure. A function that reads a file, transforms it and writes another file can only be tested by giving it files, which is slow and awkward. Splitting the code into pure transformations — GeoDataFrame in, GeoDataFrame out — and a thin I/O shell means the interesting logic can be tested in milliseconds with fixtures built from a few squares. The shell then needs only one or two tests of its own.

Assertions should target properties rather than coordinates. Coordinate equality is brittle: GEOS and PROJ can change vertex output between versions, buffer approximations depend on segment counts, and reprojection is not exactly invertible — the round trip above moved every vertex by about a millimetre in Rutland, and by 86 m at the grid origin. Meanwhile the things that actually break — a lost CRS, a dropped column, a doubled row count after a join, an invalid geometry, an area that changed by 40% — are all cheap to assert and stable across library upgrades. When you must compare shapes, compare the area of their symmetric difference against a tolerance.

Fixtures should be mostly constructed, not committed. A 100 m square has an area of exactly 10,000 m², so an assertion about it reads clearly and a failure is immediately interpretable. Real data earns its place in the end-to-end test, cut down to the smallest piece that still exercises the real formats and the real quirks — here Oakham's school sites and parcels, about 3 MB rather than the county's 21,253 parcels.

Finally, speed determines whether the suite gets run. Marking database and network tests so they can be excluded keeps the everyday loop to a few seconds, which is the difference between tests that guard the pipeline and tests that everyone skips.

Strict equality, approximate areas and drifting sums

• gdf.crs == "EPSG:27700" compares a CRS object to a string: It works via pyproj's __eq__, but gdf.crs.to_epsg() == 27700 is clearer and less surprising.
• assert gdf1.equals(gdf2) is strict: It compares dtypes and index too. Prefer targeted property assertions.
• Buffer areas are approximations: A buffered point is a polygon with quad_segs × 4 vertices, so its area is slightly under πr². Use a relative tolerance.
• tmp_path is per test: Do not share written files between tests through module-level state; make a fixture instead.
• Floating-point sums drift: Comparing area.sum() across large frames needs rel=1e-9, not exact equality.
• Empty GeoDataFrames are worth an explicit CRS assertion: buffer_metres() on an empty slice kept EPSG:27700 here, but a frame built empty from scratch has no CRS unless you pass one.
• Round coordinates sit off the map: British National Grid's (0, 0) is south-west of the Isles of Scilly. Reprojecting a fixture there and back missed by 86 m; move it into the area your data covers before testing a transformation.
• Mock the network, not the geometry: Patch the download step so tests do not depend on a remote service, but let the geometry code run for real.

Internal links

• How to Validate Pipeline Inputs and Outputs Automatically in Python
• How to Chain GIS Processing Steps into a Reusable Pipeline in Python
• How to Make a GIS Workflow Reproducible in Python
• How to Run a Python GIS Pipeline in CI with GitHub Actions
• My GIS Pipeline Fails Silently: Fixing Swallowed Errors and Wrong Exit Codes
• How to Validate a GeoDataFrame Against a Schema Before Analysis

FAQ

What should I actually assert on for spatial output?

CRS, feature count, geometry types, validity, total area or length, bounding box, and the column schema. Those seven properties catch nearly every silent regression and survive library upgrades.

How do I compare two geometries in a test?

Use shapely.equals_exact(a, b, tolerance=...), or assert that the area of their symmetric difference is approximately zero. Exact == compares structure and fails on trivially different vertex ordering.

Where do I get test data?

Build it with Shapely inside fixtures — a few squares and points with round numbers. Keep at most one small real file for a single end-to-end test.

How do I test code that reads and writes files?

Use pytest's tmp_path fixture to write into a per-test directory, then read the result back and assert on it. That covers driver behaviour and schema round-tripping.

Should I test against a live PostGIS database?

Yes, but mark those tests (@pytest.mark.db) so they can be skipped locally and run in CI with a service container. Keep the fast suite free of anything that needs a network.

How do I catch the "same code, different GEOS" problem?

Assert on properties with tolerances rather than exact coordinates, and run the suite in the same container image the pipeline uses so library versions match production.

What is a golden test and is it worth it?

It compares a compact summary of the output — counts, area, bounds, schema — against a committed baseline. It is worth it for pipelines whose output should be stable, because it catches unintended behaviour changes that no unit test anticipated.