Spatial Diff Algorithms for Polygon Data

Compute topology-aware Boolean set operations that classify every polygon change between dataset revisions — part of the Branching & Merge Strategies for Spatial Datasets discipline.

Unlike tabular diffs that compare rows and cells, polygon diffing must reason about shared boundaries, coordinate precision drift, topological adjacency, and area conservation. A naive row-by-row equality check fails the moment a coordinate shifts by a single floating-point ULP; a geometry-unaware diff tool silently corrupts topology. The algorithms below treat two polygon layers as geometric sets, apply Boolean partitioning to classify change types, and attach structured metadata to every output feature so that downstream automated conflict detection in merge requests can gate commits on objective quality thresholds.

Prerequisites & Environment Setup

Before writing any diff logic, verify that your environment and input datasets satisfy strict geometric invariants. Violations cause silent errors or degenerate GEOS outputs that are harder to debug than upfront validation failures.

• Python 3.10+ with Shapely 2.x and GeoPandas 0.14+ — Shapely 2 ships with the GEOS C library bundled and exposes vectorised operations; the 1.x series is end-of-life and lacks set_precision.
• Consistent projected CRS on all layers — reproject to a metric CRS (e.g., EPSG:32633, EPSG:3857) before any operation. Geographic CRS (EPSG:4326) distorts area calculations and makes tolerance thresholds latitude-dependent.
• Geometry validity verified and repaired — run GeoDataFrame.is_valid on import; fix failures with shapely.make_valid before proceeding. Self-intersecting rings and duplicate consecutive vertices corrupt Boolean operations.
• Precision grid applied — call GeoDataFrame.set_precision(grid_size) with a project-wide standard (e.g., 1e-6 decimal degrees or 0.001 metres) on both layers before diffing. This is the single most effective sliver-prevention measure.
• Shared feature ID column — both baseline and target layers must carry an immutable identifier column (e.g., feature_id) so that geometric change can be separated from feature lifecycle events (add/remove).
• Shapely GEOS backend confirmed — run shapely.geos_version to confirm version ≥ 3.11; earlier releases have precision-model bugs that produce incorrect symmetric difference outputs.

Core Algorithmic Patterns

1. Boolean Set Partition (BSP) Diff

The foundational algorithm treats two polygon layers as mathematical sets A (baseline) and B (target) and computes four spatial regions:

BSP diff has O(n log n) spatial complexity when backed by an R-tree index (the default in GeoPandas overlay). Without spatial indexing, naive pairwise comparison is O(n²) and fails on any layer above ~50 k features.

2. Feature-Identity Diff

When each polygon carries a stable ID, geometric change detection is cheaper than overlay — compare geometries by ID using .equals() or symmetric-difference area before falling back to GEOS Boolean operations. This two-pass strategy reduces the number of expensive GEOS calls to only the truly modified subset.

The output for each feature is one of four change states: added, removed, modified, or unchanged. Attaching area_delta (signed area change in m²) and centroid_shift (Euclidean distance in metres) to each modified feature gives reviewers quantitative metrics for deciding whether a change is a legitimate edit or a coordinate-precision artefact.

3. Topological Overlay Diff

For datasets without a stable ID — common in OpenStreetMap extracts or third-party basemaps — a full topological overlay is required. geopandas.overlay(baseline, target, how="symmetric_difference") partitions the combined geometry into non-overlapping fragments, each labelled with source membership. Fragments that belong to only one source are classified as added or removed; fragments shared by both indicate topology-preserving edits such as boundary adjustments.

This approach is more expensive (full planar graph construction) but handles the resolution of overlapping polygons in collaborative editing correctly — it does not rely on IDs that may have been reassigned during merge.

Production Workflow Implementation

The diagram below shows the end-to-end diff pipeline, from raw inputs through validation to versioned output:

Step 1 — Load and repair inputs

import geopandas as gpd
from shapely.make_valid import make_valid   # Shapely 2.x

def load_and_repair(path: str, target_crs: str = "EPSG:32633") -> gpd.GeoDataFrame:
    gdf = gpd.read_file(path)
    gdf = gdf.to_crs(target_crs)
    invalid_mask = ~gdf.is_valid
    if invalid_mask.any():
        gdf.loc[invalid_mask, "geometry"] = gdf.loc[invalid_mask, "geometry"].apply(make_valid)
    return gdf

Step 2 — Apply precision snapping

PRECISION_GRID = 0.001  # 1 mm in metres for a metric CRS

