Spatial Epi Public Health GIS

Controlling Floating-Point Drift in Coordinate Pipelines

This guide, part of Precision Standards in Epi-Mapping, deals with the error that accumulates not from any one step but from many: the slow wander of coordinates as they pass through reprojection, text serialisation, database round-trips and geometry operations.

Problem Context & Constraints

A double-precision float represents a projected coordinate in the millions of metres to around a nanometre. No single operation in a spatial pipeline loses meaningful precision. The problem is composition. A case point reprojected from geographic to projected coordinates, written to a text format with eight decimals, read back, transformed to a web projection for a map tile, transformed back for analysis, snapped during an overlay, and written again has been through a chain of approximations. Each contributes a tiny error, and the errors do not cancel reliably.

For most analyses the result is irrelevant — a few millimetres. It becomes relevant in three situations that surveillance pipelines hit routinely. Points on or near boundaries. A case geocoded to a parcel centroid that lies exactly on a tract line can flip between tracts from one run to the next if its coordinates wander by a micrometre across the line, and a count moves with it. Equality tests. Deduplication by exact coordinate, cache lookups keyed on geometry, and joins on coordinate values all break when two representations of the same point differ in the last bit. Idempotency. A pipeline whose output geometry changes slightly on every re-run cannot be tested by comparing outputs, which defeats the regression discipline described in Testing Geospatial Analysis Code.

Three constraints shape the fix. Establish a canonical precision well below any meaningful spatial scale and well above floating-point noise — a millimetre or a centimetre in projected coordinates — and snap to it at defined points. Minimise round-trips between coordinate systems and between binary and text. And measure the drift, because the only way to know whether it matters is to run the chain and compare.

Small Errors Compose Along the Chain A case point passing through six steps: reprojection to the working system, serialisation to text with limited decimals, reading back, transformation to web projection for a tile, transformation back for analysis, and an overlay that snaps to a grid. The cumulative displacement from the original position is drawn above each step as a rising bar. No single step moves the point meaningfully, but the total after six steps is enough to move a point that lies on a tract boundary to the other side. A canonical snap after the first step and binary serialisation hold the cumulative displacement flat. No step is to blame; the chain is distance to the nearest tract line reproject to text read back web tile back again overlay snap crosses the line cumulative displacement from the original coordinate after each step shown exaggerated; for a point sitting on a boundary, any crossing moves a count

Prerequisites

  • python 3.11
  • shapely 2.0.4, geopandas 0.14.4, pyproj 3.6.1, numpy 1.26.4
  • A declared working coordinate reference system in metres
  • A chosen canonical precision — a grid size in metres
  • A representative sample of case points, including some known to fall on area boundaries

Step-by-Step Solution

Snap geometry to the canonical grid immediately after entering the working system, serialise in binary, avoid unnecessary transformations, and measure drift by running the chain twice.

# Canonical precision, drift measurement and boundary-stable assignment.
# Pinned: shapely==2.0.4, geopandas==0.14.4, pyproj==3.6.1, numpy==1.26.4
import logging

import numpy as np
import geopandas as gpd
import shapely

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("precision.drift")

GRID_M = 0.01          # canonical precision: one centimetre in the working CRS


def to_canonical(gdf: gpd.GeoDataFrame, work_crs: str, grid: float = GRID_M) -> gpd.GeoDataFrame:
    """Enter the working CRS once, then snap every vertex to the canonical grid."""
    g = gdf.to_crs(work_crs) if gdf.crs != work_crs else gdf.copy()
    g["geometry"] = shapely.set_precision(g.geometry.values, grid_size=grid)
    invalid = int((~g.geometry.is_valid).sum())
    if invalid:
        log.warning("%d geometries became invalid on snapping — repair before use", invalid)
    log.info("canonicalised %d geometries to a %.3f m grid in %s", len(g), grid, work_crs)
    return g


def measure_drift(points: gpd.GeoDataFrame, chain) -> dict:
    """Run a transformation chain and report how far points moved."""
    before = points.geometry
    after = chain(points).geometry
    d = before.distance(after)
    out = {"max_m": float(d.max()), "p99_m": float(np.percentile(d, 99)),
           "moved_over_1mm": int((d > 0.001).sum()), "n": int(len(d))}
    log.info("drift over chain: max=%.6f m p99=%.6f m, %d of %d moved more than 1 mm",
             out["max_m"], out["p99_m"], out["moved_over_1mm"], out["n"])
    return out


def boundary_stable_assign(points: gpd.GeoDataFrame, areas: gpd.GeoDataFrame,
                           id_col: str = "geoid", tolerance: float = 0.05) -> gpd.GeoDataFrame:
    """Point-in-polygon with a deterministic rule for points on or near a boundary.

    Points within `tolerance` of an area boundary are assigned to the candidate
    with the smallest identifier, so a micrometre of drift cannot flip them.
    """
    near = gpd.sjoin(points, areas[[id_col, "geometry"]].assign(
        geometry=areas.geometry.buffer(tolerance)), how="left", predicate="within")
    counts = near.groupby(level=0)[id_col].nunique()
    ambiguous = counts[counts > 1].index
    resolved = (near.sort_values(id_col).groupby(level=0)[id_col].first())
    out = points.copy()
    out[id_col] = resolved.reindex(out.index)
    out["boundary_ambiguous"] = out.index.isin(ambiguous)
    log.info("assigned %d points; %d within %.2f m of a boundary resolved by rule",
             len(out), len(ambiguous), tolerance)
    return out

