This guide is part of Automated Vector & Raster Cleaning Workflows.

Handling Precision and Coordinate Rounding in Python Geospatial ETL

Why Floating-Point Precision Breaks Geospatial Pipelines

Raw coordinate data rarely arrives in a clean, deterministic state. GPS traces, CAD exports, and legacy GIS dumps frequently carry 12–15 significant decimal digits — far beyond any instrument’s actual measurement accuracy. That excess precision is not signal; it is floating-point noise that propagates into every downstream operation.

The consequences are concrete and costly. Spatial joins produce false negatives when shared boundary vertices differ by 1e-12 degrees. Topology checks flag valid polygons as self-intersecting due to rounding artifacts introduced during reprojection. Storage footprints balloon because binary geometry encodings retain every significant digit. Spatial index (R-tree) build times increase as the engine must accommodate a wider bounding box per vertex cluster.

Naive solutions — truncating with Python’s round(), casting coordinates to strings and back, or applying np.round() independently to X and Y arrays — each fail in different ways. The correct approach is topology-aware grid snapping, which collapses near-identical vertices to a uniform lattice rather than rounding them independently.

Coordinate Precision Pipeline Five-stage pipeline diagram showing the flow from raw coordinate data through precision audit, CRS projection, grid snapping with shapely.set_precision, post-round validation, and optimised output. Raw Input GPS / CAD / legacy GIS 12–15 decimal digits Precision Audit shapely. get_coordinates() measure decimal depth CRS Projection to_crs(target) metric units required before snapping Grid Snapping set_precision( grid_size=0.001) topology preserved Validate & Load is_valid() check GeoPackage / GeoParquet Stage 1 Stage 2 Stage 3 Stage 4 Stage 5

Prerequisites & Environment

Ensure your environment meets the following baseline before implementing coordinate rounding:

  • Python 3.9+
  • shapely>=2.0.0 — the set_precision() vectorized API is Shapely 2.0+; earlier versions lack grid-snapping support
  • geopandas>=0.13.0 with its bundled GEOS ≥ 3.9
  • numpy for vectorized numeric operations
  • Understanding of IEEE 754 — Python stores coordinates as 64-bit floats (≈15–17 significant digits). Excess precision reflects measurement noise, not accuracy

Install dependencies:

pip install "geopandas>=0.13.0" "shapely>=2.0.0" numpy pyproj

Check your GEOS version at runtime:

import shapely
print(shapely.geos_version_string)  # e.g. "3.12.1"

Always apply CRS Normalization Across Mixed Datasets before precision transformations — rounding in unprojected decimal degrees (EPSG:4326) introduces variable ground distortion across latitudes.

Version and Compatibility Matrix

Shapely GeoPandas GEOS Recommended rounding method Caveats
≥ 2.0 ≥ 0.13 ≥ 3.9 shapely.set_precision() vectorized None — preferred path
≥ 2.0 0.12 ≥ 3.9 shapely.set_precision() with .values .to_crs() returns GeoSeries, extract .values manually
1.8 0.11 ≥ 3.8 .apply(lambda g: snap(g, snap_geom, tol)) No vectorized path; slow on large datasets
< 1.7 < 0.10 < 3.8 Manual coordinate mapping via mapping() set_precision unavailable; topology not guaranteed

Step 1: Baseline Precision Audit

Before applying any transformation, quantify existing coordinate precision. Auditing prevents silent degradation of survey-grade data and identifies datasets that genuinely need intervention.

import geopandas as gpd
import numpy as np
import shapely

def audit_coordinate_precision(
    gdf: gpd.GeoDataFrame,
    sample_size: int = 1000,
) -> dict:
    """
    Samples vertices across all geometry types and estimates decimal depth.
    Uses string formatting for audit accuracy — float repr can hide trailing zeros.
    """
    if gdf.empty:
        return {"error": "Empty GeoDataFrame"}

    sampled = gdf.geometry.sample(min(sample_size, len(gdf)), random_state=42)

    # get_coordinates() handles Point, LineString, Polygon, Multi*, GeometryCollection
    # Returns an (N, 2) float64 array of all vertices
    coords = shapely.get_coordinates(sampled.values)
    x_coords = coords[:, 0]
    y_coords = coords[:, 1]

    def count_decimals(val: float) -> int:
        s = f"{val:.16f}"
        return len(s.split(".")[1].rstrip("0"))

    x_dec = np.array([count_decimals(v) for v in x_coords])
    y_dec = np.array([count_decimals(v) for v in y_coords])

    return {
        "x_max_decimals": int(np.max(x_dec)),
        "y_max_decimals": int(np.max(y_dec)),
        "x_mean_decimals": round(float(np.mean(x_dec)), 2),
        "y_mean_decimals": round(float(np.mean(y_dec)), 2),
        "total_vertices_sampled": len(coords),
        "total_geometries": len(gdf),
    }

