Spatial Epi Public Health GIS

Writing Geometry Validity Assertions for a Pipeline

This guide, part of Spatial Data Quality Monitoring, implements the validity family of checks — the assertions that stop structurally broken geometry before it reaches an algorithm that will either crash on it or, worse, quietly compute something slightly wrong.

Problem Context & Constraints

Invalid geometry is common in public health reference data. Boundary files edited by hand acquire self-intersections where two edits met; polygons converted between formats lose ring orientation; tiny slivers appear where two jurisdictions’ boundaries were digitised separately; multipart features arrive where single parts were expected.

The downstream consequences split into two groups. Loud failures — an overlay library raising a topology exception — are annoying and safe, because the run stops. Silent failures are the dangerous ones: a spatial weights build that treats a bow-tie polygon as two disconnected parts and reports an island that is not one; an area computation that returns a negative value for a reversed ring and quietly subtracts from a sum; a point-in-polygon test that gives inconsistent answers near a self-intersection.

Three constraints shape the assertions. Validity is necessary and not sufficient — a valid geometry of the wrong type, or outside the jurisdiction, or overlapping its neighbour, passes a validity check and still breaks the analysis. Repair must be explicit — automatic repair functions change geometry, sometimes substantially, and a pipeline that repairs silently has replaced the problem with an unrecorded modification. And the check must run where geometry enters, including reference layers, not only on case points.

Four Geometry Defects and What They Silently Break Four panels. A bow-tie polygon, self-intersecting at a single point, causes a spatial weights build to report two disconnected parts. A sliver gap between two adjacent polygons causes contiguity weights to miss a neighbour relationship. An overlapping pair of polygons causes an areal aggregation to count a region twice. A polygon with a reversed ring causes an area computation to return a negative value that subtracts from a total. Each panel names the check that catches it. Each defect passes some checks and silently breaks something downstream Bow-tie weights see two parts caught by: is_valid Sliver gap neighbours not adjacent caught by: gap check Overlap area counted twice caught by: overlap check Reversed ring signed area negative caught by: orientation A validity check alone catches only the first of the four

Prerequisites

  • python 3.11
  • geopandas 0.14.4, shapely 2.0.4, pandas 2.2.2, numpy 1.26.4
  • The layer to be checked, in the pipeline’s projected metric coordinate reference system
  • The jurisdiction polygon for containment
  • A declared expected geometry type per layer and, for areal layers, whether they are meant to tessellate

Step-by-Step Solution

Assert in layers — validity, type, emptiness, containment, then topology for areal layers — and repair only by rules that record what they changed.

# Layered geometry assertions with explicit, logged repair.
# Pinned: geopandas==0.14.4, shapely==2.0.4, pandas==2.2.2, numpy==1.26.4
import logging

import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
from shapely.validation import explain_validity, make_valid

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("dq.geometry")

SLIVER_M2 = 25.0             # overlaps or gaps smaller than this are digitising noise
MAX_REPAIR_AREA_CHANGE = 0.01


def assert_geometry(gdf: gpd.GeoDataFrame, expected_types: set[str],
                    jurisdiction: gpd.GeoDataFrame | None = None,
                    id_col: str = "geoid") -> dict:
    """Run the layered checks and return counts per defect class."""
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError("geometry checks require a projected metric CRS")
    out = {"n": int(len(gdf))}

    out["null"] = int(gdf.geometry.isna().sum())
    out["empty"] = int(gdf.geometry.is_empty.sum())
    invalid_mask = ~gdf.geometry.is_valid
    out["invalid"] = int(invalid_mask.sum())
    if out["invalid"]:
        reasons = gdf.loc[invalid_mask, "geometry"].map(explain_validity)
        out["invalid_reasons"] = (reasons.str.extract(r"^([A-Za-z ]+)")[0]
                                  .value_counts().to_dict())
    types = set(gdf.geometry.geom_type.dropna().unique())
    out["unexpected_types"] = sorted(types - expected_types)

    if jurisdiction is not None:
        j = jurisdiction.union_all().buffer(1.0)     # one metre border tolerance
        out["outside_jurisdiction"] = int((~gdf.geometry.within(j)).sum())

    if types & {"Polygon", "MultiPolygon"}:
        # Orientation: shapely normalises most operations, but some writers do not.
        ext = gdf.geometry.map(lambda g: g if g.geom_type == "Polygon"
                               else (max(g.geoms, key=lambda p: p.area) if not g.is_empty else g))
        out["clockwise_exterior"] = int(sum(
            1 for g in ext if g is not None and not g.is_empty
            and g.geom_type == "Polygon" and not g.exterior.is_ccw))

    for k, v in out.items():
        if k not in ("n",) and v not in (0, [], {}):
            log.warning("geometry check %-22s %s", k, v)
    return out


def assert_tessellation(gdf: gpd.GeoDataFrame, id_col: str = "geoid") -> dict:
    """For layers meant to tile the plane: find overlaps and gaps above noise."""
    g = gdf[[id_col, "geometry"]].reset_index(drop=True)
    pairs = gpd.sjoin(g, g, predicate="intersects", how="inner")
    pairs = pairs[pairs[f"{id_col}_left"] < pairs[f"{id_col}_right"]]
    overlap_area = []
    for r in pairs.itertuples():
        a = g.geometry.iloc[r.Index]
        b = g.geometry.iloc[r.index_right]
        overlap_area.append(a.intersection(b).area)
    overlap_area = np.asarray(overlap_area)
    n_overlaps = int((overlap_area > SLIVER_M2).sum())

    union = g.geometry.union_all()
    hull = shapely.concave_hull(union, ratio=0.02) if not union.is_empty else union
    gaps = hull.difference(union)
    gap_parts = [p for p in getattr(gaps, "geoms", [gaps]) if p.area > SLIVER_M2]
    result = {"overlaps_above_noise": n_overlaps,
              "max_overlap_m2": float(overlap_area.max()) if overlap_area.size else 0.0,
              "gaps_above_noise": len(gap_parts)}
    log.info("tessellation: %s", result)
    return result


