Topology Rule Enforcement

11 min read

A topology rule is a deterministic check that interrogates a relationship between two or more features — whether parcels overlap, whether utility lines terminate at a node, whether a building falls inside its declared zone — and returns a pass/fail verdict under an explicit tolerance and snapping model. It matters because relationship defects are invisible to single-feature checks: every geometry in a layer can be individually valid yet collectively wrong, with slivers between adjacent polygons, dangling network ends a few micrometres short of their junctions, or holes that should have been filled by a neighbour. This page sits directly beneath Spatial Test Pattern Design & Implementation, the implementation layer of the core geospatial QA architecture, and catalogues the relationship-level rules that run after geometry is structurally clean but before attribute and service checks. It is written for GIS QA engineers, data engineers, and platform teams who need topology checks that behave identically on a laptop and a memory-constrained CI runner.

Each branch answers a different relational question. Adjacency and coverage rules ask do these features tile space correctly — no two parcels overlapping, no gap between them. Connectivity rules ask is this network wired — every line endpoint coincident with a node within tolerance, no dangling stubs. Containment rules ask is this feature where its parent says it is — a point inside its polygon, a polygon inside its administrative boundary. Exclusion rules ask are these features kept apart — a building disjoint from a flood mask. A mature suite layers all four and, exactly as the GIS test pyramid prescribes, runs the cheapest index-prefilterable predicates before the expensive boundary-touching ones.

Rule Taxonomy and Tolerance Reference

The table below maps each topology rule to the spatial predicate that evaluates it, the tolerance strategy it requires, a typical threshold range, and the CRS unit that threshold is expressed in. A threshold is meaningless without its unit: a snap of 1e-6 is sub-micrometre in a projected CRS measured in metres but roughly 0.1 m in EPSG:4326 degrees at the equator. Read this table together with the rule for setting up spatial tolerance thresholds in assertions, which derives the numbers below from the CRS unit and the operation that produced the relationship.

Two columns deserve emphasis. First, almost every metric threshold is projected-only: signed area, snap radius, and setback distances are linear or areal quantities that only make sense in a CRS whose units are metres. Evaluate them in a geographic CRS and a 0.01 m snap silently becomes a 0.01° ≈ 1.1 km snap. Second, the boolean rules (within, intersects) carry no tolerance at all — adding an epsilon there masks an ETL defect rather than absorbing float noise.

DE-9IM: the Model Behind Every Predicate

Every predicate in the table is shorthand for a constraint on the Dimensionally Extended 9-Intersection Model (DE-9IM) matrix — the 3×3 matrix recording the dimension of the intersection between the interior, boundary, and exterior of two geometries. Knowing the matrix is what lets you distinguish a true overlap from a features-merely-touching case at a shared edge. The intersection matrix is

where each cell is -1 (empty, F), 0, 1, or 2. Named predicates are pattern matches over this matrix, so a precise rule can target a pattern Shapely’s helpers do not expose directly:

from shapely import relate, relate_pattern

# Two parcels share only an edge — valid adjacency, NOT an overlap.
relate(parcel_a, parcel_b)            # e.g. 'FF2F11212'
relate_pattern(parcel_a, parcel_b, "F***1****")  # interiors disjoint, boundaries meet on a line

The F***1**** mask above is exactly the “shares a boundary line but interiors do not intersect” rule that a must-not-overlap coverage layer needs — a constraint overlaps() alone cannot express, because two polygons sharing only an edge are not “overlapping” yet must still be distinguished from disjoint neighbours.

Adjacency and Gap-Free Coverage

A coverage layer (parcels, census blocks, land-use polygons) must tile its extent with no overlaps and no gaps. The overlap half is a pairwise predicate; the gap half is an areal identity — the union of the parts must reconstruct the whole within an area delta $\varepsilon_A$:

from shapely import overlaps, union_all

assert not any(overlaps(a, b) for a, b in candidate_pairs)   # no double coverage
gap = extent.area - union_all(parts).area                    # residual uncovered area
assert abs(gap) <= 1e-3                                       # m² floor, projected CRS

The area floor exists because union_all accumulates floating-point error proportional to vertex count; a non-zero $\varepsilon_A$ absorbs that noise without hiding a real hole. Compute it in a projected CRS, never in degrees.

Connectivity and Dangle Detection

A network is connected when every line endpoint coincides with a node — another endpoint or a permitted junction — inside a snap tolerance $\tau$. Two endpoints are the “same” node when $\lVert p_i - p_j \rVert \le \tau$. Endpoints with no match within $\tau$ are dangles, the classic source of routing failures:

from shapely import touches, snap

snapped = snap(line_a, line_b, tolerance=0.01)   # metres, projected CRS
assert touches(snapped, line_b)                  # endpoints now coincide at a node