Interpret results as follows: mean decimals above 6 in a metric CRS (meters) means sub-millimeter noise with no analytical value. Mean decimals above 10 in a geographic CRS (degrees) means sub-nanometer-scale noise. Both indicate rounding is appropriate. If max decimals are low (≤ 3) and mean is low (≤ 2), the dataset may already be pre-rounded — verify before applying a second pass.

Step 2: CRS Projection to a Metric System

Always reproject to an appropriate projected CRS before rounding. A grid_size=0.001 in degrees (EPSG:4326) corresponds to roughly 111 meters at the equator — far too coarse for most use cases. The same grid_size in EPSG:3857 (meters) equals 1 mm.

def project_for_rounding(
    gdf: gpd.GeoDataFrame,
    target_crs: str,
) -> tuple[gpd.GeoDataFrame, bool]:
    """
    Projects to target_crs if needed. Returns (projected_gdf, was_reprojected).
    Raises if source CRS is undefined — undefined CRS cannot be safely reprojected.
    """
    if gdf.crs is None:
        raise ValueError(
            "GeoDataFrame has no CRS. Assign a source CRS before reprojecting."
        )

    if gdf.crs.to_epsg() == int(target_crs.split(":")[-1]):
        return gdf, False

    return gdf.to_crs(target_crs), True

For datasets covering multiple UTM zones, reproject into a continental equal-area projection (e.g. EPSG:3035 for Europe, EPSG:5070 for the contiguous US) rather than a single UTM zone. This keeps the ground resolution of grid_size consistent across all features.

Step 3: Topology-Safe Grid Snapping with shapely.set_precision

Shapely 2.0’s set_precision() implements a robust grid-snapping algorithm that collapses near-identical coordinates uniformly. The function operates on NumPy arrays of geometries — pass gdf.geometry.values (a GeometryArray) rather than iterating with .apply().

import shapely
import geopandas as gpd
import logging

log = logging.getLogger(__name__)

def snap_to_grid(
    gdf: gpd.GeoDataFrame,
    grid_size: float = 0.001,
) -> gpd.GeoDataFrame:
    """
    Snaps all geometry coordinates to a uniform grid.

    Args:
        gdf: Input GeoDataFrame in a projected (metric) CRS.
        grid_size: Tolerance in CRS units. For metric CRS: 0.001 = 1 mm.

    Returns:
        New GeoDataFrame with snapped geometries.

    Raises:
        ValueError: If input contains invalid geometries before snapping.
    """
    invalid_before = (~shapely.is_valid(gdf.geometry.values)).sum()
    if invalid_before > 0:
        raise ValueError(
            f"{invalid_before} invalid geometries detected. "
            "Run geometry repair before coordinate rounding."
        )

    gdf_snapped = gdf.copy()
    # Vectorized: single C-level call across the entire geometry array
    gdf_snapped.geometry = shapely.set_precision(
        gdf.geometry.values,
        grid_size=grid_size,
    )

    invalid_after = (~shapely.is_valid(gdf_snapped.geometry.values)).sum()
    if invalid_after > 0:
        log.warning(
            "%d geometries became invalid after grid snapping at grid_size=%s. "
            "Consider reducing grid_size or running make_valid() on the output.",
            invalid_after,
            grid_size,
        )

    log.info(
        "Grid snapping complete: %d features, grid_size=%s, "
        "%d invalid post-snap.",
        len(gdf_snapped),
        grid_size,
        invalid_after,
    )
    return gdf_snapped

The mode parameter controls how set_precision handles topological conflicts. The default "valid_output" mode attempts to preserve validity. Use "pointwise" only when you need coordinates snapped without any topological reorganisation — this can produce self-intersections in complex polygons. See the Shapely set_precision reference for the full mode list.