def snap_precision(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    gdf = gdf.copy()
    gdf["geometry"] = gdf.geometry.set_precision(PRECISION_GRID)
    # Drop any zero-area artefacts produced by snapping
    gdf = gdf[gdf.geometry.area > PRECISION_GRID ** 2]
    return gdf

Step 3 — Feature-identity diff

import pandas as pd

def compute_polygon_diff(
    baseline: gpd.GeoDataFrame,
    target: gpd.GeoDataFrame,
    id_col: str = "feature_id",
) -> gpd.GeoDataFrame:
    """
    Returns a GeoDataFrame with change_type, area_delta, and centroid_shift
    columns for every feature in baseline and target.
    Both inputs must already share the same CRS and precision grid.
    """
    base_ids = set(baseline[id_col])
    tgt_ids = set(target[id_col])

    added_ids   = tgt_ids - base_ids
    removed_ids = base_ids - tgt_ids
    shared_ids  = base_ids & tgt_ids

    base_idx = baseline.set_index(id_col)
    tgt_idx  = target.set_index(id_col)

    modified_ids, unchanged_ids = set(), set()
    for fid in shared_ids:
        if not base_idx.loc[fid, "geometry"].equals(tgt_idx.loc[fid, "geometry"]):
            modified_ids.add(fid)
        else:
            unchanged_ids.add(fid)

    def _label(gdf: gpd.GeoDataFrame, ids, label: str) -> gpd.GeoDataFrame:
        out = gdf[gdf[id_col].isin(ids)].copy()
        out["change_type"] = label
        return out

    added     = _label(target,   added_ids,     "added")
    removed   = _label(baseline, removed_ids,   "removed")
    unchanged = _label(target,   unchanged_ids, "unchanged")

    # Annotate modified features with quantitative deltas
    mod_base = base_idx.loc[list(modified_ids)]
    mod_tgt  = tgt_idx.loc[list(modified_ids)].copy()
    mod_tgt["change_type"]    = "modified"
    mod_tgt["area_delta"]     = mod_tgt.geometry.area - mod_base.geometry.area
    mod_tgt["centroid_shift"] = mod_tgt.geometry.centroid.distance(
        mod_base.geometry.centroid
    )
    mod_tgt = mod_tgt.reset_index()

    result = pd.concat([added, removed, unchanged, mod_tgt], ignore_index=True)
    return gpd.GeoDataFrame(result, geometry="geometry", crs=baseline.crs)

For deeper implementation details — chunked processing, tolerance snapping edge cases, and attribute reconciliation — see Implementing spatial diff algorithms in Python.

Step 4 — Topology re-validation and serialisation

def validate_and_save(diff: gpd.GeoDataFrame, output_path: str, commit_hash: str) -> None:
    invalid = diff[~diff.is_valid]
    if not invalid.empty:
        raise ValueError(f"{len(invalid)} invalid geometries in diff output — fix before saving")

    diff = diff.copy()
    diff["commit_hash"] = commit_hash
    diff.to_parquet(output_path, engine="pyarrow", index=False)

Code Reliability Patterns

Defensive programming at the geometry layer prevents silent corruption:

• Wrap GEOS calls in try/except. shapely.errors.TopologicalError surfaces when two geometries are genuinely incompatible after snapping. Log the offending feature ID, skip the pair, and emit a warning — do not crash the pipeline.
• Re-check validity after every mutation. set_precision, buffer(0), and make_valid can each produce new degenerate geometries. Run is_valid immediately after each step.
• Normalise ring orientation before equality checks. Some sources use clockwise exterior rings; Shapely normalises counter-clockwise. Call shapely.normalize() on both geometries before .equals() comparisons to avoid false positives.
• Use equals_exact(other, tolerance) for near-identical geometries. Plain .equals() requires bit-identical coordinates; equals_exact with a small tolerance correctly identifies features that have drifted by less than the precision grid.
• Emit structured JSON logs. Each diff run should produce a summary {"total": n, "added": a, "removed": r, "modified": m, "unchanged": u, "commit": "..."} — this feeds the automated patching workflow for minor geometry shifts downstream.

Performance & Scale Considerations

Key tuning practices:

• Build R-tree indexes before overlay operations. GeoPandas overlay() uses STRtree internally; ensure you pass keep_geom_type=True and explicitly build GeoDataFrame.sindex before joining to avoid lazy index construction inside the loop.
• Partition by bounding envelope, not row count. Spatial partitions minimise cross-tile edge features that require duplicate processing; row-count partitions produce uneven geometric complexity per tile.
• Stream outputs to GeoParquet row groups. Write each tile’s diff results as a separate Parquet row group rather than accumulating all results in memory before writing. pyarrow.parquet.ParquetWriter supports incremental appends.
• Pin GEOS version in CI. GEOS patch releases occasionally change floating-point rounding in set operations, producing different (but equally valid) boundary geometries. Pinning prevents non-deterministic diffs across pipeline runs.

The delta tracking algorithms for vector data article covers the complementary problem of encoding and compressing the deltas once they have been computed.

Troubleshooting & Failure Modes

FAQ

Related

• Implementing Spatial Diff Algorithms in Python — step-by-step code walkthrough with chunking, attribute reconciliation, and CLI integration
• Automated Conflict Detection in Merge Requests — how diff outputs feed the merge-request validation gate
• Geometry Overlap Resolution Techniques — resolving the overlapping-polygon conflicts that diffs surface
• Delta Tracking Algorithms for Vector Data — encoding and compressing the deltas once computed
• Release Tagging Strategies for Spatial Basemaps — promoting a validated diff result into a versioned release

Back to Branching & Merge Strategies for Spatial Datasets