Coordinate Precision and Floating Point in GIS Explained
On this page
Problem statement
Two things that should be equal are not:
from shapely.geometry import Point
print(Point(0.1 + 0.2, 1).equals(Point(0.3, 1)), 0.1 + 0.2)False 0.30000000000000004
And a coordinate that goes to longitude and latitude and back does not come home — here the first vertex of a Land Registry parcel in Rutland:
from shapely.geometry import Point as P
import geopandas as gpd
vertex = parcels.geometry.iloc[0].exterior.coords[0]
back = gpd.GeoSeries([P(vertex)], crs=27700).to_crs(4326).to_crs(27700).iloc[0]
print(vertex, "->", (back.x, back.y))(495533.11, 299678.86) -> (495533.110881145, 299678.8595994822)
The second is not even floating-point residue: the vertex moved by nearly a millimetre, because PROJ's route between British National Grid and WGS 84 (without the OSTN15 grid) does not invert exactly — every parcel in the file moves 0.96 to 0.97 mm per round trip. Neither is a bug in GeoPandas, GEOS or your data. Both are consequences of storing coordinates as IEEE-754 double-precision floating-point numbers, which is what every part of the GIS stack does. Understanding what that guarantees — and what it does not — explains sliver polygons, failed equality tests, TopologyException, and why a file can shrink by a third without losing anything real.
Quick answer
Coordinates are doubles, so:
- never compare coordinates or geometries with
==; use a tolerance - a double holds ~15–17 significant digits, so precision available depends on magnitude
- in EPSG:4326, 6 decimal places ≈ 0.1 m; in a metric CRS, 3 decimals ≈ 1 mm
- store what your data actually measures — extra digits are noise that costs space and causes slivers
- snap to a shared grid with
set_precision()before overlays that must agree
from shapely.geometry import Point
from shapely import equals_exact, set_precision
a, b = Point(0.1 + 0.2, 1), Point(0.3, 1)
print(a.equals(b)) # False — exact comparison
print(equals_exact(a, b, tolerance=1e-9)) # True — within a sane tolerance
snapped_a = set_precision(a, grid_size=1e-6)
snapped_b = set_precision(b, grid_size=1e-6)
print(snapped_a.equals(snapped_b)) # True — both land on the same grid pointFalse
True
True
The mental shift is to stop thinking of a coordinate as an exact number and start thinking of it as a measurement with a tolerance — which is what it always was on the ground.
What a double can hold
From floating-point error to metres on the ground
Why 0.1 + 0.2 is not 0.3
A double stores a number as a sign, an 11-bit exponent and a 52-bit fraction — a binary approximation. Numbers that are simple in decimal are often infinitely repeating in binary, exactly as 1/3 is in decimal.
from decimal import Decimal
print(Decimal(0.1))
# 0.1000000000000000055511151231257827021181583404541015625
print(Decimal(0.1) + Decimal(0.2) == Decimal(0.3)) # False
print(format(0.1 + 0.2, ".20f")) # 0.300000000000000044410.1000000000000000055511151231257827021181583404541015625
False
0.30000000000000004441
The error is around 1 part in 10¹⁶. That is irrelevant for a distance and decisive for an equality test, which is why every geometric comparison in a robust workflow carries a tolerance.
Precision depends on magnitude
This is the part that catches GIS specifically. A double has a fixed number of significant digits, so the absolute precision available shrinks as coordinates get larger.
import numpy as np
for value in [1.0, 100.0, 55.95, 325_000.0, 6_378_137.0, 30_000_000.0]:
print(f"{value:>12,.1f} smallest representable step: {np.spacing(value):.2e}") 1.0 smallest representable step: 2.22e-16
100.0 smallest representable step: 1.42e-14
56.0 smallest representable step: 7.11e-15
325,000.0 smallest representable step: 5.82e-11
6,378,137.0 smallest representable step: 9.31e-10
30,000,000.0 smallest representable step: 3.73e-09
Even at Earth-radius magnitudes a double resolves a nanometre, so raw precision is never the limit in GIS. What is a limit is accumulated error: every reprojection, buffer and intersection performs arithmetic, and each operation can move a coordinate in the last few bits. Do that a dozen times to two copies of the same boundary and they diverge.
What decimal places mean on the ground
def ground_resolution(decimals: int, geographic: bool = True) -> float:
"""Approximate ground distance of one unit in the last decimal place."""
return (111_320 if geographic else 1.0) / (10 ** decimals)
print("EPSG:4326 (degrees)")
for dp in range(0, 9):
print(f" {dp} dp → {ground_resolution(dp):>12,.4f} m")
print("\nProjected CRS (metres)")
for dp in range(0, 5):
print(f" {dp} dp → {ground_resolution(dp, geographic=False):>12,.4f} m")EPSG:4326 (degrees)
0 dp → 111,320.0000 m
1 dp → 11,132.0000 m
2 dp → 1,113.2000 m
3 dp → 111.3200 m
4 dp → 11.1320 m
5 dp → 1.1132 m
6 dp → 0.1113 m
7 dp → 0.0111 m
8 dp → 0.0011 m
Projected CRS (metres)
0 dp → 1.0000 m
1 dp → 0.1000 m
2 dp → 0.0100 m
3 dp → 0.0010 m
4 dp → 0.0001 m
The useful landmarks: in degrees, 5 dp is about a metre, 6 dp about 10 cm, 7 dp about 1 cm. In metres, 2 dp is a centimetre and 3 dp a millimetre. Anything finer than the survey that produced the data is noise.
Match stored precision to real accuracy
import geopandas as gpd
gdf = gpd.read_file("data/raw/parcels.gpkg")
coords = gdf.geometry.get_coordinates()
decimals = coords["x"].astype(str).str.split(".").str[1].fillna("").str.len()
print("decimal places actually present:")
print(decimals.value_counts().sort_index().tail(6).to_string())
print(f"\ncoordinates: {len(coords):,}")
round_trip = gdf.to_crs(4326).to_crs(27700).geometry.get_coordinates()
print("after one trip through EPSG:4326:",
round_trip["x"].astype(str).str.split(".").str[1].str.len().value_counts().sort_index().tail(4).to_dict())decimal places actually present:
x
1 111901
2 568875
3 170045
coordinates: 850,821
after one trip through EPSG:4326: {8: 4304, 9: 44533, 10: 446425, 11: 355074}
The Land Registry publishes these parcels to at most three decimal places — a millimetre. One trip through longitude and latitude and back, and most coordinates carry ten or eleven, and have moved a millimetre besides.
A layer whose coordinates carry ten decimal places did not come from a survey with sub-nanometre accuracy; it came from a reprojection. The extra digits are computational residue, and on the original data rounding them away loses nothing:
from shapely import set_precision
metric = gdf.to_crs(gdf.estimate_utm_crs())
before = metric.geometry.area.sum()
metric["geometry"] = set_precision(metric.geometry.values, grid_size=0.01) # 1 cm
after = metric.geometry.area.sum()
print(f"area change after snapping to 1 cm: {(after - before) / before:+.6%}")area change after snapping to 1 cm: -0.000005%
Comparing geometries safely
from shapely import equals_exact
from shapely.geometry import Point
import geopandas as gpd
a, b = Point(0.1 + 0.2, 1), Point(0.3, 1)
print(a == b) # False — object identity/structure
print(a.equals(b)) # False — exact coordinate equality
print(equals_exact(a, b, tolerance=1e-9)) # True — the one to use
print(a.dwithin(b, 1e-9)) # True — distance-based
# for whole layers: the parcels against their own round trip through EPSG:4326,
# which moves every vertex about 1 mm
gdf_a = gdf
gdf_b = gdf.to_crs(4326).to_crs(gdf.crs)
print(f"{gdf_a.geometry.geom_equals_exact(gdf_b.geometry, tolerance=0).sum()} of {len(gdf_a)} identical")
for tolerance in (1e-6, 2e-3):
same = gdf_a.geometry.geom_equals_exact(gdf_b.geometry, tolerance=tolerance)
print(f"{same.sum()} of {len(same)} geometries match within {tolerance:g}")False
False
True
True
0 of 21253 identical
0 of 21253 geometries match within 1e-06
21253 of 21253 geometries match within 0.002
equals() is a topological test that still relies on exact coordinates for its predicates; equals_exact(tolerance=…) is the tolerant comparison you almost always want in a test or a diff.
Snap before operations that must agree
import geopandas as gpd
from shapely import set_precision
def prepare(gdf, grid_size=0.001):
metric = gdf.to_crs(gdf.estimate_utm_crs())
metric["geometry"] = set_precision(metric.geometry.values, grid_size=grid_size)
metric["geometry"] = metric.geometry.make_valid()
return metric[metric.geometry.notna() & ~metric.geometry.is_empty]
parcels_prepared = prepare(gpd.read_file("data/raw/parcels.gpkg"))
zones = prepare(gpd.read_file("data/ref/zones.gpkg").to_crs(parcels_prepared.crs))
overlay = gpd.overlay(parcels_prepared, zones, how="intersection")
print(len(overlay), "result features")21870 result features
Snapping both inputs to the same grid is what turns "these boundaries are nearly the same" into "these boundaries are the same", which is the precondition GEOS's overlay algorithms actually want. It is not a cure for every sliver, though. Neighbouring Rutland parcels overlap each other — not by float noise, but by what the digitising left:
def neighbour_overlaps(frame):
pairs = gpd.sjoin(frame, frame[["geometry"]], predicate="intersects")
pairs = pairs[pairs.index < pairs["index_right"]]
areas = frame.geometry.loc[pairs.index].intersection(
frame.geometry.loc[pairs["index_right"]].set_axis(pairs.index), align=False).area
return len(pairs), int((areas > 0).sum()), float(areas[areas > 0].median())
sample = gdf.iloc[:3000]
for grid in (None, 0.01, 0.1):
frame = sample if grid is None else sample.set_geometry(set_precision(sample.geometry.values, grid_size=grid))
n, overlapping, median = neighbour_overlaps(frame)
print(f"grid {str(grid):5}: {overlapping:,} of {n:,} touching pairs overlap, median {median:.4f} m²")grid None : 1,798 of 6,229 touching pairs overlap, median 0.0071 m²
grid 0.01 : 764 of 6,234 touching pairs overlap, median 2.9086 m²
grid 0.1 : 640 of 6,258 touching pairs overlap, median 23.2000 m²
A centimetre grid removes more than half of them — those were disagreements smaller than a centimetre — and what is left is larger, median 2.9 m²: two surveys of the same fence that genuinely disagree. A coarser grid keeps removing small overlaps while its larger moves create new ones, so the median climbs to 23 m². Snapping fixes disagreement below the grid size. Disagreement in the data needs the slivers-and-gaps treatment.
Know when precision is not the problem
import geopandas as gpd
gdf = gpd.read_file("data/raw/parcels.gpkg")
print("CRS:", gdf.crs)
print("bounds:", [round(float(v), 2) for v in gdf.total_bounds])CRS: EPSG:27700
bounds: [479000.61, 292380.57, 506878.39, 319200.7]
If features are metres out of place, that is precision. If they are hundreds of metres out, it is a datum or transformation issue; if they are hundreds of kilometres out, it is the wrong CRS entirely. Floating point is a millimetre-scale phenomenon — anything bigger has a different cause.
A precision audit of Rutland's parcels
The data used on this page
HM Land Registry's INSPIRE Index Polygons for Rutland (21,253 freehold parcels, Open Government Licence, © Crown copyright and database rights), in British National Grid, and the ONS wards of Rutland as the zones for the overlay example. The download service sets a cookie before it serves the file, so it is fetched with a requests.Session.
import io
import zipfile
from pathlib import Path
import geopandas as gpd
import requests
for d in ("data/raw", "data/ref"):
Path(d).mkdir(parents=True, exist_ok=True)
if not Path("data/raw/parcels.gpkg").exists():
reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip",
timeout=300)
reply.raise_for_status()
zipfile.ZipFile(io.BytesIO(reply.content)).extract("Land_Registry_Cadastral_Parcels.gml")
gpd.read_file("Land_Registry_Cadastral_Parcels.gml")[["INSPIREID", "geometry"]].to_file("data/raw/parcels.gpkg")
gpd.read_file("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
"Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query"
"?where=LAD24CD%3D%27E06000017%27&outFields=WD24NM&outSR=27700&f=geojson").to_file("data/ref/zones.gpkg")
parcels = gpd.read_file("data/raw/parcels.gpkg")
print(len(parcels), "parcels,", len(gpd.read_file("data/ref/zones.gpkg")), "wards,", parcels.crs)21253 parcels, 15 wards, EPSG:27700
Example 1: a precision audit for a layer
import geopandas as gpd
import numpy as np
def precision_audit(gdf: gpd.GeoDataFrame) -> dict:
coords = gdf.geometry.get_coordinates()
geographic = bool(gdf.crs and gdf.crs.is_geographic)
decimals = (
coords["x"].astype(str).str.split(".").str[1].fillna("").str.len()
)
magnitude = float(np.abs(coords[["x", "y"]].to_numpy()).max())
unit_m = 111_320 if geographic else 1.0
implied_m = unit_m / (10 ** float(decimals.median()))
return {
"crs": str(gdf.crs),
"units": "degrees" if geographic else "linear",
"vertices": len(coords),
"max_magnitude": round(magnitude, 2),
"double_resolution_at_magnitude": float(np.spacing(magnitude)),
"median_decimal_places": int(decimals.median()),
"implied_ground_precision_m": round(implied_m, 6),
"storage_bytes_as_wkb": int(len(coords) * 16),
}
for k, v in precision_audit(gpd.read_file("data/raw/parcels.gpkg")).items():
print(f"{k:32} {v}")crs EPSG:27700
units linear
vertices 850821
max_magnitude 506878.39
double_resolution_at_magnitude 5.820766091346741e-11
median_decimal_places 2
implied_ground_precision_m 0.01
storage_bytes_as_wkb 13613136
Example 2: choose a grid size from the data's real accuracy
import geopandas as gpd
from shapely import set_precision
def grid_trial(gdf, grids=(0.0001, 0.001, 0.01, 0.1, 1.0)):
metric = gdf.to_crs(gdf.estimate_utm_crs())
base_area = metric.geometry.area.sum()
base_vertices = len(metric.geometry.get_coordinates())
for grid in grids:
snapped = set_precision(metric.geometry.values, grid_size=grid)
series = gpd.GeoSeries(snapped, crs=metric.crs)
valid = series[~series.is_empty & series.notna()]
print(
f"grid {grid:>7} m → "
f"area {(valid.area.sum() / base_area - 1):+.6%}, "
f"vertices {len(valid.get_coordinates()) / base_vertices:6.1%}, "
f"empty {len(series) - len(valid)}"
)
grid_trial(gpd.read_file("data/raw/parcels.gpkg"))grid 0.0001 m → area -0.000000%, vertices 100.0%, empty 0
grid 0.001 m → area +0.000001%, vertices 100.0%, empty 0
grid 0.01 m → area -0.000005%, vertices 99.9%, empty 4
grid 0.1 m → area -0.000074%, vertices 99.2%, empty 12
grid 1.0 m → area -0.000683%, vertices 91.3%, empty 50
Pick the coarsest grid whose area change is still negligible — that is the precision your data actually carries.
Example 3: tolerant assertions in tests
import pytest
from shapely import equals_exact
from shapely.geometry import Point
def test_reprojection_round_trips(parcels):
there_and_back = parcels.to_crs(4326).to_crs(parcels.crs)
for a, b in zip(parcels.geometry, there_and_back.geometry):
assert equals_exact(a, b, tolerance=2e-3) # 2 mm: this route moves ~1 mm per trip
def test_centroid_is_where_expected(square):
c = square.centroid
assert c.x == pytest.approx(50.0, abs=1e-9)
assert c.y == pytest.approx(50.0, abs=1e-9)
def test_area_is_stable_after_cleaning(raw, cleaned):
assert cleaned.geometry.area.sum() == pytest.approx(
raw.geometry.area.sum(), rel=1e-4) # 0.01% toleranceWith real inputs standing in for the fixtures:
from shapely.geometry import box
sample = gpd.read_file("data/raw/parcels.gpkg").head(500)
snapped = sample.set_geometry(set_precision(sample.geometry.values, grid_size=0.01))
for test, args in ((test_reprojection_round_trips, (sample,)),
(test_centroid_is_where_expected, (box(0, 0, 100, 100),)),
(test_area_is_stable_after_cleaning, (sample, snapped))):
try:
test(*args)
print(test.__name__, "passed")
except AssertionError as exc:
print(test.__name__, "FAILED")test_reprojection_round_trips passed
test_centroid_is_where_expected passed
test_area_is_stable_after_cleaning passed
Choosing the tolerance is part of the test, and the first version of this one chose wrong: 1e-6 — a micrometre — assumed a reprojection round trip only leaves floating-point residue, and on the Rutland parcels it failed, because this BNG–WGS 84 route moves every vertex about a millimetre. Measure what an operation actually does before you pin a tolerance to it: 1e-9 for a coordinate you computed yourself, what you measured for a reprojection, a relative tolerance for an area sum over many features.
Example 4: what accumulates through a chain of operations
import geopandas as gpd
from shapely.geometry import Polygon
from shapely import equals_exact
original = Polygon([(325000, 674000), (325100, 674000), (325100, 674100), (325000, 674100)])
series = gpd.GeoSeries([original], crs=27700)
drifted = series
for _ in range(10):
drifted = drifted.to_crs(4326).to_crs(27700) # ten reprojection round trips
a = original.exterior.coords[0]
b = drifted.geometry.iloc[0].exterior.coords[0]
print("start :", a)
print("after :", b)
print("moved :", ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** 0.5, "m")
print("equal :", equals_exact(original, drifted.geometry.iloc[0], tolerance=1e-6))start : (325000.0, 674000.0)
after : (325000.0078662433, 673999.9993289249)
moved : 0.00789481637126392 m
equal : False
Ten round trips moved the corner 7.9 mm — about the millimetre per trip the parcels showed, accumulating. That is not floating point; it is a transformation whose inverse is approximate, and it is exactly enough to break an exact comparison or to turn a shared boundary into a sliver once two copies drift differently.
What a 64-bit double can and cannot hold
Every coordinate in the GIS stack — Shapely, GEOS, GDAL, PostGIS, GeoPackage, shapefile — is an IEEE-754 double. That is 64 bits: one sign, eleven exponent, fifty-two fraction, giving roughly 15 to 17 significant decimal digits. Two properties follow, and almost every precision surprise in GIS comes from one of them.
The first is that most decimal fractions are not exactly representable in binary, so arithmetic carries a tiny residue. That residue is around 10⁻¹⁶ relative to the magnitude of the number, which is utterly irrelevant to any measurement and completely fatal to an equality test. Hence the rule that geometric comparisons take a tolerance: not because the data is imprecise, but because exact equality of computed floats is the wrong question.
The second is that precision is relative, not absolute. Near zero a double resolves 10⁻¹⁶; at 6,378,137 — the Earth's radius, and the magnitude of geocentric coordinates — it resolves about a nanometre. That is still far finer than any survey, which is why the "doubles are not precise enough for GIS" worry is misplaced. The real issue is accumulated error through chains of operations, and the fact that two copies of a boundary accumulate different errors.
This is where precision meets topology. Two adjacent parcels digitised from the same fence line should share a boundary exactly. After each has been reprojected, simplified and re-exported by different processes, their shared vertices differ in the last few bits — and GEOS, comparing exactly, correctly reports a sliver of overlap instead of a clean adjacency. Snapping both to a shared grid with set_precision() collapses those differences and restores the relationship, which is why it appears in every serious cleaning workflow.
The last piece is storage. Because precision is relative, writing 15 digits per ordinate stores about nine digits of noise for typical GIS data. Those digits cost real bytes — the GeoJSON above is 34% smaller at seven decimal places, with nothing lost that a survey recorded — and they actively cause the sliver problem by preserving differences that are not real. Deciding the precision your data genuinely carries, snapping to it, and recording that decision is one of the cheapest quality improvements available.
Equality tests and very large coordinates
==on geometries is not a coordinate test: It compares structure. Useequals_exact(tolerance=…)orgeom_equals_exactfor layers.- Very large coordinates lose absolute precision: Geocentric or Web Mercator coordinates near the poles have far coarser steps than local grid coordinates.
- Web Mercator distorts, it does not lose precision: A metre near the pole is a different number of Mercator units than at the equator — that is projection, not floats.
- Shapefile stores doubles too: Its precision limits come from the format's field definitions for attributes, not from geometry.
- GeoJSON precision is a writer option: GDAL writes 15 significant digits by default, and 7 decimal places (about 1 cm) only when asked, with
COORDINATE_PRECISION=7orRFC7946=YES:
import os
wgs84 = gpd.read_file("data/raw/parcels.gpkg").head(2000).to_crs(4326)
wgs84.to_file("default.geojson")
wgs84.to_file("seven.geojson", COORDINATE_PRECISION=7)
first = open("default.geojson").read().split('"coordinates": [ [ [ ')[1].split(" ]")[0]
print("default writes:", first)
print(f"default {os.path.getsize('default.geojson') / 1e6:.2f} MB, 7 dp {os.path.getsize('seven.geojson') / 1e6:.2f} MB")default writes: -0.591348602878956, 52.586531042083593
default 4.25 MB, 7 dp 2.79 MB
set_precisionre-nodes: It does not merely round; it rebuilds the geometry so snapped edges are genuinely shared. That is why it can create empty geometries from slivers.- Decimal or fixed-point types do not help: No mainstream GIS stack uses them, and they would not fix accumulated error in projections, which is trigonometric.
Internal links
- Topological Coverage: What "Clean" Means for Adjacent Polygons
- How to Snap and Align Geometries to Fix Slivers and Gaps in Python
- How to Reduce GIS File Size in Python Without Wrecking the Data
- WKT, WKB and GeoJSON: How Geometry Is Actually Stored
- Shapely TopologyException: Found Non-Noded Intersection (How to Fix)
- Coordinate Reference Systems (CRS) Explained for Python GIS
FAQ
Why does 0.1 + 0.2 != 0.3 matter in GIS?
Because geometric equality tests compare coordinates exactly. Any coordinate that has been through arithmetic carries a residue of about 1 part in 10¹⁶, so exact comparison fails on values that are, for every practical purpose, identical.
How many decimal places should I store?
Match your data's real accuracy. In degrees, 6 decimal places is about 0.1 m and covers almost everything; in metres, 2 or 3 places is centimetre to millimetre. More than that is stored noise.
Are doubles precise enough for GIS?
Yes, easily. Even at Earth-radius magnitudes a double resolves about a nanometre. Problems come from accumulated error and exact comparisons, not from the underlying precision.
What is the difference between precision and accuracy?
Precision is how many digits you store; accuracy is how close the value is to the truth. A badly surveyed point recorded to fifteen decimal places is precise and inaccurate.
How do I compare two geometries in a test?
shapely.equals_exact(a, b, tolerance=…), or geom_equals_exact for a whole GeoSeries. Pick the tolerance to match what the operation should have preserved.
Does snapping lose data?
It loses the digits below the grid size, which for a well-chosen grid are computational noise. Verify by checking that total area changes by a negligible fraction, and record the grid size you used.
Why does the same overlay work in one CRS and fail in another?
Coordinate magnitude changes the absolute size of floating-point steps, and different projections put the data at very different magnitudes. Working in a local projected CRS keeps numbers small and behaviour predictable.