Geospatial Schema Standardization & Mapping Automation

Setting Tolerance Thresholds for Automated Coordinate Snapping

Automated coordinate snapping fails predictably when the tolerance is hardcoded, ignored across a projection change, or decoupled from the source dataset’s declared precision. A single literal such as tolerance=0.001 means one millimetre in a metric projection but roughly 111 metres in unprojected WGS84 — apply it to the wrong dataset and you collapse valid micro-features, weld unrelated vertices together, and trigger silent topology drift that only surfaces when a downstream spatial join returns nonsense. This page gives the exact steps to derive a snapping tolerance from the active coordinate system, apply it deterministically with shapely.set_precision, and gate the result in CI.

This task is one operation inside Unit Conversion & Tolerance Thresholds, the tolerance-gating stage of Coordinate Reference System (CRS) Normalization & Sync. Snapping must run after projection so that the tolerance is expressed in the target CRS’s native unit; pre-transformation thresholds are mathematically invalid the moment a reproject changes the axis unit. All routines assume shapely >= 2.0, geopandas >= 0.14, and pyproj >= 3.6.

CRS-Aware Coordinate Snapping Decision Flow A GeoDataFrame enters the stage; resolve_tolerance reads the active CRS axis unit and branches three ways — geographic selects a 1e-6 degree grid, metre selects 0.001 m, and US survey foot selects 0.00328 ft. The chosen grid size feeds shapely.set_precision, which snaps every vertex. A topology-and-bounds-drift gate then asks whether geometry is valid and the extent moved no more than two grid cells: passing geometry is committed and written to the audit log, while failing geometry is sent to make_valid or rejected and the tolerance is tightened. Input GeoDataFrame gdf.crs is not None Resolve CRS Unit crs.axis_info[0].unit_name axis unit? geographic / m / ft degree → 1e-6 ~0.11 m at equator metre → 0.001 survey-grade band US foot → 0.00328 legacy State Plane shapely.set_precision grid_size = tolerance grid size valid & drift OK? drift ≤ 2 cells Commit + Audit Log epsg, tolerance, drift pass make_valid / Reject tighten tolerance fail

Prerequisites

Confirm the environment before running any snapping routine:

Step 1 — Configure Baseline Tolerances Per CRS Unit

Tolerance must scale with the coordinate system’s native unit and the dataset’s declared precision, so it belongs in configuration rather than buried in code. Externalising it keeps the value auditable and lets the same routine run across mixed-CRS inputs without edits. Use these defaults as the starting point and tighten them per data class:

python
# snap_config.py
# Requires: Python >= 3.10
# Snapping tolerance keyed by CRS axis unit, expressed in that unit.
SNAP_TOLERANCE: dict[str, float] = {
    # Decimal degrees (EPSG:4326 and other geographic CRSs).
    # 1e-6 deg ~= 0.11 m at the equator; 5e-6 ~= 0.55 m.
    "degree": 1e-6,       # cadastral / administrative boundaries
    # Metric projections (UTM, State Plane metric).
    # 0.001 m preserves survey-grade precision; 0.01 m suits GPS field data.
    "metre": 0.001,
    # US Survey Foot projections (legacy State Plane).
    # Direct metre-to-foot conversion of the metric band.
    "foot": 0.00328,
}

# Bands you may widen to per dataset class, in the unit above.
SNAP_TOLERANCE_RELAXED: dict[str, float] = {
    "degree": 5e-6,       # environmental / hydrographic linework
    "metre": 0.01,        # GPS-collected field data
    "foot": 0.0328,
}

# Reject if the post-snap bounding box shifts by more than this many grid cells.
MAX_BOUNDS_DRIFT_CELLS: float = 2.0

Keep the geographic band tight: 1e-6 for boundary data that must not move, 5e-6 only for environmental linework where sub-metre wobble is acceptable. For State Plane datasets encoded in US Survey Feet, note that the unit string reported by PROJ is distinct from the International Foot — Step 2 matches on the substring "foot" precisely to avoid silently mixing the two. The conversion factor that produced 0.00328 comes from the same canonical table documented in Unit Conversion & Tolerance Thresholds.

Step 2 — Derive the Tolerance and Apply Grid Snapping

Resolve the active CRS, read its axis unit, look up the matching tolerance, and snap. shapely.set_precision aligns every coordinate to a regular grid of the given size, which both removes floating-point drift and welds coincident vertices at shared boundaries — exactly what manual rounding cannot do safely. The routine copies the frame first so non-geometry attributes are never mutated.

python
# snap.py
# Requires: Python >= 3.10, shapely >= 2.0, geopandas >= 0.14, pyproj >= 3.6
import logging

import geopandas as gpd
import shapely
from pyproj import CRS

from snap_config import SNAP_TOLERANCE

logger = logging.getLogger(__name__)


