How to Assert on Geometry in pytest Without Flaky Comparisons

Problem statement

The test is right and it fails:

import math

from shapely.geometry import Point, Polygon


def test_buffer_area():
    square = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
    assert square.buffer(5).area == 100 + 4 * 50 + math.pi * 25


try:
    test_buffer_area()
except AssertionError:
    square = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
    print(f"assert {square.buffer(5).area!r} == {100 + 4 * 50 + math.pi * 25!r}")

assert 378.4137122636484 == 378.53981633974485

Or the geometry is correct and the equality is not:

from pyproj import Transformer
from shapely.ops import transform


def reproject(geom, src, dst):
    return transform(Transformer.from_crs(src, dst, always_xy=True).transform, geom)


def test_reproject_round_trip():
    original = Point(-1.5474, 53.7952)                  # Leeds station
    there_and_back = reproject(reproject(original, 4326, 27700), 27700, 4326)
    assert there_and_back == original, f"{there_and_back} != {original}"


try:
    test_reproject_round_trip()
except AssertionError as err:
    print(f"AssertionError: {err}")

AssertionError: POINT (-1.5473999870586006 53.79519999733061) != POINT (-1.5474 53.7952)

Or the shapes are identical and the objects are not:

a = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
b = Polygon([(10, 10), (0, 10), (0, 0), (10, 0)])
print(a == b, "— same square, different starting vertex")

False — same square, different starting vertex

Geometry is floating-point data with several valid representations of the same shape. Exact equality is almost never the assertion you want, and the failure message from a naive comparison — two thousand-character WKT strings — tells you nothing about what differed.

Quick answer

Assert on properties, with tolerance, and use the right equality for the question:

import pytest
from shapely import equals_exact
from shapely.testing import assert_geometries_equal

# a measurement — always with tolerance
assert result.area == pytest.approx(10_000, rel=1e-9)

# "the same shape", tolerant of vertex order and float noise
assert result.equals(expected)

# "the same coordinates", within a tolerance
assert equals_exact(result, expected, tolerance=1e-6)

# a GeoSeries or array of geometries
assert_geometries_equal(result.geometry.values, expected.geometry.values, tolerance=1e-6)

Prefer .equals() unless the test is specifically about vertex representation.

False
378.4137122636484

No float ==, the right equality, area tolerances

1. Never compare floats with ==

print(0.1 + 0.2 == 0.3)
print(Polygon([(0, 0), (10, 0), (10, 10), (0, 10)]).buffer(5).area)

False
378.4137122636484

buffer approximates a circle with straight segments — 16 per quadrant by default in Shapely 2 — so the area is not π r². The number depends on quad_segs, and that is a documented, deliberate approximation rather than an error.

import pytest

square = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
buffered = square.buffer(5)

# ✅ relative tolerance for a measurement
assert buffered.area == pytest.approx(378.45, rel=1e-3)

# ✅ or assert the property you actually care about
assert buffered.area > square.area
assert buffered.contains(square)

pytest.approx takes rel (relative) or abs (absolute). Use rel for areas and distances, whose magnitude varies; abs for coordinates, where you want "within a millimetre" regardless of the coordinate values:

assert point.x == pytest.approx(429_911.504, abs=0.001)     # within 1 mm
assert gdf.area.sum() == pytest.approx(551_657_398.7, rel=1e-6)

A relative tolerance on a coordinate near zero is a trap: pytest.approx(0.0, rel=1e-9) accepts only exactly zero, because the tolerance is a fraction of the expected value.

2. Choose the right notion of equality

from shapely.geometry import Polygon
from shapely import equals_exact

a = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
b = Polygon([(10, 10), (0, 10), (0, 0), (10, 0)])       # rotated start point
c = Polygon([(0, 0), (0, 10), (10, 10), (10, 0)])       # reversed winding
d = Polygon([(0, 0), (5, 0), (10, 0), (10, 10), (0, 10)])  # extra collinear vertex

for name, g in [("rotated", b), ("reversed", c), ("collinear", d)]:
    print(f"{name:<10} ==  {a == g}   equals  {a.equals(g)}   "
          f"exact  {equals_exact(a, g, 1e-9)}")

