Testing Coordinate Precision Loss During Conversion

11 min read

Coordinate precision loss during format, CRS, or serialization conversion is a silent failure mode that routinely corrupts spatial joins, invalidates topology, and introduces sub-centimeter drift in high-accuracy pipelines. The tools at the centre of this problem are concrete: numpy.float32, the PROJ transformation pipeline behind pyproj.Transformer, and the GDAL/OGR COORDINATE_PRECISION creation option. This page sits within Attribute & Metadata Checks because coordinate precision is, fundamentally, a metadata contract — the declared dtype and CRS of a coordinate array — that must be asserted, not assumed. For GIS QA engineers, data engineers, and platform/DevOps teams, detecting and preventing this drift means moving beyond naive equality assertions toward tolerance-aware tests aligned with the broader core pipeline architecture and its deterministic predicates.

Root-Cause Framing: Why Precision Bleeds

Precision degradation rarely stems from a single operation. At the engineering level it emerges from the intersection of binary storage width, transformation engines, and implicit serialization defaults — each of which can independently strip significant digits before any validation runs:

  1. Float32 vs Float64 storage. GeoJSON, Parquet, and certain PostGIS configurations default to 32-bit floats for coordinate arrays. A 32-bit IEEE 754 float carries a 23-bit mantissa, retaining roughly 7 decimal significant digits — about 1 metre of positional resolution at the equator. High-accuracy survey data, cadastral boundaries, and engineering-grade LiDAR require 64-bit floats (52-bit mantissa, ~15 digits, ~1 cm). The loss is not noise; it is a deterministic truncation of the mantissa.
  2. PROJ transformation pipelines. Default PROJ pipelines often omit high-accuracy datum-shift grids (e.g., conus, ntv2, or ntf_r93). Without explicit grid parameters or +towgs84 overrides, transformations introduce systematic offsets that compound across multi-step conversions. See the official PROJ transformation documentation for pipeline configuration standards.
  3. GDAL/OGR export rounding. Vector drivers apply implicit coordinate truncation during serialization. The COORDINATE_PRECISION creation option defaults to 7 for GeoJSON (per RFC 7946) and 15 for Shapefile, but downstream parsers — especially JavaScript-based or lightweight ETL tools — frequently re-serialize at lower precision. Refer to the GDAL GeoJSON driver specifications for creation-option overrides.
  4. Implicit ETL type coercion. GeoPandas, Dask, and PyArrow may silently downcast coordinate arrays during concatenation, partitioning, or schema inference to optimize memory bandwidth. This strips precision before validation occurs, making the drift invisible until topology checks fail downstream.

Because the loss is deterministic, it is also testable: given a known input coordinate and a known conversion, the maximum tolerable deviation can be expressed as a fixed bound and asserted in CI.

Precision Reference and Tolerance Bounds

The table below maps the storage and conversion knobs that govern precision to their effective ground resolution and a recommended assertion threshold. Use it to pick the tolerance constant your pytest checks compare against.

Mechanism Controlling parameter Mantissa / digits Effective resolution Recommended test tolerance
32-bit float storage numpy.dtype('float32') 23-bit / ~7 digits ~1 m at equator reject for cadastral; flag if dtype != float64
64-bit float storage numpy.dtype('float64') 52-bit / ~15 digits ~1 cm baseline reference
GeoJSON serialization COORDINATE_PRECISION (RFC 7946 default 7) 7 decimal places ~1.1 cm at equator >= 7 for survey-grade
Shapefile serialization COORDINATE_PRECISION (default 15) full double ~1 cm parity with source
Datum shift omission PROJ grid / +towgs84 n/a 0.1–2 m systematic < 0.05 m for engineering work

For a coordinate value xx stored as float64 and downcast to float32, the absolute drift is bounded by the machine epsilon of the narrower type:

Δx=∣x−f32(x)∣≤12 ε32 ∣x∣,ε32=2−23≈1.19×10−7\Delta x = |x - \mathrm{f32}(x)| \le \tfrac{1}{2}\,\varepsilon_{32}\,|x|, \qquad \varepsilon_{32} = 2^{-23} \approx 1.19\times10^{-7}

In Web Mercator (EPSG:3857) a mid-latitude easting is on the order of 10710^{7} metres, so the relative epsilon multiplies up to roughly 107×12ε32≈0.610^{7}\times\tfrac{1}{2}\varepsilon_{32}\approx 0.6 m of worst-case drift — which is why projected coordinates are far more sensitive to float32 storage than geographic degrees. This is the core reason the verification step below asserts against an absolute metre tolerance rather than a digit count.