If snap_to_grid raises because of pre-existing invalid geometries, resolve those first with Geometry Repair with Shapely & GeoPandas before returning to this step.

Step 4: Post-Rounding Validation

Rounding is a destructive operation. Validation must confirm spatial relationships remain intact, area drift is within tolerance, and output geometry counts match inputs.

def validate_rounding_impact(
    original: gpd.GeoDataFrame,
    rounded: gpd.GeoDataFrame,
    max_area_drift_pct: float = 0.01,
) -> dict:
    """
    Compares spatial metrics before and after rounding.
    Returns a result dict; raises AssertionError if thresholds are breached.
    """
    orig_area = original.geometry.area.sum()
    round_area = rounded.geometry.area.sum()
    area_drift_pct = (
        abs(orig_area - round_area) / orig_area * 100 if orig_area > 0 else 0.0
    )

    all_valid = bool(shapely.is_valid(rounded.geometry.values).all())
    count_match = len(original) == len(rounded)
    crs_aligned = str(original.crs) == str(rounded.crs)

    if area_drift_pct > max_area_drift_pct:
        raise AssertionError(
            f"Area drift {area_drift_pct:.4f}% exceeds threshold "
            f"{max_area_drift_pct}%. Reduce grid_size or investigate outliers."
        )

    return {
        "geometry_count_match": count_match,
        "area_drift_pct": round(area_drift_pct, 6),
        "topology_valid": all_valid,
        "crs_aligned": crs_aligned,
        "passed": count_match and all_valid and crs_aligned,
    }

Acceptable area drift is typically below 0.01% for administrative boundaries and 0.1% for coarse environmental polygons. If drift exceeds thresholds, reduce grid_size by one order of magnitude and re-run. Repeated excess drift signals either pathological geometry or a CRS mismatch — check that the area comparison happens in the same projected CRS.

Step 5: ETL Pipeline Integration and Logging

Embedding precision control into batch ETL requires deterministic logging, ordered chaining with other cleaning steps, and graceful fallback handling.

import logging
import geopandas as gpd
from pathlib import Path

log = logging.getLogger(__name__)

def process_coordinate_precision(
    input_path: Path,
    output_path: Path,
    target_crs: str,
    grid_size: float = 0.001,
    max_area_drift_pct: float = 0.01,
) -> dict:
    """
    Full pipeline: read → audit → project → snap → validate → write.
    Returns metrics dict for upstream orchestrators (Airflow, Prefect, Dagster).
    """
    metrics: dict = {
        "input_path": str(input_path),
        "output_path": str(output_path),
        "grid_size": grid_size,
        "processed_features": 0,
        "invalid_post_snap": 0,
        "area_drift_pct": None,
        "storage_before_mb": round(input_path.stat().st_size / (1024 ** 2), 3),
        "storage_after_mb": None,
        "passed_validation": False,
    }

    gdf = gpd.read_file(input_path)
    log.info("Loaded %d features from %s", len(gdf), input_path)

    # 1. Audit current precision
    audit = audit_coordinate_precision(gdf)
    log.info("Precision audit: %s", audit)

    # 2. Project to metric CRS
    gdf, reprojected = project_for_rounding(gdf, target_crs)
    if reprojected:
        log.info("Reprojected to %s", target_crs)

    # 3. Grid snap
    gdf_snapped = snap_to_grid(gdf, grid_size=grid_size)

    # 4. Validate
    result = validate_rounding_impact(gdf, gdf_snapped, max_area_drift_pct)
    metrics["area_drift_pct"] = result["area_drift_pct"]
    metrics["passed_validation"] = result["passed"]
    metrics["invalid_post_snap"] = int(
        (~shapely.is_valid(gdf_snapped.geometry.values)).sum()
    )

    # 5. Write
    gdf_snapped.to_file(output_path, driver="GPKG")
    metrics["storage_after_mb"] = round(output_path.stat().st_size / (1024 ** 2), 3)
    metrics["processed_features"] = len(gdf_snapped)

    log.info("Pipeline complete: %s", metrics)
    return metrics

Within a broader cleaning chain, coordinate rounding must occur after attribute normalisation and CRS alignment but before spatial deduplication and indexing. Running deduplication before rounding allows floating-point duplicates to slip through as distinct features. For spatial deduplication logic that integrates cleanly after this step, see Spatial Deduplication & Topology Simplification.