rotated    ==  False   equals  True   exact  False
reversed   ==  False   equals  True   exact  False
collinear  ==  False   equals  True   exact  False

All three are the same square. == and equals_exact compare representation; .equals() compares the point set — the actual region of the plane covered.

Any operation that goes through GEOS may return a different but equivalent representation. union, intersection, buffer, make_valid and a GeoPackage round-trip can all reorder vertices or change the start point. So a test asserting == on the output of any of those is asserting an implementation detail.

Use equals_exact only when the test is about vertices — that simplification removed the right ones, that a snapping operation moved a coordinate to a specific place.

3. Use area-based tolerance for shapes that should nearly match

.equals() is exact about the point set, so a coordinate moved by one nanometre makes it False. For "these should be the same shape, to within a millimetre", compare the symmetric difference:

def assert_shapes_match(actual, expected, *, tolerance_m2=1e-6):
    """The two shapes agree except for at most `tolerance_m2` of area."""
    diff = actual.symmetric_difference(expected)
    assert diff.area <= tolerance_m2, (
        f"shapes differ by {diff.area:.9f} m² "
        f"(actual {actual.area:.4f}, expected {expected.area:.4f})")

def assert_shapes_match_relative(actual, expected, *, max_fraction=1e-9):
    diff = actual.symmetric_difference(expected).area
    base = max(actual.area, expected.area, 1e-12)
    assert diff / base <= max_fraction, (
        f"shapes differ by {100 * diff / base:.6f}% of area")

The symmetric difference is the area in one shape or the other but not both, so it is exactly "how much do these disagree". It is the most robust geometry assertion available: insensitive to vertex order, vertex count, winding and float noise, and it fails with a number you can interpret.

4. Compare whole GeoDataFrames on the properties that matter

import geopandas as gpd
import pandas as pd
from geopandas.testing import assert_geodataframe_equal

# strict: same rows, same order, same dtypes, same CRS, same geometry
assert_geodataframe_equal(result, expected)

# tolerant of the things that do not matter
assert_geodataframe_equal(
    result.sort_values("id").reset_index(drop=True),
    expected.sort_values("id").reset_index(drop=True),
    check_dtype=False,
    check_less_precise=True,      # coordinate comparison at reduced precision
    check_crs=True,
)

assert_geodataframe_equal gives a useful failure message that names the differing column and row, which a hand-rolled comparison does not.

For most tests, though, asserting a few properties is better than asserting the whole frame:

def test_cleaning_preserves_the_things_it_should(raw):
    result = clean(raw)
    assert result.crs == raw.crs
    assert set(result.columns) == set(raw.columns)
    assert len(result) <= len(raw), "cleaning must not invent rows"
    assert result.geometry.is_valid.all()
    assert result.geometry.notna().all()
    assert result.area.sum() <= raw.area.sum() * 1.001

Six assertions, each naming a property of the function. None of them breaks when the input data is updated — the failure described in test fixtures for GIS code.

5. Make failures readable

A failed geometry assertion that prints two WKT strings is useless. Write helpers that report the difference:

def describe_difference(actual, expected, name=""):
    """A human-readable summary of how two geometries differ."""
    lines = [f"geometry mismatch{f' ({name})' if name else ''}:"]
    lines.append(f"  type      {actual.geom_type} vs {expected.geom_type}")
    lines.append(f"  valid     {actual.is_valid} vs {expected.is_valid}")
    lines.append(f"  empty     {actual.is_empty} vs {expected.is_empty}")
    if not (actual.is_empty or expected.is_empty):
        lines.append(f"  area      {actual.area:.6f} vs {expected.area:.6f} "
                     f"(Δ {actual.area - expected.area:+.6f})")
        lines.append(f"  vertices  {count_vertices(actual)} vs {count_vertices(expected)}")
        lines.append(f"  bounds    {tuple(round(v, 4) for v in actual.bounds)}")
        lines.append(f"            {tuple(round(v, 4) for v in expected.bounds)}")
        sym = actual.symmetric_difference(expected)
        lines.append(f"  sym diff  {sym.area:.9f} "
                     f"({100 * sym.area / max(expected.area, 1e-12):.6f}% of expected)")
        if not sym.is_empty:
            lines.append(f"  first diff at {sym.representative_point().wkt}")
    return "\n".join(lines)