Choosing $\tau$ is a trade-off: too small and legitimate junctions read as dangles; too large and two distinct nodes collapse into one, silently rewiring the graph. Derive $\tau$ from the capture tolerance of the source survey, not from a round number.

Containment and Exclusion

Containment asserts a feature lies inside a parent (within / covered_by); exclusion asserts two layers stay apart (disjoint, optionally with a buffered setback). These are largely boolean and tolerance-free, but covered_by is the correct choice over within whenever a child may legitimately touch its parent’s boundary — within is False for a feature flush against the edge, which produces false failures on coastlines and parcel fits.

from shapely import covered_by, disjoint

assert covered_by(building, zone)            # building inside its zone (boundary-inclusive)
assert disjoint(building, flood_mask.buffer(0))  # buffer(0) cleans the mask first

Pre-cleaning matters: these predicates assume valid inputs, which is why topology checks run after the geometry validation patterns have already repaired ring orientation, removed self-intersections, and snapped precision — a self-intersecting parent geometry makes covered_by undefined.

Production-Grade Python Implementation

Production topology validation must operate on a declarative rule schema rather than imperative spatial queries scattered through test code. Hardcoded predicates drift across coordinate reference systems and rot under maintenance. Drive a rule engine from config, resolving tolerance to the native linear unit of the target CRS before any predicate runs:

# topology_rules.yaml
topology_rules:
  - id: rule_parcel_no_overlap
    description: "Parcel boundaries must not overlap"
    predicate: must_not_overlap
    tolerance_m: 0.01          # metres in the working projected CRS
    snap_mode: shared_boundary
    severity: critical
    action: fail_pipeline
  - id: rule_network_connectivity
    description: "Utility lines must terminate at network nodes"
    predicate: endpoint_touches
    tolerance_m: 0.005
    snap_mode: endpoint_only
    severity: warning
    action: quarantine

# test_topology_rules.py — pytest 7+, Shapely 2.x, GeoPandas 0.14+
import yaml
import geopandas as gpd
from shapely import STRtree, overlaps, snap, touches

PREDICATES = {
    "must_not_overlap": lambda a, b, tol: not overlaps(snap(a, b, tol), b),
    "endpoint_touches": lambda a, b, tol: touches(snap(a, b, tol), b),
}

def load_rules(path="topology_rules.yaml"):
    return yaml.safe_load(open(path))["topology_rules"]

def to_metric_tolerance(gdf, tol_m):
    # Reject ambiguous degree-based tolerances up front.
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError("topology tolerances require a projected CRS (metres)")
    return tol_m  # already in the layer's linear unit

def test_topology(rule, gdf: gpd.GeoDataFrame):
    tol = to_metric_tolerance(gdf, rule["tolerance_m"])
    check = PREDICATES[rule["predicate"]]
    tree = STRtree(gdf.geometry.values)          # R-tree prefilter, Shapely 2.x
    violations = []
    for i, geom in enumerate(gdf.geometry.values):
        for j in tree.query(geom):               # only bbox-overlapping candidates
            if j <= i:
                continue                          # deterministic, deduplicated pairs
            if not check(geom, gdf.geometry.values[j], tol):
                violations.append((rule["id"], i, j))
    assert not violations, f"{rule['id']}: {len(violations)} violations e.g. {violations[:3]}"

Three details make this deterministic rather than flaky. The STRtree index turns an $O(n^2)$ pairwise scan into a candidate prefilter, so only bounding-box-overlapping features reach the expensive predicate. The j <= i guard yields each unordered pair exactly once and in a stable order, guaranteeing identical output across parallel runners. And to_metric_tolerance hard-rejects a geographic CRS instead of silently evaluating a metre threshold in degrees.

PostGIS Database-Side Counterparts

For data already in PostGIS, push the same rules server-side so the index and the predicate live next to the data. The && bounding-box operator engages the GiST index before the exact predicate runs, mirroring the STRtree prefilter above:

-- Must-not-overlap: any pair of parcels whose interiors intersect
SELECT a.id, b.id
FROM parcels a
JOIN parcels b
  ON a.id < b.id              -- deterministic, deduplicated pairs
 AND a.geom && b.geom         -- GiST bbox prefilter
WHERE ST_Overlaps(a.geom, b.geom);

-- Dangle nodes: line endpoints not snapped to any other endpoint within 0.01 m
SELECT l.id
FROM lines l
WHERE NOT EXISTS (
  SELECT 1 FROM lines o
  WHERE o.id <> l.id
    AND ST_DWithin(ST_EndPoint(l.geom), o.geom, 0.01)  -- metres in a projected SRID
);