Advanced Patterns and Edge Cases

Handling Multipart and Mixed Geometry Collections

shapely.set_precision() operates correctly on MultiPolygon, MultiLineString, and GeometryCollection types — no special handling is needed. However, grid snapping can convert a MultiPolygon into a Polygon when component parts collapse onto a shared boundary. Check for type changes post-snap if your schema enforces a fixed geometry type:

import shapely
import numpy as np

def check_geometry_type_drift(
    original: gpd.GeoDataFrame,
    snapped: gpd.GeoDataFrame,
) -> int:
    """Returns count of features whose geometry type changed after snapping."""
    orig_types = np.array([g.geom_type for g in original.geometry.values])
    snap_types = np.array([g.geom_type for g in snapped.geometry.values])
    changed = int((orig_types != snap_types).sum())
    if changed:
        log.warning("%d features changed geometry type after snapping.", changed)
    return changed

Attribute Mapping Before Rounding

Coordinate rounding should always follow schema harmonisation, not precede it. Joining on string-formatted coordinate keys (a common workaround for float equality) breaks if the schema renames or drops the key column mid-pipeline. See Attribute Mapping & Schema Harmonization for the correct ordering of these transforms.

Precision Snapping in Chunked I/O for Large Datasets

Datasets exceeding available RAM require chunked processing. GeoPandas’ read_file supports the rows parameter for slice-based iteration; pair it with pyogrio for faster I/O:

import geopandas as gpd
import shapely
from pathlib import Path

def snap_chunked(
    input_path: Path,
    output_path: Path,
    target_crs: str,
    grid_size: float = 0.001,
    chunk_size: int = 50_000,
) -> int:
    """Processes large vector files in row-sliced chunks to manage memory."""
    import fiona

    total_written = 0
    schema = None

    with fiona.open(input_path) as src:
        total_rows = len(src)
        offsets = range(0, total_rows, chunk_size)

    for i, offset in enumerate(offsets):
        chunk = gpd.read_file(
            input_path,
            rows=slice(offset, offset + chunk_size),
            engine="pyogrio",
        )
        chunk, _ = project_for_rounding(chunk, target_crs)
        chunk = snap_to_grid(chunk, grid_size=grid_size)

        write_mode = "w" if i == 0 else "a"
        chunk.to_file(output_path, driver="GPKG", mode=write_mode)
        total_written += len(chunk)
        log.info("Processed chunk %d/%d (%d features)", i + 1, len(list(offsets)), len(chunk))

    return total_written

Performance Optimization: Vectorized API vs. .apply()

The performance gap between Shapely 2.0’s vectorized API and Shapely 1.x’s .apply() approach is substantial. shapely.set_precision() dispatches to a single C-level GEOS operation across the full geometry array, while .apply(lambda g: ...) incurs Python function-call overhead per feature.

import time
import numpy as np
import shapely
import geopandas as gpd

def benchmark_precision_methods(gdf: gpd.GeoDataFrame, grid_size: float = 0.001) -> dict:
    """Compares vectorized vs. apply-based rounding on the same dataset."""
    geom_array = gdf.geometry.values

    # Vectorized (Shapely 2.0+)
    t0 = time.perf_counter()
    _ = shapely.set_precision(geom_array, grid_size=grid_size)
    vectorized_ms = (time.perf_counter() - t0) * 1000

    # Element-wise via apply (legacy pattern — for comparison only)
    snap_geom = shapely.from_wkt("POINT (0 0)")  # dummy snap target
    t1 = time.perf_counter()
    _ = gdf.geometry.apply(lambda g: shapely.snap(g, snap_geom, grid_size))
    apply_ms = (time.perf_counter() - t1) * 1000

    return {
        "features": len(gdf),
        "vectorized_ms": round(vectorized_ms, 2),
        "apply_ms": round(apply_ms, 2),
        "speedup_x": round(apply_ms / vectorized_ms, 1) if vectorized_ms > 0 else None,
    }

In practice, the vectorized path is 20–80x faster on datasets of 10,000+ features. The speedup compounds when rounding is embedded in a loop over many input files.

Integration into ETL Pipelines

Schema Enforcement Hooks