def count_vertices(geom):
    from shapely import get_num_coordinates
    return int(get_num_coordinates(geom))


# a dissolve that returned a MultiPolygon where a Polygon was expected
from shapely.geometry import MultiPolygon, box

print(describe_difference(MultiPolygon([box(0, 0, 100, 100)]), box(0, 0, 100, 100), "dissolved ward"))

geometry mismatch (dissolved ward):
  type      MultiPolygon vs Polygon
  valid     True vs True
  empty     False vs False
  area      10000.000000 vs 10000.000000 (Δ +0.000000)
  vertices  5 vs 5
  bounds    (0.0, 0.0, 100.0, 100.0)
            (0.0, 0.0, 100.0, 100.0)
  sym diff  0.000000000 (0.000000% of expected)

Identical area, identical bounds, empty symmetric difference — the shapes are the same and only the type differs, MultiPolygon versus Polygon. That is a real and common difference after a dissolve, and it takes seconds to diagnose from this output and a long time from two WKT dumps.

An assertion module and property-based tests

The data used on this page

None — deliberately. Tests should build the geometry they assert on, so the examples here construct squares and boxes whose true areas are known exactly. The only real coordinate is Leeds station's, in the round-trip test. The setup makes the tests/ folder the examples write into.

import sys
from pathlib import Path

Path("tests").mkdir(exist_ok=True)
Path("tests/__init__.py").touch()

Example 1: a reusable assertion module

%%writefile tests/geometry_assertions.py
import pytest
from shapely import get_num_coordinates
from shapely.geometry.base import BaseGeometry
from shapely.validation import explain_validity


def describe_difference(actual, expected, name=""):
    lines = [f"geometry mismatch{f' ({name})' if name else ''}:",
             f"  type      {actual.geom_type} vs {expected.geom_type}",
             f"  area      {actual.area:.6f} vs {expected.area:.6f}",
             f"  vertices  {get_num_coordinates(actual)} vs {get_num_coordinates(expected)}"]
    sym = actual.symmetric_difference(expected)
    lines.append(f"  sym diff  {sym.area:.9f}")
    return "\n".join(lines)


def assert_same_shape(actual, expected, *, tolerance=1e-9, relative=True, name=""):
    """The two geometries cover the same region, within a tolerance."""
    __tracebackhide__ = True
    if actual is None or expected is None:
        assert actual is expected, f"{name}: one geometry is None"
        return
    if actual.is_empty or expected.is_empty:
        assert actual.is_empty == expected.is_empty, (
            f"{name}: one geometry is empty and the other is not")
        return

    diff = actual.symmetric_difference(expected).area
    limit = (tolerance * max(actual.area, expected.area, 1e-12)
             if relative else tolerance)
    if diff > limit:
        pytest.fail(describe_difference(actual, expected, name))

def assert_area(geom, expected_m2, *, rel=1e-9, name=""):
    __tracebackhide__ = True
    assert geom.area == pytest.approx(expected_m2, rel=rel), (
        f"{name}: area {geom.area:.6f} != expected {expected_m2:.6f}")

def assert_within_distance(a, b, max_m, *, name=""):
    __tracebackhide__ = True
    d = a.distance(b)
    assert d <= max_m, f"{name}: {d:.6f} m apart, expected at most {max_m} m"

def assert_valid(gdf_or_geom, *, name=""):
    __tracebackhide__ = True
    if isinstance(gdf_or_geom, BaseGeometry):
        assert gdf_or_geom.is_valid, f"{name}: invalid — {explain_validity(gdf_or_geom)}"
        return
    invalid = gdf_or_geom[~gdf_or_geom.geometry.is_valid]
    if len(invalid):
        detail = "; ".join(f"row {i}: {explain_validity(g)}"
                           for i, g in invalid.geometry.head(5).items())
        pytest.fail(f"{name}: {len(invalid)} of {len(gdf_or_geom)} invalid — {detail}")