Coordinate precision loss data flow with a tolerance gate A horizontal pipeline. A float64 source coordinate, the accuracy contract, flows through three boxes in series. Stage one is a PROJ transform with no datum grid, contributing 0.1 to 2 metres of systematic offset. Stage two is a float32 downcast in Arrow or Parquet, contributing roughly 0.6 metres of worst-case drift in projected coordinates. Stage three is GeoJSON serialization at 7 decimal places per RFC 7946, contributing about 1.1 centimetres of quantization. The cumulative drift arrives at a diamond-shaped tolerance gate that tests whether it is less than or equal to 0.01 metres, splitting into a PASS branch, keep float64 and stay within contract, and a FAIL branch, drift breaches the gate. float64 source coordinate declared accuracy contract Δ = 0 (reference) PROJ transform no datum grid +0.1–2 m systematic float32 downcast Arrow / Parquet +~0.6 m projected GeoJSON write 7 dp · RFC 7946 +~1.1 cm quantization Tolerance ≤ 0.01 m ? PASS within contract FAIL breaches gate

Step-by-Step Implementation

The following pattern isolates the exact deviation threshold that triggers topology failures and integrates directly into a pytest suite. It is written against Shapely 2.x, GeoPandas 0.14+, pyproj 3.6+, and pytest 7+.

Step 1 — Establish a float64 reference geometry. Pin the input as float64 explicitly so the source of truth cannot be silently downcast by NumPy’s default promotion rules.

import numpy as np
import geopandas as gpd
from shapely.geometry import Point
from pyproj import Transformer

# float64 reference — the accuracy contract anchors here
ref_coords = np.array([[-122.4194155, 37.7749295]], dtype=np.float64)
ref_gdf = gpd.GeoDataFrame(
    geometry=[Point(ref_coords[0])], crs="EPSG:4326"
)

Step 2 — Transform with an explicit, pinned PROJ pipeline. Use always_xy=True to avoid axis-order surprises and capture the full-precision projected coordinate before any storage step.

transformer = Transformer.from_crs(
    "EPSG:4326", "EPSG:3857", always_xy=True
)
x64, y64 = transformer.transform(ref_coords[0, 0], ref_coords[0, 1])

Step 3 — Simulate the lossy storage step. Downcasting to float32 reproduces what Parquet/Arrow pipelines do under default schema inference.

x32 = np.float32(x64)
y32 = np.float32(y64)

Step 4 — Quantify the drift in projection units (metres). Compute the per-axis absolute delta and take the worst case.

drift_x = abs(x64 - float(x32))
drift_y = abs(y64 - float(y32))
max_drift = max(drift_x, drift_y)

Step 5 — Compare against a configured tolerance. Load the threshold from a data contract rather than hard-coding it, mirroring how spatial tolerance thresholds are configured for every geometric assertion on this site.

TOLERANCE_METERS = 0.01  # cadastral-grade; load from YAML/JSON contract in CI
assert max_drift <= TOLERANCE_METERS, (
    f"Precision drift {max_drift:.6f} m exceeds {TOLERANCE_METERS} m"
)

Running this against a mid-latitude point reveals that float32 serialization in Web Mercator space typically introduces ~0.05–0.15 m of drift, immediately flagging datasets that claim sub-centimeter accuracy they cannot actually deliver.

How much position each representation can hold

Precision loss is not mysterious once the representable step size is written down next to the coordinate magnitude. A float32 has roughly seven significant decimal digits; a projected easting in a UTM zone is six or seven digits before the decimal point, which leaves essentially nothing after it. That single arithmetic fact explains most “the data moved slightly” reports.

Representation Significant digits Step at a 500 000 m easting Step at 50° latitude in degrees
float64 ~15–17 sub-nanometre sub-nanometre
float32 ~7 ~0.03 m ~4 m
GeoJSON, 7 decimals fixed 7 dp ~1 cm equivalent ~1.1 cm
GeoJSON, 6 decimals fixed 6 dp ~10 cm equivalent ~11 cm
Integer grid, 1 cm exact on grid 0.01 m not applicable
Where a float32's seven digits go, by CRS Two number layouts compared. The projected case shows the value five hundred thousand point one two, with six digit positions consumed before the decimal separator and only about one remaining after it within the float32 budget; the resulting representable step is roughly three centimetres and grows as the easting grows. The geographic case shows fifty point one two three four five, with two digit positions before the separator and five after; the absolute step is tiny but, because one degree spans about a hundred and eleven kilometres, it corresponds to metres on the ground. A closing note states that the same downcast produces different ground error in different coordinate reference systems, so the test must be written against the CRS actually in use. float32 budget: ~7 significant digits Projected easting 5 0 0 0 0 0 . 1 2 6 digits spent here ~1 left step ≈ 0.03 m — and grows with the easting Geographic latitude 5 0 . 1 2 3 4 5 2 spent 5 left step tiny in degrees — but ≈ 4 m on the ground The same downcast produces different ground error in different coordinate reference systems, and in opposite directions from what intuition suggests: the projected case looks worse on paper and is usually better in practice, because a degree is a very large unit. Consequence for the test: assert the drift in metres after projecting, never in the storage units, or the budget means nothing.

The table also explains why the GeoJSON rows are constant while the float32 rows are not. Fixed decimal serialisation quantises absolutely — seven decimal places is about a centimetre of latitude anywhere on Earth — whereas floating-point quantises relatively, so the error grows with the coordinate’s magnitude. A pipeline that reprojects into a coordinate system with large false eastings and then downcasts is combining the two worst cases, and it is the specific combination behind most reports of features drifting by a few centimetres for no visible reason.

