Geospatial Schema Standardization & Mapping Automation

Computing Geometry Hashes for Change Detection

A content hash is only useful if it changes exactly when the content does. For attributes that is nearly automatic; for geometry it is not, because the same polygon has many valid byte encodings. Start the ring at a different vertex, wind it the other way, order the parts of a multipolygon differently, or reproject and reproject back, and the shape is identical while the bytes are not. A hash over raw WKB therefore reports change on every re-export, which makes an incremental pipeline slower than the full reload it replaced. This procedure builds a canonical form that hashes stably, implementing the detection half of incremental change detection and delta loads inside Automated Attribute Transformation & ETL Workflows.

The steps map to configure (Steps 1–3, the canonical form), execute (Step 4, the digest), and validate (Step 5, the sensitivity pair). The whole procedure is about twenty lines of code and about four decisions, and the decisions are the part worth writing down.

Prerequisites checklist

Step 1: Snap to a precision grid

python
# etl/geomhash.py — shapely >=2.0 — Python 3.10+
from shapely import set_precision
from shapely.geometry.base import BaseGeometry

GRID = 0.001        # 1 mm in a metre-based CRS


def quantize(geom: BaseGeometry, grid_size: float = GRID) -> BaseGeometry:
    """Snap ordinates to a fixed grid so arithmetic noise cannot alter the digest."""
    return set_precision(geom, grid_size)

Reprojection, buffering and even reading through different drivers introduce differences in the last few bits of a coordinate. Those differences are not data — they are arithmetic history — and without quantization every one of them changes the hash. Choose the grid from the data’s actual precision: 1 mm for cadastral and utility data, 1 cm for topographic layers, and finer only if you can defend it. A grid coarser than the smallest real edit will hide changes, so the number is a trade-off, not a default.

For geographic coordinates the equivalent grid is about 1e-8 degrees for millimetre-scale work, which is another reason to hash in a projected CRS wherever possible.

Step 2: Normalize ring order and the start vertex

python
# etl/geomhash.py — shapely >=2.0 — Python 3.10+
def canonical(geom: BaseGeometry, grid_size: float = GRID) -> bytes:
    """A byte form that is identical for geometrically identical shapes."""
    return quantize(geom, grid_size).normalize().wkb

normalize() is the single operation that makes geometry hashing practical. It imposes shapely’s canonical ordering: exterior rings wound one way and interiors the other, each ring rotated to a canonical start vertex, and the parts of a multi-geometry sorted. Two encodings of the same polygon — one from a shapefile, one from a GeoPackage, one that went through a desktop edit and came back — collapse to the same bytes.

Without it, the most common false change in the fleet is a shapefile round-trip, which reverses ring orientation as described in enforcing ring orientation and winding order — a change to every byte of every polygon and to none of the shapes.

Three Encodings, One Canonical Form Three squares are drawn side by side, each annotated with a different WKB encoding of the same shape: the first starts its ring at the lower-left corner, the second starts at the upper-right, and the third is a multipolygon whose two parts are listed in the opposite order. Arrows from all three lead down into a single canonical form box, which applies precision snapping and shapely normalize and yields one digest. A caption notes that without normalization these three would produce three different hashes and three false change events. starts lower-left starts upper-right, reversed parts listed in reverse canonical form set_precision(1 mm) normalize() then wkb one digest no change event without normalization these are three different digests and three false "updated" events

Step 3: Settle the Z and M policy explicitly

python
# etl/geomhash.py — shapely >=2.0 — Python 3.10+
from shapely import force_2d


def canonical(geom: BaseGeometry, grid_size: float = GRID, include_z: bool = False) -> bytes:
    prepared = geom if include_z else force_2d(geom)
    return set_precision(prepared, grid_size).normalize().wkb

Whether height participates in identity is a data question, not a technical one. For a parcel fabric where Z is an artefact of the capture process and is not maintained, including it means every re-capture registers as a boundary change. For a utility network where invert levels are the point of the dataset, excluding it means a re-levelled manhole registers as unchanged. Decide per layer, record the decision in the delta manifest, and note that set_precision does not quantize Z — if Z is included, quantize it separately or accept that its noise reaches the digest.

M values almost never belong in identity: they are usually linear-referencing measures recomputed on every export.

Step 4: Digest the canonical bytes together with the SRID

python
# etl/geomhash.py — Python 3.10+
import hashlib