def assert_frames_match(actual, expected, *, key=None, tolerance=1e-9,
                        check_crs=True, columns=None):
    """Compare two GeoDataFrames row by row, ignoring row order."""
    __tracebackhide__ = True
    assert len(actual) == len(expected), (
        f"{len(actual)} rows vs {len(expected)} expected")
    if check_crs:
        assert actual.crs == expected.crs, f"CRS {actual.crs} vs {expected.crs}"

    if key:
        a = actual.sort_values(key).reset_index(drop=True)
        e = expected.sort_values(key).reset_index(drop=True)
    else:
        a, e = actual.reset_index(drop=True), expected.reset_index(drop=True)

    for col in (columns or [c for c in e.columns if c != e.geometry.name]):
        assert col in a.columns, f"missing column '{col}'"
        import pandas.testing as pdt
        pdt.assert_series_equal(a[col], e[col], check_dtype=False,
                                obj=f"column '{col}'")

    for i, (ga, ge) in enumerate(zip(a.geometry, e.geometry)):
        assert_same_shape(ga, ge, tolerance=tolerance,
                          name=f"row {i}" + (f" ({a.loc[i, key]})" if key else ""))

%%writefile is the notebook's way of saving a cell as a file, so the helpers become an importable module the tests below can share. The earlier version of this module used explain_validity before importing it and called a describe_difference it never defined — both invisible until an assertion failed, which is exactly when you need the helper to work.

__tracebackhide__ = True is the detail that makes these pleasant to use: pytest then reports the failure at the test's call site rather than inside the helper, so the traceback points at the line you wrote.

Every helper takes a name and puts it in the message, because in a loop over twenty features "row 14 (Headingley & Hyde Park)" is the difference between a two-minute diagnosis and a twenty-minute one.

Example 2: testing operations that legitimately change representation

%%writefile tests/test_representation.py
import pytest
import geopandas as gpd
from shapely.geometry import box
from tests.geometry_assertions import assert_same_shape, assert_area, assert_valid

@pytest.fixture
def two_squares():
    """Two adjacent 10×10 squares sharing an edge — total area 200 m²."""
    return gpd.GeoDataFrame(
        {"group": ["a", "a"], "id": [1, 2]},
        geometry=[box(0, 0, 10, 10), box(10, 0, 20, 10)],
        crs=27700)

def test_dissolve_merges_adjacent_squares(two_squares):
    result = two_squares.dissolve(by="group")
    assert len(result) == 1
    # area is conserved exactly — this is the assertion that matters
    assert_area(result.geometry.iloc[0], 200.0, rel=1e-12, name="dissolved")
    # the shape is one 20×10 rectangle, however GEOS chooses to represent it
    assert_same_shape(result.geometry.iloc[0], box(0, 0, 20, 10), name="dissolved")

def test_dissolve_may_return_multipolygon(two_squares):
    """Non-adjacent inputs dissolve to a MultiPolygon — do not assert Polygon."""
    apart = two_squares.copy()
    apart.loc[1, "geometry"] = box(50, 0, 60, 10)
    result = apart.dissolve(by="group")
    assert result.geometry.iloc[0].geom_type in {"Polygon", "MultiPolygon"}
    assert_area(result.geometry.iloc[0], 200.0, rel=1e-12)

def test_simplify_preserves_area_approximately(two_squares):
    merged = two_squares.dissolve(by="group").geometry.iloc[0]
    simplified = merged.simplify(0.5, preserve_topology=True)
    assert_valid(simplified, name="simplified")
    # simplification is lossy — assert a bound, not equality
    assert simplified.area == pytest.approx(merged.area, rel=0.05)
    assert simplified.area <= merged.area * 1.05

def test_round_trip_through_geopackage(two_squares, tmp_path):
    path = tmp_path / "round.gpkg"
    two_squares.to_file(path, driver="GPKG")
    back = gpd.read_file(path)
    assert back.crs == two_squares.crs
    for original, restored in zip(two_squares.geometry, back.geometry):
        # a file round trip can reorder vertices; the shape must survive
        assert_same_shape(original, restored, tolerance=1e-9)

