How to Test a GIS Pipeline with pytest

Problem statement

The pipeline has run every night for eight months. Last Tuesday someone changed a buffer distance from metres to feet, and nobody noticed until a planning decision came back wrong.

# the diff that did it
- buffered = parcels.to_crs(27700).buffer(50)
+ buffered = parcels.buffer(50)          # now degrees, not metres

No exception. No warning. Every downstream step ran happily on geometry that was 5,500 km across. GIS code is unusually good at failing silently — a CRS mistake produces valid 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:

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

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

from src.transforms import buffer_metres

@pytest.fixture
def points_27700():
    return gpd.GeoDataFrame(
        {"id": [1, 2]},
        geometry=[Point(325000, 674000), Point(326000, 675000)],
        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(-3.19, 55.95)], crs="EPSG:4326")
    with pytest.raises(ValueError, match="projected CRS"):
        buffer_metres(gdf, 50)

The second test is the one that would have caught the incident: it asserts that the function refuses to buffer in degrees at all.

What to test at each level

Layered test pyramid from pure transform tests through I/O tests to end-to-end runs. — Many fast transform tests, a few I/O tests, one end-to-end run.

Step-by-step solution

Vertical steps: extract pure functions, build fixtures, assert properties, add I/O tests, run in CI. — Five steps — the first one is what makes the other four possible.

Make the pipeline testable first

A script that reads, transforms and writes in one function can only be tested end to end. Split it.

# 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}")
    out = gdf.copy()
    out["geometry"] = gdf.geometry.buffer(distance_m)
    return out

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

# src/pipeline.py — the impure shell: reads, calls transforms, writes
def run(config):
    gdf = gpd.read_file(config["input"])
    gdf, dropped = drop_invalid(gdf)
    gdf = buffer_metres(gdf.to_crs(config["crs"]), config["buffer_m"])
    gdf.to_file(config["output"], driver="GPKG")
    return {"features": len(gdf), "dropped": dropped}

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 a 200 MB shapefile to test a clip is how test suites become unusable. Construct exactly what the test needs.

# tests/conftest.py
import geopandas as gpd
import pytest
from shapely.geometry import Point, Polygon, LineString

@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="EPSG:27700",
    )

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

@pytest.fixture
def parcels_file(parcels, tmp_path):
    path = tmp_path / "parcels.gpkg"
    parcels.to_file(path, layer="parcels", driver="GPKG")
    return path

A square of side 100 has an area of exactly 10,000 — which makes assertions readable and failures obvious. tmp_path gives each test its own directory, cleaned up automatically.

Assert on properties, not on coordinates

def test_clip_keeps_only_overlapping_parcels(parcels, square):
    boundary = gpd.GeoDataFrame(geometry=[square.buffer(10)], crs=parcels.crs)
    clipped = clip_to(parcels, boundary)

    assert len(clipped) == 1                                  # count
    assert clipped.crs == parcels.crs                         # CRS preserved
    assert clipped.geometry.is_valid.all()                    # validity
    assert (clipped.geom_type == "Polygon").all()             # type
    assert clipped.geometry.area.sum() <= parcels.geometry.area.sum()   # invariant
    assert set(clipped.columns) == set(parcels.columns)       # schema

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 changes a vertex in the fifteenth decimal place.

Compare geometry with a tolerance

import pytest
from shapely.geometry import Point
from shapely import equals_exact

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(parcels):
    there_and_back = parcels.to_crs(4326).to_crs(parcels.crs)
    for a, b in zip(parcels.geometry, there_and_back.geometry):
        assert equals_exact(a, b, tolerance=1e-6)

def test_geometry_matches_expected(parcels):
    expected = Point(50, 50).buffer(10)
    actual = Point(50, 50).buffer(10, quad_segs=8)
    # symmetric difference area is the robust way to compare shapes
    assert actual.symmetric_difference(expected).area == pytest.approx(0, abs=1e-9)

== on geometries tests structural equality, which fails on a differently ordered but identical ring. Comparing the area of the symmetric difference is the tolerant, meaningful check.

Test the failure paths explicitly

import pytest
import geopandas as gpd
from shapely.geometry import Point

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

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

def test_invalid_geometry_is_counted_not_silently_dropped(parcels, bowtie):
    mixed = gpd.GeoDataFrame(geometry=list(parcels.geometry) + [bowtie], crs=parcels.crs)
    clean, dropped = drop_invalid(mixed)
    assert dropped == 1 and len(clean) == 3

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