Precision rounding must be recorded in your dataset’s processing log so downstream consumers know what tolerance was applied. Attach metadata to the output file:

def write_with_precision_metadata(
    gdf: gpd.GeoDataFrame,
    output_path: Path,
    grid_size: float,
    source_crs: str,
    target_crs: str,
) -> None:
    """Writes GeoPackage and attaches a JSON sidecar with processing metadata."""
    import json, datetime

    gdf.to_file(output_path, driver="GPKG")
    meta = {
        "grid_size": grid_size,
        "source_crs": source_crs,
        "target_crs": target_crs,
        "processed_at": datetime.datetime.utcnow().isoformat(),
        "shapely_version": shapely.__version__,
    }
    sidecar = output_path.with_suffix(".precision.json")
    sidecar.write_text(json.dumps(meta, indent=2))

Dead-Letter Queue Pattern

Features that become invalid after snapping should be routed to a dead-letter file rather than silently dropped or written with broken geometry:

def snap_with_dead_letter(
    gdf: gpd.GeoDataFrame,
    output_path: Path,
    dead_letter_path: Path,
    grid_size: float = 0.001,
) -> tuple[int, int]:
    """
    Snaps geometries and splits valid vs. invalid outputs.
    Returns (valid_count, dead_letter_count).
    """
    gdf_snapped = snap_to_grid(gdf, grid_size=grid_size)
    valid_mask = shapely.is_valid(gdf_snapped.geometry.values)

    gdf_snapped[valid_mask].to_file(output_path, driver="GPKG")

    dead = gdf_snapped[~valid_mask]
    if len(dead) > 0:
        dead.to_file(dead_letter_path, driver="GPKG")
        log.warning("%d features sent to dead-letter: %s", len(dead), dead_letter_path)

    return int(valid_mask.sum()), int((~valid_mask).sum())

CI/CD Embedding

Add a precision gate to your CI pipeline to prevent regressions — a future import of the same source dataset should not silently reintroduce high-precision noise:

python - <<'EOF'
import geopandas as gpd, sys
gdf = gpd.read_file("outputs/parcels_cleaned.gpkg")
audit = audit_coordinate_precision(gdf)
if audit["x_max_decimals"] > 4 or audit["y_max_decimals"] > 4:
    print(f"FAIL: coordinate precision exceeds threshold: {audit}")
    sys.exit(1)
print("PASS: coordinate precision within bounds.")
EOF

Best Practices and Common Pitfalls

When to Omit Rounding

Survey-grade cadastral records, legal boundary descriptions, and high-frequency GNSS telemetry often require full 64-bit float retention. Rounding these datasets can violate regulatory tolerances or erase legitimate micro-topographic variation. Always verify data provenance and any accompanying accuracy specifications before applying grid snapping.

Float Equality and Epsilon Drift

Even after rounding, Python’s underlying float representation can introduce microscopic drift during arithmetic. Use shapely.equals_exact(a, b, tolerance=1e-9) or numpy.isclose() for post-rounding equality checks. Avoid == on float coordinates.

Pipeline Idempotency

Design rounding functions to be idempotent: running set_precision twice with the same grid_size on already-snapped data must return an identical result. This prevents double-rounding in retry-heavy orchestration frameworks and ensures audit logs are deterministic. Verify idempotency in your test suite:

def test_snap_idempotency(gdf: gpd.GeoDataFrame, grid_size: float = 0.001) -> bool:
    first = snap_to_grid(gdf, grid_size=grid_size)
    second = snap_to_grid(first, grid_size=grid_size)
    return bool(
        shapely.equals_exact(first.geometry.values, second.geometry.values, tolerance=0).all()
    )

Troubleshooting Reference

Failure Mode Root Cause Mitigation
TopologyException during set_precision Input polygon has pre-existing self-intersection Run shapely.make_valid() before grid snapping
Geometry type changes after snapping Adjacent rings collapse onto shared edge Reduce grid_size by 10×; check for degenerate inputs
Area drift exceeds threshold grid_size too coarse relative to feature size Use smaller grid_size; flag features smaller than 10× grid_size
Post-snap ST_Intersects false negatives Rounding in geographic CRS before projection Always project to metric CRS before calling set_precision
Duplicate vertices remain after snapping Using mode="pointwise" Switch to default mode="valid_output" for topology-aware snapping