def resolve_tolerance(gdf: gpd.GeoDataFrame) -> float:
    """Return the snapping tolerance for the GeoDataFrame's active CRS unit."""
    if gdf.crs is None:
        # Never assume a CRS — an unset CRS is a pipeline defect, not a default.
        raise ValueError(
            "GeoDataFrame lacks a defined CRS. Assign via .set_crs() before snapping."
        )

    crs = CRS.from_user_input(gdf.crs)
    unit_name = crs.axis_info[0].unit_name.lower() if crs.axis_info else ""

    if crs.is_geographic:
        return SNAP_TOLERANCE["degree"]
    if "metre" in unit_name or "meter" in unit_name:
        return SNAP_TOLERANCE["metre"]
    if "foot" in unit_name:
        return SNAP_TOLERANCE["foot"]

    logger.warning(
        "Unrecognised CRS unit '%s' (EPSG:%s); falling back to metric tolerance.",
        unit_name,
        crs.to_epsg(),
    )
    return SNAP_TOLERANCE["metre"]


def snap_with_dynamic_tolerance(gdf: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, float]:
    """Apply CRS-aware grid snapping; return the snapped frame and the tolerance used."""
    if gdf.empty:
        raise ValueError("Input GeoDataFrame is empty. Aborting snapping routine.")

    tolerance = resolve_tolerance(gdf)

    # set_precision snaps coordinates to a grid of `grid_size`, eliminating
    # floating-point drift and aligning coincident vertices deterministically.
    snapped = shapely.set_precision(gdf.geometry.values, grid_size=tolerance)

    out = gdf.copy()
    out["geometry"] = gpd.GeoSeries(snapped, index=gdf.index, crs=gdf.crs)
    return out, tolerance

Because the tolerance is keyed to the CRS unit, this routine slots directly behind a reproject. When you chain it after the Projection Normalization Workflows stage, the frame arriving here already carries the target CRS, so resolve_tolerance reads the post-reproject unit and never the source one. For datasets assembled from several source projections, run Multi-CRS Dataset Harmonization first so all features share one projected CRS — geographic-to-geographic snapping produces inconsistent ground distances across latitudes and must be avoided.

Step 3 — Validate Topology and Bounds Drift

Snapping can fold a thin sliver into a self-intersection, so the result is unsafe to commit until it has been checked. Two assertions catch the common failure modes: a validity check, and a bounding-box drift check that confirms the snap moved vertices only within its own grid cell rather than relocating features.

Correct vs Too-Coarse Snapping Tolerance on the Precision Grid Two panels share the same noisy input vertices of a quadrilateral parcel. Left panel, tolerance correct: each raw vertex (hollow dot) snaps a short distance to the nearest corner of a fine grid (filled dot); the original and snapped bounding boxes coincide, drift stays within two cells, and the gate passes. Right panel, tolerance too coarse: a coarse grid pulls several raw vertices onto the same cell corner, the snapped bounding box (solid) shrinks well inside the original bounding box (dashed) by more than two cells, the geometry degenerates, and the drift gate rejects the snap. Tolerance correct fine grid · drift ≤ 2 cells · PASS bbox unchanged Tolerance too coarse coarse grid · drift > 2 cells · REJECT original bbox snapped bbox shrinks raw vertex snapped vertex snap move to nearest cell corner
python
# validate.py
# Requires: Python >= 3.10, shapely >= 2.0, geopandas >= 0.14
import logging

import geopandas as gpd
import shapely

from snap_config import MAX_BOUNDS_DRIFT_CELLS

logger = logging.getLogger(__name__)


def validate_snap(
    original: gpd.GeoDataFrame,
    snapped: gpd.GeoDataFrame,
    tolerance: float,
    *,
    fail_on_invalid: bool = True,
) -> gpd.GeoDataFrame:
    """Assert topology validity and bounds stability; repair or raise on failure."""
    # 1. Topology validity.
    invalid_mask = ~snapped.geometry.is_valid
    invalid_count = int(invalid_mask.sum())
    if invalid_count:
        if fail_on_invalid:
            raise RuntimeError(
                f"Snapping produced {invalid_count} invalid geometries. "
                "Review source topology before relaxing fail_on_invalid."
            )
        logger.warning("Repairing %d invalid geometries post-snap.", invalid_count)
        snapped = snapped.copy()
        snapped.loc[invalid_mask, "geometry"] = shapely.make_valid(
            snapped.geometry[invalid_mask].values
        )

    # 2. Bounds drift — the extent must not move beyond N grid cells.
    o_bounds = original.total_bounds        # [minx, miny, maxx, maxy]
    s_bounds = snapped.total_bounds
    drift = abs(s_bounds - o_bounds).max()
    if drift > tolerance * MAX_BOUNDS_DRIFT_CELLS:
        raise RuntimeError(
            f"Post-snap bounds drifted {drift:g} (> {tolerance * MAX_BOUNDS_DRIFT_CELLS:g}). "
            "Tolerance is too coarse for this dataset's precision."
        )

    # 3. Stale spatial indexes return false-negative joins — force a rebuild.
    snapped.sindex  # noqa: B018  (accessing the property regenerates it)
    return snapped, invalid_count, float(drift)