The second test is the instructive one. dissolve returns a Polygon when the parts touch and a MultiPolygon when they do not, so assert result.geom_type == "Polygon" passes on one fixture and fails on another with no bug involved. Asserting the area — which is conserved either way — tests the behaviour that matters.

The round-trip test asserts shape rather than == because writing to GeoPackage and reading back goes through WKB and GDAL, either of which may normalise vertex order or ring direction. That normalisation is correct, and a test that fails on it is testing the file format rather than your code.

Example 3: property-based testing with Hypothesis

For operations with invariants that should hold for any input, generate the input:

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

coords = st.floats(min_value=-1000, max_value=1000,
                   allow_nan=False, allow_infinity=False, width=32)

@st.composite
def rectangles(draw):
    x0 = draw(coords)
    y0 = draw(coords)
    w = draw(st.floats(min_value=0.1, max_value=500, allow_nan=False))
    h = draw(st.floats(min_value=0.1, max_value=500, allow_nan=False))
    return box(x0, y0, x0 + w, y0 + h)

@given(a=rectangles(), b=rectangles())
@settings(max_examples=200, deadline=None)
def test_union_area_is_bounded(a, b):
    """|A ∪ B| is at least max(|A|,|B|) and at most |A| + |B|, always."""
    u = a.union(b)
    assert u.area >= max(a.area, b.area) - 1e-9
    assert u.area <= a.area + b.area + 1e-9

@given(a=rectangles(), b=rectangles())
@settings(max_examples=200, deadline=None)
def test_inclusion_exclusion(a, b):
    """|A ∪ B| + |A ∩ B| == |A| + |B| — a law, for any two shapes."""
    assert a.union(b).area + a.intersection(b).area == pytest.approx(
        a.area + b.area, rel=1e-9, abs=1e-9)

@given(g=rectangles(), d=st.floats(min_value=0.01, max_value=100, allow_nan=False))
@settings(max_examples=100, deadline=None)
def test_buffer_contains_original(g, d):
    """A positive buffer always contains what it buffered."""
    assert g.buffer(d).contains(g)

Run all three files the way CI would:

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

.......                                                                  [100%]
7 passed in 0.90s

All seven tests pass, and the three properties held for every one of the few hundred rectangles Hypothesis generated. That is the useful outcome, and it is not guaranteed: when a property does fail, Hypothesis shrinks the input to the smallest case that still fails — usually a degenerate one you would never have written by hand, such as rectangles sharing exactly one vertex or coordinates one float step apart. The tolerances in test_inclusion_exclusion are what keep float noise in union and intersection from being reported as a broken law.

Inclusion–exclusion is a mathematical identity, so it holds for every pair of shapes. Assertions of that kind are ideal for property-based testing, because there is no expected value to compute. Note both rel and abs in the tolerance: near-zero areas need the absolute term, since a relative tolerance around zero admits almost nothing.

deadline=None is necessary because GEOS operations on large shapes can exceed Hypothesis's default per-example time limit, which would otherwise produce flaky failures unrelated to correctness.

Floating point and representation, compounded

Geometry assertions are hard for two reasons, and they compound.

The first is floating point. Coordinates are doubles, and every geometric operation is arithmetic on them. Intersections compute new coordinates that are almost never exactly representable; a union of two shapes with a shared edge may produce a vertex a few ulps off that edge; a round trip through a projection and back does not land on the original number. This is not a defect in GEOS but a consequence of computing with finite precision — the subject of coordinate precision and floating point in GIS. Exact equality on any computed geometry is therefore an assertion about arithmetic luck.

The second is representational freedom. A polygon is defined by the region it covers, but it is stored as an ordered ring of vertices, and many rings describe one region. The start vertex can be any vertex; the winding can be clockwise or counter-clockwise; collinear vertices may be present or absent; a shape may be a Polygon or a single-part MultiPolygon. GEOS makes no promise about which representation an operation returns, and it is free to change between versions. So ==, which compares representation, asserts something your code did not decide.

