Classifying Topology Errors by Severity

You have a spatial dataset — municipal parcel boundaries, utility networks, or environmental monitoring zones — and your validation run has returned hundreds of flagged features. The critical engineering question is not just what broke, but how badly: a self-intersecting polygon that causes a database join to crash is a completely different operational problem from a polygon with 8 000 vertices that merely degrades tile-server performance. This page shows you how to build a deterministic severity classifier so that your categorizing and prioritizing spatial errors workflow routes each topology violation to the right handler — pipeline halt, automated repair queue, or scheduled batch cleanup — rather than dumping everything into a single error log.

The approach uses Python 3.11+, GeoPandas 1.0+, and Shapely 2.0+. All tolerance thresholds are expressed in metric CRS units. Coordinate reference system (CRS) precision standards must be enforced before any topology check; applying degree-unit tolerances to unprojected WGS 84 data produces results that vary by up to a factor of 111 across latitudes.

Severity Tier Overview

Prerequisites

Before running the classifier, confirm:

• Python 3.11+ with geopandas>=1.0, shapely>=2.0, pandas>=2.0, and numpy>=1.26
• Metric CRS applied — project all input layers to UTM, State Plane, or an equivalent equal-area projection. Topology tolerances in this guide are in metres.
• Structured validation output — your rule engine (see building rule engines with GeoPandas for serialisation patterns) must produce a GeoDataFrame with a valid .geometry column; missing or null geometries must be handled upstream.
• OGC topology rules understood — if you need a reference on what constitutes a topologically valid geometry per the Open Geospatial Consortium (OGC) Simple Features specification, read understanding OGC topology rules first.
• Overlap detection strategy chosen — for datasets under ~50 000 features, the pairwise self-join below is sufficient. For larger datasets, replace it with a PostGIS ST_Relate query or an R-tree-indexed shapely.STRtree sweep.

Step-by-Step Procedure

1. Load and project your dataset

import geopandas as gpd

# Load source data — replace with your file or DB connection
gdf = gpd.read_file("input_boundaries.gpkg")

# Confirm geometry column is populated
assert gdf.geometry.notna().all(), "Null geometries must be filtered before classification"

# Project to a metric CRS — use the appropriate UTM zone for your region
if gdf.crs is None or gdf.crs.is_geographic:
    gdf = gdf.to_crs("EPSG:32633")   # UTM zone 33N — adjust to your area

Verification: print(gdf.crs) should return a projected CRS (units: metre, not degree).

2. Define configurable severity thresholds

Store thresholds in a dict or YAML config so they can be version-controlled and adjusted per dataset class:

THRESHOLDS = {
    "overlap_area_m2": 1.0,      # P1 overlap: area in sq metres above which a shared zone is flagged
    "gap_area_m2": 0.25,         # P1 gap: minimum enclosed gap area to report
    "max_vertices": 5000,        # P2 vertex count ceiling per polygon exterior
    "sliver_ratio": 50.0,        # P1 sliver: perimeter / sqrt(area) > this value
}

Anchoring thresholds to the ISO 19157 Geographic Information Data Quality standard’s logical consistency elements ensures your classification maps directly onto conformance metrics that compliance officers can audit.

3. Run the severity classifier

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

