Coordinate Precision and Floating Point in GIS Explained

Problem statement

Two things that should be equal are not:

>>> from shapely.geometry import Point
>>> Point(0.1 + 0.2, 1).equals(Point(0.3, 1))
False
>>> 0.1 + 0.2
0.30000000000000004

And two boundaries that should coincide, do not:

>>> a.exterior.coords[2]
(325104.20000000001, 673992.69999999995)
>>> b.exterior.coords[0]
(325104.19999999995, 673992.70000000007)

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 60% 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 point

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

Bar chart of representable precision at different coordinate magnitudes for a 64-bit double. — The same 64 bits buy nanometres near the origin and millimetres at a projected-grid magnitude.

Step-by-step solution

Vertical steps: know the data accuracy, pick a grid, snap, verify area change, store. — Five steps that turn a vague worry about precision into a chosen, documented number.

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.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")

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):,}")

A layer whose coordinates all carry 12 decimal places did not come from a survey with picometre accuracy; it came from a reprojection. The extra digits are computational residue, and 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%}")

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
same = gdf_a.geometry.geom_equals_exact(gdf_b.geometry, tolerance=1e-6)
print(f"{same.sum()} of {len(same)} geometries match within 1e-6")

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 = prepare(gpd.read_file("data/raw/parcels.gpkg"))
zones = prepare(gpd.read_file("data/ref/zones.gpkg").to_crs(parcels.crs))

overlay = gpd.overlay(parcels, zones, how="intersection")
print(len(overlay), "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.

Know when precision is not the problem

import geopandas as gpd

gdf = gpd.read_file("data/raw/sites.gpkg")
print("CRS:", gdf.crs)
print("bounds:", [round(v, 2) for v in gdf.total_bounds])

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.

Code examples

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}")

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"))

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=1e-6)      # 1 µm in a metric CRS

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% tolerance

Choosing the tolerance is part of the test: 1e-9 for a coordinate you computed yourself, 1e-6 after 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))

The drift is typically sub-millimetre — harmless on its own, and exactly enough to turn a shared boundary into a sliver once two copies drift differently.

Explanation

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.

Panels contrasting precision (digits stored) with accuracy (agreement with reality). — Fifteen decimal places of a badly surveyed point is precision without accuracy.

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 an overlap of 0.003 m² 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 — a WKT file can be a third larger for nothing — 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.

Edge cases or notes

== on geometries is not a coordinate test: It compares structure. Use equals_exact(tolerance=…) or geom_equals_exact for 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 writers often truncate: GDAL defaults to 7 decimal places (COORDINATE_PRECISION), which is about 1 cm — usually fine, occasionally not.
set_precision re-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

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.

Coordinate Precision and Floating Point in GIS Explained #

Problem statement #

Quick answer #

What a double can hold #

Step-by-step solution #

Why 0.1 + 0.2 is not 0.3 #

Precision depends on magnitude #

What decimal places mean on the ground #

Match stored precision to real accuracy #

Comparing geometries safely #

Snap before operations that must agree #

Know when precision is not the problem #

Code examples #

Example 1: a precision audit for a layer #

Example 2: choose a grid size from the data's real accuracy #

Example 3: tolerant assertions in tests #

Example 4: what accumulates through a chain of operations #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

Why does 0.1 + 0.2 != 0.3 matter in GIS? #

How many decimal places should I store? #

Are doubles precise enough for GIS? #

What is the difference between precision and accuracy? #

How do I compare two geometries in a test? #

Does snapping lose data? #

Why does the same overlay work in one CRS and fail in another? #