.equals() resolves the second problem by comparing point sets — it asks whether the two geometries cover the same region, and answers correctly regardless of representation. It does not resolve the first, because a coordinate off by one nanometre genuinely changes the region by a nanometre. This is why the symmetric difference area is the most robust assertion available: it is representation-independent and it expresses the disagreement as a magnitude you can set a threshold on. "These shapes differ by less than a square millimetre" is a claim about the world; "these shapes are identical objects" is a claim about memory.

Property assertions are stronger than value assertions, and the reason is what they survive. assert result.area == approx(8412.99) encodes one number from one run; it breaks when the fixture changes, when GEOS improves, when the buffer's segment count changes. assert result.area <= input.area encodes a fact about the function — that clipping cannot enlarge — and stays true across all of those. Conservation laws, monotonicity, containment and idempotency are all available as assertions, and they usually catch more bugs than an exact value would, because they hold for every input rather than one.

Finally, note what a good failure message is worth. Two WKT strings tell you the geometries differ. A message reporting that the areas match to twelve digits, the bounds are identical, and only the geometry type differs tells you the answer immediately — and the difference between those two experiences, multiplied over a test suite's lifetime, is most of what makes a suite something people trust rather than something they re-run.

approx near zero, equals against ==, vertex order

• pytest.approx(0.0, rel=...) accepts only exactly zero. Use abs= for values near zero.
• .equals() compares point sets; == compares structure. Almost always you want the first.
• equals_exact(a, b, tolerance) still requires the same vertex order and count — it is a tolerant ==, not a tolerant .equals().
• dissolve and unary_union return Polygon or MultiPolygon depending on adjacency. Do not assert the type.
• buffer approximates circles with quad_segs segments per quadrant, so its area is not π r².
• assert_geodataframe_equal checks row order unless you sort first, and check_less_precise=True relaxes coordinates.
• shapely.testing.assert_geometries_equal handles arrays of geometries with a tolerance.
• __tracebackhide__ = True in a helper makes pytest report at the caller, which matters a great deal in practice.
• explain_validity(geom) turns "invalid" into "Self-intersection at or near point (5, 5)".
• Hypothesis needs deadline=None for GEOS work, or slow examples fail spuriously.

Internal links

• How to test a GIS pipeline with pytest — the surrounding setup
• Test fixtures for GIS code — the data these assertions run against
• What to test in a GIS pipeline (and what you cannot) — choosing what to assert
• Coordinate precision and floating point in GIS explained — why exact equality fails
• What makes a geometry valid? The OGC rules explained — what is_valid is checking
• How to test a GIS script that reads and writes files — round-trip assertions
• How to fix invalid geometries in Python — what a validity assertion should point you to
• GIS tests pass locally but fail in CI — when a tolerance is too tight for a different GEOS

FAQ

Why does my geometry equality test fail when the shapes look identical?

== compares structure — vertex order, start point, winding — not the region covered. Use .equals(), which compares point sets, or compare the symmetric difference area.

What tolerance should I use?

For coordinates in metres, abs=1e-6 (a micrometre) is safely below any real difference. For areas, rel=1e-9. For anything that has been through a projection round trip, loosen to about a millimetre.

Should I use equals or equals_exact?

.equals() for "the same shape", which is almost always the question. equals_exact(a, b, tol) only when the test is specifically about vertex positions, such as verifying a snap or a simplification.

How do I compare two GeoDataFrames?

assert_geodataframe_equal from geopandas.testing for a strict comparison, after sorting to remove row-order dependence. For most tests, asserting a few properties is more durable.

Why does buffer(5).area not equal π r² plus the square?

buffer approximates the circular arcs with straight segments — 16 per quadrant by default in Shapely 2, which is 378.41 m² for the square above against 378.54 for the exact shape. The result is deliberately slightly smaller. Increase quad_segs or assert with tolerance.

Should I assert exact areas?

Assert properties instead where you can: area conserved by a dissolve, area not increased by a clip, area within 5% after a simplify. Those survive data changes and library updates; a copied number does not.

Is property-based testing worth it for geometry?

For operations with mathematical invariants — inclusion–exclusion, containment after buffering, idempotency of cleaning — yes. Hypothesis shrinks failures to minimal cases that expose degenerate inputs you would not have thought to write.