βš™οΈ How-To Testing, CI & Docker Geometry & topology Updated

How to Assert on Geometry in pytest Without Flaky Comparisons

Problem statement

The test is right and it fails:

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

Or the geometry is correct and the equality is not:

def test_reproject_round_trip():
    original = Point(-2.2426, 53.4808)
    there_and_back = reproject(reproject(original, 4326, 27700), 27700, 4326)
    assert there_and_back == original
assert POINT (-2.2426000000000004 53.480799999999995) == POINT (-2.2426 53.4808)

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)])
assert a == b          # False β€” same polygon, different vertex order

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)
Grid comparing ==, equals, equals_exact and area comparison against what each is sensitive to.
Four kinds of "the same". Each answers a different question.
Assertion True when Sensitive to
a == b identical structure and coordinates vertex order, float noise, ring direction
a.equals(b) the same point set nothing but the shape
equals_exact(a, b, tol) same vertices within tol vertex order and count
a.symmetric_difference(b).area < tol the shapes agree to within an area neither
a.area == approx(x) one measured property only that property

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

Step-by-step solution

1. Never compare floats with == 

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

buffer approximates a circle with straight segments β€” 8 per quadrant by default β€” so the area is not Ο€ rΒ². The number depends on quad_segs, and that is a documented, deliberate approximation rather than an error.

# βœ… 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(383_618.507, abs=0.001)     # within 1 mm
assert gdf.area.sum() == pytest.approx(8_412_993.4, 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

Panels showing two polygons that are the same shape but differ in vertex order, start point and ring direction.
Three representations of one square. `==` says they are different; `.equals()` says they are not.
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))
geometry mismatch (dissolved ward):
  type      Polygon vs MultiPolygon
  valid     True vs True
  empty     False vs False
  area      8412.993211 vs 8412.993211 (Ξ” +0.000000)
  vertices  84 vs 86
  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)
  first diff at POINT EMPTY

Identical area, identical bounds, empty symmetric difference β€” the shapes are the same and only the type differs, Polygon versus MultiPolygon. 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.

Code examples

Example 1: a reusable assertion module 

# tests/geometry_assertions.py
import math
import pytest
import geopandas as gpd
from shapely import equals_exact, get_num_coordinates
from shapely.geometry.base import BaseGeometry

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):
        from shapely.validation import explain_validity
        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 ""))

__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 (Ancoats)" is the difference between a two-minute diagnosis and a twenty-minute one.

Example 2: testing operations that legitimately change representation 

import pytest
import geopandas as gpd
from shapely.geometry import Polygon, MultiPolygon, 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:

import pytest
from hypothesis import given, settings, strategies as st, assume
from shapely.geometry import Polygon, box
from shapely.ops import unary_union

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)
tests/test_properties.py::test_inclusion_exclusion FAILED

Falsifying example: test_inclusion_exclusion(
    a=<POLYGON ((0.1 0.1, ...))>,
    b=<POLYGON ((0.1 0.1, ...))>,
)

Hypothesis generates hundreds of cases and, on failure, shrinks them to the smallest input that still fails β€” which is usually a degenerate case you would never have written by hand: a zero-width rectangle, coordinates differing by one float ulp, shapes sharing exactly one vertex.

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.

Explanation

Scene showing two nearly identical polygons and the thin symmetric difference between them.
The symmetric difference turns "are these the same?" into a number you can set a threshold on.

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.

Edge cases or notes

  • 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.

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 β€” 8 per quadrant by default. 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.