def repair(gdf: gpd.GeoDataFrame, id_col: str = "geoid") -> tuple[gpd.GeoDataFrame, pd.DataFrame]:
    """Repair invalid geometry and REPORT how much each repair changed it."""
    g = gdf.copy()
    bad = ~g.geometry.is_valid
    before = g.loc[bad, "geometry"].area
    g.loc[bad, "geometry"] = g.loc[bad, "geometry"].map(make_valid)
    # make_valid can return a collection; keep only the areal parts.
    g.loc[bad, "geometry"] = g.loc[bad, "geometry"].map(
        lambda x: shapely.union_all([p for p in getattr(x, "geoms", [x])
                                     if p.geom_type in ("Polygon", "MultiPolygon")]))
    after = g.loc[bad, "geometry"].area
    change = ((after - before).abs() / before.replace(0, np.nan)).fillna(0)
    report = pd.DataFrame({id_col: g.loc[bad, id_col], "area_change_frac": change})
    big = report[report["area_change_frac"] > MAX_REPAIR_AREA_CHANGE]
    log.info("repaired %d geometries; %d changed area by more than %.0f%%",
             int(bad.sum()), len(big), 100 * MAX_REPAIR_AREA_CHANGE)
    if len(big):
        log.error("repairs that materially changed geometry need review before use")
    return g, report

The repair report is the part that turns a convenience function into a controlled operation. A validity repair on a bow-tie polygon can drop half its area, and a pipeline that repairs and continues has published counts for an area that no longer exists in the shape anybody expects. Reporting the fractional area change per repaired feature, and refusing to proceed when any repair is material, keeps the modification visible.

The tessellation check applies only to layers intended to tile the study area — tracts, counties, service areas — and it catches the defects validity cannot: overlaps that double-count and gaps that break contiguity. The sliver threshold filters digitising noise so the check reports defects that matter rather than every floating-point disagreement along a shared edge.

Most Repairs Are Harmless; the Rest Must Stop the Run A distribution of repaired features by the fractional change in area caused by repair. The large majority change by a fraction of a per cent, typically from snapping a ring closed or removing a zero-width spike, and are accepted. A small tail changes by more than one per cent, and a few by tens of per cent, typically bow-tie polygons where the repair keeps only one lobe. Those in the tail are flagged for human review and the run does not proceed with them, because the repaired shape is no longer the area the counts refer to. A repair that changes area materially has changed the analysis 1% area change accepted, logged review before use bow-ties that lost a lobe <0.1% 1–10% >30% fractional area change caused by repair

Validation & Edge Cases

Checking case points but not boundaries. Point layers rarely fail validity; boundary layers regularly do. Run the full suite whenever a reference vintage is registered, not only on nightly case data.

Precision-induced invalidity. Reprojecting or rounding coordinates can make a valid polygon invalid by collapsing a narrow neck. Check validity after every transformation that changes coordinates, not only at ingestion.

Multipart features where single parts were expected. A tract that arrives as a MultiPolygon because of an offshore island is legitimate; one that arrives as two parts because of an editing error is not. The type check flags both, and the resolution is a human decision recorded against the identifier.

A jurisdiction polygon with its own defects. Containment against an invalid jurisdiction polygon produces nonsense. Validate the reference before using it to validate anything else.

Tolerance that hides a real problem. A sliver threshold large enough to suppress every warning also suppresses real overlaps. Derive it from the digitising precision of the source, and review the largest suppressed overlap periodically.

WARNING dq.geometry geometry check invalid                 3
WARNING dq.geometry geometry check invalid_reasons         {'Self-intersection': 2, 'Ring Self-intersection': 1}
INFO dq.geometry tessellation: {'overlaps_above_noise': 1, 'max_overlap_m2': 18422.6, 'gaps_above_noise': 0}
INFO dq.geometry repaired 3 geometries; 1 changed area by more than 1%
ERROR dq.geometry repairs that materially changed geometry need review before use

Precision-induced invalidity is the failure that appears after a step that should have been harmless, and a single narrow neck shows why:

Rounding Can Make a Valid Polygon Invalid A valid polygon shaped like an hourglass with a very narrow neck, where the two sides of the neck are a fraction of a metre apart. After coordinates are rounded to a coarser precision for storage, the vertices on either side of the neck snap to the same or crossing positions, and the ring now touches or crosses itself. The polygon was valid before the transformation and invalid after it, which is why validity must be rechecked after any step that changes coordinates. Check validity again after every step that moves coordinates Before rounding: valid neck 0.4 m wide After rounding: self-intersecting neck vertices crossed round to 1 m the transformation raised no error; only a re-check finds it

Compliance Notes

Record, per layer and per vintage: the counts from each check, the invalidity reasons, the tessellation results with the sliver threshold used, every repair with its fractional area change, and the review decision for any material repair.

Repaired geometry is a modified version of an authoritative source. When a repaired boundary is used to publish counts, the product’s metadata should state that the source geometry was repaired, by what method, and how many features were affected. An interagency recipient comparing your tract counts with their own will otherwise see small discrepancies at exactly the features you repaired and have no way to explain them.