def geometry_hash(geom: BaseGeometry, srid: int, grid_size: float = GRID,
                  include_z: bool = False) -> str:
    digest = hashlib.sha256()
    digest.update(srid.to_bytes(4, "big"))          # a reprojection IS a change
    digest.update(b"\x1e")                          # separator between srid and geometry
    digest.update(canonical(geom, grid_size, include_z))
    return digest.hexdigest()

Including the SRID makes a reprojection visible. Without it, a layer that starts arriving in a different CRS but describing identical ground hashes identically, and the delta reports no change — while every consumer sees coordinates move. Since the normalization stage should have brought everything to one CRS before this point, an SRID change here is a signal that something upstream broke, and the hash is a good place to notice it.

SHA-256 is the right default. It is fast enough that hashing is never the bottleneck — hundreds of megabytes per second per core — and collision-free for this purpose in any practical sense. Truncating the digest to 16 hex characters for storage is acceptable at layer sizes below a few hundred million features, but store the full digest if there is any chance of the data being used as evidence.

Step 5: Prove sensitivity and insensitivity together

python
# tests/test_geomhash.py — pytest >=7, shapely >=2.0
from shapely import affinity
from shapely.geometry import MultiPolygon, Polygon

from etl.geomhash import geometry_hash

SQUARE = Polygon([(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)])


def test_rotated_start_vertex_hashes_identically():
    rotated_coords = Polygon([(10, 10), (0, 10), (0, 0), (10, 0), (10, 10)])
    assert geometry_hash(SQUARE, 25832) == geometry_hash(rotated_coords, 25832)


def test_reversed_winding_hashes_identically():
    reversed_ring = Polygon(list(SQUARE.exterior.coords)[::-1])
    assert geometry_hash(SQUARE, 25832) == geometry_hash(reversed_ring, 25832)


def test_reordered_multipolygon_parts_hash_identically():
    a, b = SQUARE, affinity.translate(SQUARE, xoff=20)
    assert geometry_hash(MultiPolygon([a, b]), 25832) == geometry_hash(MultiPolygon([b, a]), 25832)


def test_one_centimetre_move_changes_the_hash():
    moved = affinity.translate(SQUARE, xoff=0.01)
    assert geometry_hash(SQUARE, 25832) != geometry_hash(moved, 25832)


def test_sub_grid_jitter_does_not_change_the_hash():
    jittered = affinity.translate(SQUARE, xoff=0.0001)      # 0.1 mm, below the 1 mm grid
    assert geometry_hash(SQUARE, 25832) == geometry_hash(jittered, 25832)


def test_srid_participates_in_identity():
    # Negative control: the same coordinates in a different CRS are NOT the same feature.
    assert geometry_hash(SQUARE, 25832) != geometry_hash(SQUARE, 25833)

Read the last three as a set. A hash tuned only for insensitivity passes the first three tests by ignoring everything; the 1 cm test is what forces it to remain a change detector. The jitter test pins the boundary between the two, and it is the one to revisit when the grid size changes.

Verification

Across a real layer, the diagnostic is the distribution rather than any single hash:

python
import geopandas as gpd

gdf = gpd.read_file("input/parcels.gpkg")
gdf["hash"] = [geometry_hash(g, 25832) for g in gdf.geometry]

print(gdf["hash"].nunique(), "distinct of", len(gdf))
# 214,338 distinct of 214,338     → no accidental collisions

Re-running the same hash over a re-export of the identical layer must produce a set-equal result:

text
INFO etl.delta hashes: 214,338 unchanged, 0 updated, 0 inserted, 0 deleted

Anything else after a pure re-export means the canonical form is still leaking a serialization detail.

Troubleshooting

Symptom Likely cause Fix
Every feature changes after a shapefile round-trip normalize() not applied, so reversed winding reaches the digest Add .normalize() to the canonical form
Every feature changes after a reprojection round-trip Grid too fine for the arithmetic noise the reprojection introduces Coarsen the grid to the data’s real precision, typically 1 mm
Two genuinely different parcels share a hash Both quantized to the same shape by a grid far coarser than the difference The grid is too coarse; it must be finer than the smallest meaningful edit
Hash changes when nothing did, only for 3D layers Z values carry noise and set_precision does not quantize them Set include_z=False, or round Z explicitly before hashing
Hashing is a measurable share of run time Hashing row by row in Python over tens of millions of features Hash in chunks and parallelize by partition; the digest itself is not the cost
Hash matches but the feature clearly moved The move is below the precision grid Confirm the grid against the accuracy of the source; a 1 m grid hides a 90 cm correction