The drift gate is what distinguishes a correct tolerance from one that is too coarse. If the extent shifts by more than two grid cells, the chosen radius is collapsing real geometry rather than just cleaning floating-point noise — the fix is to tighten the tolerance, not to widen the gate.

Step 4 — Log the Snapping Operation for Audit

Every transformation must record what it did so the lineage trail can reproduce it. Emit the resolved CRS, the tolerance applied, the invalid-geometry count, and the bounds delta as a single structured record. These fields satisfy ISO 19115 lineage expectations and let an auditor confirm the snap was deterministic. When this stage runs inside a larger run, the same record is the unit of lineage that Batch Schema Processing Pipelines aggregate per layer, so keep the field names stable across stages.

python
# pipeline.py
# Requires: Python >= 3.10, geopandas >= 0.14
import json
import logging

import geopandas as gpd

from snap import snap_with_dynamic_tolerance
from validate import validate_snap

logger = logging.getLogger("geo.snap.audit")


def snap_stage(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Run the full configure -> execute -> validate -> log snapping stage."""
    snapped, tolerance = snap_with_dynamic_tolerance(gdf)
    snapped, invalid_count, drift = validate_snap(gdf, snapped, tolerance)

    logger.info(
        json.dumps(
            {
                "stage": "coordinate_snapping",
                "epsg": gdf.crs.to_epsg(),
                "tolerance": tolerance,
                "feature_count": len(snapped),
                "repaired_geometries": invalid_count,
                "bounds_drift": drift,
            }
        )
    )
    return snapped

Verification

Confirm the stage behaves deterministically before merging. The test snaps a frame whose vertices carry sub-tolerance noise, then asserts validity, attribute preservation, and a sane bounds delta:

python
# test_snap.py
# Requires: Python >= 3.10, shapely >= 2.0, geopandas >= 0.14, pytest >= 7
import geopandas as gpd
from shapely.geometry import Polygon

from snap import snap_with_dynamic_tolerance
from validate import validate_snap


def test_metric_snap_is_valid_and_preserves_attributes():
    # A square with 1e-7 m noise on one vertex — well below the 0.001 m grid.
    poly = Polygon([(0, 0), (1, 0.0000001), (1, 1), (0, 1)])
    gdf = gpd.GeoDataFrame({"parcel_id": ["A1"]}, geometry=[poly], crs="EPSG:32633")

    snapped, tol = snap_with_dynamic_tolerance(gdf)
    snapped, invalid_count, drift = validate_snap(gdf, snapped, tol)

    assert tol == 0.001                       # metre band selected from CRS unit
    assert snapped.geometry.is_valid.all()    # topology intact
    assert invalid_count == 0
    assert list(snapped.columns) == list(gdf.columns)  # schema untouched
    assert drift <= tol * 2                    # extent stable

Run it and watch for the audit line in the log sink:

bash
pytest test_snap.py -v

Expected output:

text
test_snap.py::test_metric_snap_is_valid_and_preserves_attributes PASSED

In production, the structured record confirms the tolerance matched the CRS — a metric layer should never log a degree-band tolerance:

text
INFO geo.snap.audit {"stage": "coordinate_snapping", "epsg": 32633, "tolerance": 0.001, "feature_count": 1, "repaired_geometries": 0, "bounds_drift": 0.0}

If tolerance reads 1e-06 for a projected layer, the CRS unit was misread — check that the reproject ran before this stage and that gdf.crs is the target, not the source.

Troubleshooting

Symptom Likely cause Fix
1e-6 tolerance applied to a metric layer, welding distant vertices Snapping ran before the reproject, so the axis unit was still degree Move the snap stage after Projection Normalization Workflows; confirm gdf.crs is the target CRS at entry
RuntimeError: bounds drifted ... on otherwise clean data Tolerance too coarse for the dataset’s precision; real geometry is being collapsed Drop to the tight band (0.001 m or 1e-6 deg); never widen MAX_BOUNDS_DRIFT_CELLS to pass the gate
ValueError: GeoDataFrame lacks a defined CRS Upstream reader dropped the CRS or the source file had none Assign explicitly with .set_crs(epsg=..., allow_override=False); do not default to EPSG:4326
Mixed metre/foot outputs across one merged dataset Sources in different projections snapped independently before harmonisation Run Multi-CRS Dataset Harmonization to a single projected CRS, then snap once
Spatial join misses features that clearly overlap after snapping Stale sindex left over from the pre-snap geometry Ensure validate_snap ran — it touches snapped.sindex to force an index rebuild
make_valid changes geometry type (polygon → collection) on repair Self-intersection introduced by an over-coarse grid Tighten the tolerance so snapping no longer folds slivers, rather than relying on fail_on_invalid=False