def test_pipeline_end_to_end(parcels_file, tmp_path):
    out = tmp_path / "out.gpkg"
    result = run({"input": str(parcels_file), "output": str(out),
                  "crs": "EPSG:27700", "buffer_m": 25})

    assert out.exists()
    written = gpd.read_file(out)
    assert len(written) == result["features"] == 3
    assert written.crs.to_epsg() == 27700
    assert written.geometry.is_valid.all()
    assert written.geometry.area.min() > 0

One test that exercises reading, transforming, writing and reading back catches driver problems, schema surprises and permission errors that pure-function tests cannot.

Keep the slow tests separate

# 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"

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

pytest -m "not db and not slow"      # the fast suite you run constantly
pytest                               # everything, in CI

Code examples

Example 1: a fixture library worth having

# tests/conftest.py
import geopandas as gpd
import numpy as np
import pytest
from shapely.geometry import Point, Polygon, LineString, MultiPolygon

BNG = "EPSG:27700"

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

@pytest.fixture
def random_points():
    rng = np.random.default_rng(42)               # seeded: reproducible failures
    xs = rng.uniform(0, 300, 50)
    ys = 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."""
    rows = list(grid_polygons.geometry) + [bowtie, None, Polygon()]
    return gpd.GeoDataFrame({"id": range(len(rows))}, geometry=rows, crs=BNG)

@pytest.fixture
def bowtie():
    return Polygon([(0, 0), (10, 10), (10, 0), (0, 10)])

A seeded random generator matters: an unseeded fixture produces a test that fails once a fortnight and cannot be reproduced.

Example 2: property-based tests for geometric invariants

from hypothesis import given, strategies as st
import geopandas as gpd
from shapely.geometry import Point

@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²

Property-based testing suits geometry unusually well, because so many operations have laws — buffering scales area quadratically, a reprojection round trip is the identity, a dissolve never increases total area.

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

def test_run_summary_is_honest(parcels_file, tmp_path):
    result = run({"input": str(parcels_file), "output": str(tmp_path / "o.gpkg"),
                  "crs": "EPSG:27700", "buffer_m": 10})

    assert set(result) >= {"features", "dropped"}
    assert result["features"] == 3
    assert result["dropped"] == 0

    written = gpd.read_file(tmp_path / "o.gpkg")
    assert len(written) == result["features"], "summary disagrees with the output file"

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

Example 4: a regression test against a golden result

import json
from pathlib import Path
import geopandas as gpd
import pytest

GOLDEN = Path(__file__).parent / "golden" / "parcels_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()), 3),
        "bounds": [round(float(v), 3) for v in gdf.total_bounds],
        "columns": sorted(c for c in gdf.columns if c != "geometry"),
    }

def test_output_matches_golden(parcels_file, tmp_path):
    out = tmp_path / "out.gpkg"
    run({"input": str(parcels_file), "output": str(out),
         "crs": "EPSG:27700", "buffer_m": 25})
    actual = summarise(gpd.read_file(out))

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

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.

Explanation

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.

Checklist of geometric assertions: CRS, count, type, validity, area, bounds, schema. — Seven properties that catch nearly every silent spatial regression.

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 four points. 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 change vertex output between versions, buffer approximations depend on segment counts, and reprojection is not exactly invertible. 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 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 datasets bring bloat, licensing questions and irrelevant complexity; keep one small real file for the end-to-end test and build everything else with Shapely.

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

Edge cases or notes

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 lose their CRS easily: Assert on it explicitly; several operations return an empty frame with crs=None.
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

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.

How to Test a GIS Pipeline with pytest #

Problem statement #

Quick answer #

What to test at each level #

Step-by-step solution #

Make the pipeline testable first #

Build small fixtures in code #

Assert on properties, not on coordinates #

Compare geometry with a tolerance #

Test the failure paths explicitly #

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

Keep the slow tests separate #

Code examples #

Example 1: a fixture library worth having #

Example 2: property-based tests for geometric invariants #

Example 3: assert on the run summary, not just the files #

Example 4: a regression test against a golden result #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

What should I actually assert on for spatial output? #

How do I compare two geometries in a test? #

Where do I get test data? #

How do I test code that reads and writes files? #

Should I test against a live PostGIS database? #

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

What is a golden test and is it worth it? #