Verification Pattern

Wrap the logic as a parameterized test so it runs as a deterministic CI gate. The following is copy-button-ready and exits non-zero the moment any sampled coordinate breaches the contract:

import numpy as np
import pytest
from pyproj import Transformer

TRANSFORMER = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
TOLERANCE_METERS = 0.01

@pytest.mark.parametrize("lon,lat", [
    (-122.4194155, 37.7749295),   # San Francisco
    (151.2092955, -33.8688197),   # Sydney
    (139.6917060, 35.6894875),    # Tokyo
])
def test_float32_storage_preserves_tolerance(lon, lat):
    x64, y64 = TRANSFORMER.transform(lon, lat)
    drift = max(abs(x64 - float(np.float32(x64))),
                abs(y64 - float(np.float32(y64))))
    assert drift <= TOLERANCE_METERS, f"{drift:.6f} m > {TOLERANCE_METERS} m"

Run it directly with pytest -q test_precision.py, or as a one-liner smoke check in a pre-commit hook:

python -c "import numpy as np; from pyproj import Transformer; t=Transformer.from_crs('EPSG:4326','EPSG:3857',always_xy=True); x,y=t.transform(-122.4194155,37.7749295); print(max(abs(x-float(np.float32(x))),abs(y-float(np.float32(y)))))"

A non-zero exit (or a printed value above your tolerance) is the signal to keep coordinates in float64 end-to-end or to widen the accuracy contract deliberately rather than by accident.

Attributing the loss to a stage

Knowing that a coordinate moved is not actionable; knowing which conversion moved it is. Because the losses are independent and roughly additive, a single instrumented pass through the pipeline attributes each contribution without any guesswork — measure once per stage against the float64 reference rather than only at the end.

Per-stage contribution to total coordinate drift A horizontal stacked bar decomposing total drift by stage. The largest segment is the datum transform performed without its grid shift file, contributing roughly one metre and described as systematic, meaning every feature moves by the same vector. A much smaller segment is the float32 downcast at roughly three centimetres, described as magnitude-scaled. The smallest segment is serialisation at seven decimal places, roughly one centimetre, described as fixed quantisation. A dashed vertical line marks the contract tolerance, which the total exceeds. Annotations note that the contributions are independent and additive, so per-stage measurement identifies the dominant term at once, and that only the dominant term is worth addressing. Total drift, decomposed by stage datum transform, no grid file ≈ 1 m · systematic · same vector for every feature float32 ≈ 3 cm · scales 7 dp ≈ 1 cm contract tolerance The total already exceeds the budget before the second and third stages contribute anything at all. Attribution rule: measure against the float64 reference after every stage, not once at the end. Because the terms are independent and additive, one instrumented run names the dominant one — and only the dominant one is worth optimising. Removing the float32 downcast here would change the total by three per cent; installing the grid file removes ninety-six per cent of it. A team that starts with the cheapest fix rather than the largest term spends a sprint moving the number by nothing anyone can measure.

Failure Modes and Edge Cases

Precision testing has spatial corner cases that a naive metre-tolerance check will miss. Cover these explicitly:

  1. Anti-meridian wrap (±180° longitude). Float32 storage near the ±180° seam can flip a coordinate across the meridian after rounding, turning a 1 cm drift into a globe-spanning artifact once the geometry is reprojected. Test points at 179.9999999 and -179.9999999 and assert the sign is preserved, not just the magnitude.
  2. Polar CRS and high northings. In polar stereographic or UTM zones at extreme latitudes, the projected coordinate magnitude grows, and the relative epsilon analysis above means absolute drift scales with it. A tolerance tuned for mid-latitude Web Mercator may silently pass dangerous drift near the poles — parameterize the tolerance per CRS unit and extent.
  3. Empty and null geometries. A Point() with no coordinates, or a null geometry slot in a GeoDataFrame, raises on .transform() or returns inf/nan. Guard with geom.is_empty and assert that null handling is explicit so a missing coordinate is never mistaken for a zero-drift pass.
  4. Mixed Z/M coordinates. Downcasting a 3D or measured geometry can corrupt the Z or M ordinate while leaving X/Y intact, so an X/Y-only drift check reports a false pass. Extend the delta computation across all ordinates returned by shapely.get_coordinates(geom, include_z=True).
  5. Repeated round-trips. A single conversion may sit within tolerance, but ETL stages that reproject and re-serialize repeatedly accumulate drift. Assert against the original float64 reference after the full pipeline, not against the previous stage, so error compounding is caught. This is the same parity discipline applied when comparing GeoJSON vs Shapefile outputs across formats.

Conclusion

Coordinate precision loss is not a theoretical edge case; it is a deterministic consequence of mantissa width, transformation defaults, and implicit type coercion — and because it is deterministic, it is fully testable with fixed tolerance bounds. By pinning float64, configuring explicit PROJ pipelines, and gating drift in CI, engineering teams eliminate silent spatial degradation before it corrupts analytics or breaks production topology. For the broader contract layer this check belongs to, return to Attribute & Metadata Checks.