Snapping once, immediately after entering the working system, fixes the representation that every later step starts from. Two runs that read the same input then begin from bit-identical coordinates, and any later operation that is itself deterministic produces identical output — which restores the ability to regression-test geometry by exact comparison.

The boundary-stable assignment addresses the case that matters epidemiologically. A point within a few centimetres of a tract line is on that line for every practical purpose, and whether it lands in one tract or the other is a numerical accident. Resolving such points by a fixed rule — smallest identifier — makes the assignment deterministic, and flagging them makes the count of boundary cases visible, which is useful in its own right because parcel centroids on boundaries are more common than one expects.

Points on the Line Need a Rule, Not Luck Two adjacent tracts share a boundary. A case point sits on the boundary. Across three runs with slightly different floating-point paths, the point falls into the left tract, then the right, then the left, and the tract counts change each time. With a tolerance band around the boundary and a rule that assigns any point inside the band to the tract with the smaller identifier, the point is assigned to the same tract on every run and is flagged as a boundary case. Inside the tolerance band, the rule decides and the run cannot tract 101 tract 102 runs: 101, 102, 101 tract 101 tract 102 every run: 101, flagged the dashed band is the tolerance; the smaller identifier wins inside it

Validation & Edge Cases

Snapping that invalidates polygons. Snapping a polygon with a narrow neck to a coarse grid can make it self-intersect. Choose the grid well below the narrowest real feature — a centimetre is safe for administrative boundaries — and check validity after snapping.

Text formats with too few decimals. Writing geographic coordinates with six decimals truncates to about ten centimetres, which is coarser than the canonical grid and re-introduces drift on every text round-trip. Prefer binary formats; where text is unavoidable, write enough decimals to exceed the canonical precision.

Web projection round-trips. Transforming to a web mapping projection and back for analysis is a common and entirely avoidable round-trip. Keep analysis geometry in the working system and derive web geometry as a one-way output.

Coordinate equality in caches and joins. Keying a cache on raw coordinates fails when two representations differ in the last bit. Key on the canonical snapped coordinates, or better, on a stable identifier.

Different grids in different stages. Snapping to a millimetre in one stage and a centimetre in another produces a pipeline whose precision is the coarser of the two and whose intermediate geometry is inconsistent. Declare one canonical grid for the pipeline.

INFO precision.drift canonicalised 48900 geometries to a 0.010 m grid in EPSG:5070
INFO precision.drift drift over chain: max=0.000000 m p99=0.000000 m, 0 of 48900 moved more than 1 mm
INFO precision.drift assigned 48900 points; 37 within 0.05 m of a boundary resolved by rule

A grid chosen too coarse. A snapping grid of a metre looks harmless and silently merges distinct parcel centroids in dense housing, collapsing records that should stay separate and inflating apparent duplicates. The canonical grid is a computational precision, not an analytic one, and should sit two or three orders of magnitude below the finest distance any analysis uses.

Legacy data already carrying drift. Archives written before a canonical grid existed contain coordinates that differ in their last digits from the same points written today. Snap the archive once, record that it was done and with which grid, and treat the snapped archive as the new reference rather than mixing snapped and unsnapped partitions.

Choosing the grid is a matter of placing it in the wide gap between floating-point noise and the smallest distance any analysis cares about:

Where the Canonical Grid Belongs A logarithmic distance scale running from nanometres to kilometres. Floating-point noise in projected coordinates sits near the nanometre end. The smallest distances analyses care about — geocode positional error, buffer radii and parcel spacing — sit between a metre and a kilometre. The canonical grid of one centimetre is placed in the wide gap between them: far above noise, so snapping removes drift, and far below analytic scales, so snapping changes nothing any analysis can detect. A one-metre grid would intrude on parcel spacing in dense housing. Far above the noise, far below anything an analysis measures 1 nm 1 µm 1 mm 1 cm 1 m 1 km float noise analytic distances canonical grid a wide safe gap: snapping here removes drift and changes nothing measurable a one-metre grid would sit inside the analytic range and merge dense parcels

Compliance Notes

Record the canonical grid, the working coordinate reference system, the point in the pipeline where snapping occurs, the serialisation formats used, the boundary tolerance and assignment rule, and the count of boundary-ambiguous points per run.

Precision decisions interact with disclosure. A canonical grid is an internal representation choice and is far finer than any published precision; publishing coordinates at the canonical precision would be a disclosure regardless of drift. The published precision is governed separately by Setting Decimal Precision for Disease Coordinate Mapping, and keeping the two decisions distinct — one for computational stability, one for release — avoids the error of treating an internal precision as safe to publish.