ST_ReducePrecision (PostGIS 3.1+) is the server-side equivalent of Shapely’s set_precision; pin the PostGIS GEOS build to match your Shapely wheel so the two engines snap identically, or the same dataset will pass in the database and fail in Python.

Memory-Safe Execution

Large-scale topology validation triggers out-of-memory failures when a naïve spatial join materialises a full Cartesian product. Bound the work explicitly:

1. Spatial partitioning. Split inputs into grid-aligned or quadtree tiles and process only intersecting tile pairs, so complexity scales with local feature density rather than total count.
2. Index-aware filtering. Always prefilter with an R-tree (STRtree) or GiST index before invoking overlaps/touches; the predicate is orders of magnitude costlier than a bounding-box test.
3. Streaming aggregation. Yield violations as a stream and write them to a Parquet or GeoJSON error layer incrementally rather than holding a full result set in RAM.
4. Deterministic ordering. Sort candidate pairs by tile ID and the i < j index so output is reproducible across distributed runners — a prerequisite for stable CI diffs.

For datasets exceeding millions of features, decouple the predicate evaluation from blocking I/O by running async spatial tests with pytest-asyncio, which parallelises tile processing across a worker pool without loading every tile at once.

Pipeline Integration and Observability

Topology rules earn their keep only when they run automatically, deterministically, and observably. Wire them into CI as severity-routed gates, and serialise every verdict as structured JSON so platform teams can track relationship-defect rates over time:

• Pre-merge: must-not-overlap and connectivity on the PR’s changed features — fast, fail-fast, blocks the merge.
• Nightly / full sync: gap-free coverage and containment across production-scale data on distributed runners.
• Artifact generation: emit a machine-readable verdict log plus a topology_violations.geojson error layer for QA triage.

Pin libgeos and PROJ in the container image so the relationship engine is bitwise reproducible; an unpinned PROJ grid shift moves coordinates by metres between builds and surfaces as a flaky connectivity failure. Because the same relationships must hold regardless of storage backend, pair these gates with cross-format parity testing so a rule that passes on GeoParquet also passes on FlatGeobuf, Shapefile, and PostGIS, and with attribute and metadata checks so a topologically valid network is also referentially intact — no orphaned segments or mismatched CRS metadata.

Common Failure Modes and Gotchas

1. Confusing touch with overlap. Two polygons sharing only an edge are not overlapping, but a careless intersects check flags them. Use the DE-9IM mask F***1**** or overlaps() to target true interior intersection.
2. Snap tolerance too large collapses nodes. A snap radius bigger than the smallest real gap between distinct junctions merges them, silently rewiring the network graph. Derive $\tau$ from the survey’s capture tolerance.
3. CRS unit mismatch on every metric threshold. Snap radii, setbacks, and area floors expressed in metres but evaluated in EPSG:4326 are off by a factor of ~100,000. Reproject to a projected CRS before any distance or area predicate.
4. within versus covered_by boundary case. A child feature flush against its parent’s boundary is within == False but covered_by == True. Use covered_by for containment unless you genuinely require strict interior containment.
5. Gap floor applied in degrees. A coverage area delta of 1e-3 m² evaluated against a geographic CRS either tolerates kilometre holes or rejects everything. Compute union and extent area in a projected CRS.
6. Running topology before geometry repair. overlaps, covered_by, and touches are undefined on invalid geometry; a self-intersecting parcel produces a meaningless verdict. Gate topology behind structural validity.
7. Non-deterministic pair ordering. Iterating pairs without an i < j (or a.id < b.id) guard yields duplicates and unstable output, making CI diffs noisy and masking regressions. Enforce a stable, deduplicated pair order.
8. Unpinned GEOS between Shapely and PostGIS. Divergent engine builds disagree at the snap boundary and produce “passes locally, fails in CI” reports — pin both and validate cross-engine.

Conclusion

Topology rule enforcement gives an engineer a decision procedure: identify which relationship a rule actually constrains — adjacency, connectivity, containment, or exclusion — express it as a DE-9IM pattern or named predicate, then apply that rule’s tolerance strategy in the correct projected CRS unit. Driven from a declarative schema, prefiltered through a spatial index, executed in memory-safe tiles, and routed by severity into CI gates with serialised verdicts, topology checks become deterministic gates rather than flaky scripts, and relationship defects are rejected at the earliest boundary instead of discovered in a broken routing graph downstream. For how these checks compose with geometry, attribute, and parity patterns into a full validation architecture, return to Spatial Test Pattern Design & Implementation.

Related Pages

• Geometry Validation Patterns
• Attribute & Metadata Checks
• Cross-Format Parity Testing
• Async Execution for Large Datasets
• Running async spatial tests with pytest-asyncio