def classify_topology_errors(
    gdf: gpd.GeoDataFrame,
    thresholds: dict,
) -> pd.DataFrame:
    """
    Classify topology errors by severity (P0, P1, P2).

    Returns a DataFrame with columns:
        row_index, severity, violation_type, details, geometry_wkt
    """
    if gdf.empty:
        return pd.DataFrame(
            columns=["row_index", "severity", "violation_type", "details", "geometry_wkt"]
        )

    reports: list[dict] = []

    # ── P0: Invalid geometries (self-intersections, ring orientation, etc.) ──
    invalid_mask = ~gdf.geometry.is_valid
    for idx in gdf.index[invalid_mask]:
        geom = gdf.loc[idx, "geometry"]
        explanation = explain_validity(geom)
        try:
            fixed = make_valid(geom)
            detail = f"{explanation} → repaired as {fixed.geom_type}"
        except Exception as exc:
            detail = f"{explanation} → unrepairable: {exc}"
        reports.append({
            "row_index": idx,
            "severity": "P0",
            "violation_type": "Invalid Geometry",
            "details": detail,
            "geometry_wkt": geom.wkt[:200],
        })

    # ── P1: Overlapping areas exceeding threshold ──
    if len(gdf) > 1:
        tree = gdf.sindex
        for idx, geom in gdf.geometry.items():
            candidate_idxs = list(tree.intersection(geom.bounds))
            for cidx in candidate_idxs:
                if cidx <= idx:   # avoid duplicate pairs
                    continue
                other = gdf.geometry.iloc[cidx]
                if geom.intersects(other):
                    overlap = geom.intersection(other)
                    area = overlap.area
                    if area > thresholds["overlap_area_m2"]:
                        reports.append({
                            "row_index": idx,
                            "severity": "P1",
                            "violation_type": "Polygon Overlap",
                            "details": (
                                f"Overlaps feature at positional index {cidx}; "
                                f"shared area = {area:.2f} m²"
                            ),
                            "geometry_wkt": overlap.wkt[:200],
                        })

    # ── P1: Sliver polygons (high perimeter-to-area ratio) ──
    def sliver_ratio(geom) -> float:
        if geom is None or geom.is_empty or geom.area == 0:
            return 0.0
        return geom.length / (geom.area ** 0.5)

    ratios = gdf.geometry.apply(sliver_ratio)
    for idx in gdf.index[ratios > thresholds["sliver_ratio"]]:
        geom = gdf.loc[idx, "geometry"]
        reports.append({
            "row_index": idx,
            "severity": "P1",
            "violation_type": "Sliver Polygon",
            "details": f"Sliver ratio = {ratios.loc[idx]:.1f} (threshold: {thresholds['sliver_ratio']})",
            "geometry_wkt": geom.wkt[:200],
        })

    # ── P2: Excessive vertex density ──
    def vertex_count(geom) -> int:
        if geom is None or geom.is_empty:
            return 0
        if hasattr(geom, "exterior"):
            return len(geom.exterior.coords)
        if hasattr(geom, "geoms"):   # MultiPolygon
            return sum(len(p.exterior.coords) for p in geom.geoms)
        return 0

    counts = gdf.geometry.apply(vertex_count)
    for idx in gdf.index[counts > thresholds["max_vertices"]]:
        geom = gdf.loc[idx, "geometry"]
        reports.append({
            "row_index": idx,
            "severity": "P2",
            "violation_type": "High Vertex Density",
            "details": f"{counts.loc[idx]} vertices (threshold: {thresholds['max_vertices']})",
            "geometry_wkt": "",
        })

    return pd.DataFrame(reports)

Verification: After running, confirm the report has no rows with severity outside {"P0", "P1", "P2"}:

report = classify_topology_errors(gdf, THRESHOLDS)
assert report["severity"].isin({"P0", "P1", "P2"}).all()
print(report.groupby("severity").size())

4. Route errors by severity tier

p0 = report[report["severity"] == "P0"]
p1 = report[report["severity"] == "P1"]
p2 = report[report["severity"] == "P2"]

if not p0.empty:
    p0.to_csv("quarantine/p0_violations.csv", index=False)
    raise RuntimeError(
        f"Pipeline halted: {len(p0)} P0 topology violation(s) detected. "
        "Review quarantine/p0_violations.csv before reprocessing."
    )

if not p1.empty:
    p1.to_csv("repair_queue/p1_candidates.csv", index=False)
    print(f"[WARN] {len(p1)} P1 violation(s) queued for automated repair.")

if not p2.empty:
    p2.to_csv("logs/p2_performance.csv", index=False)
    print(f"[INFO] {len(p2)} P2 performance issue(s) logged for batch cleanup.")

5. Attempt automated P1 repair

For overlap and sliver violations, shapely.buffer(0) closes minor ring defects, while make_valid handles more complex breakage:

from shapely.validation import make_valid

def repair_p1(gdf: gpd.GeoDataFrame, p1_report: pd.DataFrame) -> gpd.GeoDataFrame:
    """Apply make_valid to P1-flagged rows; return patched GeoDataFrame."""
    repaired = gdf.copy()
    for idx in p1_report["row_index"].unique():
        if idx in repaired.index:
            original = repaired.loc[idx, "geometry"]
            fixed = make_valid(original)
            repaired.loc[idx, "geometry"] = fixed
    return repaired

Always re-run the classifier on the repaired output before writing to production — make_valid can convert a self-intersecting polygon to a GeometryCollection that downstream tools cannot process without further normalisation.

Interpreting Results

The report columns map to actionable next steps:

For P0 Invalid Geometry rows, the details field contains the Shapely explain_validity string (e.g., Self-intersection[402142.31 5789034.12]). Copy the coordinate pair and use ST_IsValidDetail in PostGIS or the QGIS Geometry Checker plugin to locate the precise defect on the map.

Gotchas & Edge Cases

1. make_valid produces type-shifted output. A self-intersecting polygon repaired with make_valid may become a MultiPolygon or even a GeometryCollection. If downstream tools expect a single Polygon, add a geometry-type assertion after repair and reject features that cannot be normalised back.

2. Degree-unit tolerances produce inconsistent results. If gdf.crs.is_geographic returns True when you run the classifier, all area and length calculations are in degrees squared or degrees, not square metres. A 1.0 sq degree overlap near the equator is roughly 12 350 km² — far above any practical threshold. Always project first.

3. The self-join overlap check is O(n²) without spatial indexing. The implementation above uses sindex to reduce candidate pairs, but on datasets with many spatially clustered features the inner loop can still be slow. For feature counts above 100 000, partition by a grid key or switch to a PostGIS batch query using ST_Intersects with a spatial index.

4. Slivers that span tile boundaries are invisible to single-tile checks. A sliver 0.1 m wide but 2 km long may not be detected if your input is split by tile and the sliver straddles the edge. Run topology checks on the full unclipped layer, or add a tile-boundary buffer before splitting.

5. explain_validity output differs across Shapely versions. Shapely 2.0 uses GEOS 3.12+ which reformatted some error strings. If you are parsing explain_validity output in downstream logic (rather than just logging it), pin shapely>=2.0,<3.0 and test your parser against the current string format.

When to Escalate

The GeoPandas-based approach above is suited for datasets up to approximately 500 000 features on a single node with 16 GB RAM. Escalate to a more powerful method when:

• Feature count exceeds 500 000 — switch to PostGIS ST_IsValid, ST_MakeValid, and ST_Intersects with a GiST index. Set-based SQL topology checks run in the database without serialising geometries to Python.
• Pairwise overlap detection runs out of memory — implement the overlap check using a spatial database or an asynchronous validation workflow with partitioned workers so each worker only checks a spatial tile.
• Your dataset spans multiple CRS zones — a global dataset whose features span multiple UTM zones cannot be checked with a single projected CRS. Use a custom equal-area projection (e.g., EPSG:6933 for global datasets) or split by zone and merge reports.
• Compliance reporting requires formal conformance statements — the classifier above produces operational severity flags. If you need a formal data quality report conforming to ISO 19157-1 (with DQ_TopologicalConsistency elements and conformance class URIs), the report schema must be extended or generated by a dedicated quality reporting framework.

Related

• Categorizing and Prioritizing Spatial Errors — parent cluster covering the full error triage workflow
• Understanding OGC Topology Rules — the geometric validity constraints this classifier enforces
• Building Rule Engines with GeoPandas — how to author the upstream rules that feed this classifier
• Asynchronous Validation Workflows — for running topology checks at scale with partitioned workers

Back to Categorizing and Prioritizing Spatial Errors