#Geometry Encoding and Coordinate Precision

Open any national vector dataset and read a coordinate: -2.4783610000000002, 51.4544829999999987. Sixteen significant digits, describing a fence line that was walked with a GPS receiver accurate to about two metres. Fourteen of those digits are not measurement. They are the residue of a float64 round-trip through a projection library, and every one of them is stored, compressed, transferred, decompressed, and re-serialised on every read for the lifetime of the dataset.

This page is about the gap between the precision a file records and the accuracy the data possesses, and about closing that gap deliberately. It belongs to the Compression, Chunking & Spatial Indexing section, and it operates one layer below the tuning that section usually discusses: ZSTD level selection and row group sizing decide how well the compressor does its job; coordinate precision decides how much genuinely compressible data you hand it in the first place. Getting the encoding right routinely halves a geometry column, and it does so before the compressor runs, which means every downstream saving compounds.

#Prerequisites

• Python 3.10+ with shapely>=2.0, geopandas>=1.0, pyarrow>=16, and numpy>=1.26
• A dataset whose positional accuracy you can state — from survey metadata, sensor specification, or the imagery resolution it was digitised from. If nobody can answer this, that is the first problem to fix, not the encoding.
• A fixed CRS across the whole dataset, established before any of this begins; quantization grids are expressed in CRS units and mean nothing if the units vary
• A topology validation step you can run before and after, so the change is provably safe rather than plausibly safe

#Where the Bytes Actually Go

A geometry column in GeoParquet is a run of Well-Known Binary blobs, and a WKB polygon is mostly coordinates: a small header, a ring count, a point count, and then eight bytes per ordinate, forever. A 200-vertex polygon carries roughly 3,200 bytes of coordinate payload and about 20 bytes of everything else. Whatever you do to those eight-byte doubles is what happens to the file.

Compression does not see numbers; it sees byte patterns. A float64 is a sign bit, an 11-bit exponent, and a 52-bit mantissa. For coordinates within one country, the sign and exponent are nearly constant across the entire column and the high mantissa bits vary slowly — all of that compresses beautifully. The low mantissa bits, however, are close to uniformly random, and random bits are incompressible by construction. Roughly the last three bytes of every ordinate are pure entropy that no compression level can touch.

Three encoding levers act on this, and they are independent — you can apply one, two, or all three.

Quantization snaps every ordinate to a fixed grid, zeroing the noisy low bits. It is lossy by exactly the grid size and nothing more, it is applied uniformly across the dataset, and it is the lever with the best ratio of saving to risk. This is the one to reach for first.

Delta encoding stores each vertex as an offset from the previous one rather than as an absolute position. Consecutive vertices of a coastline differ by metres while the absolute values are hundreds of thousands of metres, so the deltas need far fewer significant bits. Parquet applies delta encoding to integer columns automatically, which is why quantized-to-integer geometry encodings compress so much better than float ones.

Narrowing the type replaces float64 with float32 or with a scaled integer. This is the biggest single win and the easiest way to silently destroy a dataset — a float32 has about seven significant decimal digits, and a longitude of -122.4194 has already consumed six of them before the fractional part begins.

#The Workflow

#1. Establish the accuracy the data actually has

Ask what produced the geometry. A cadastral survey delivers centimetres. A differential GPS track delivers decimetres. A boundary digitised from 10 m Sentinel-2 imagery delivers, at best, several metres, no matter how many decimals the shapefile contains. Write the answer into the dataset’s metadata, because every subsequent decision depends on it and it will otherwise be re-litigated forever.

#2. Convert accuracy into a grid in CRS units

A grid coarser than the source accuracy discards real signal; a grid much finer than it stores noise. Match them. In a projected CRS this is direct — 0.01 for centimetres in a metre-based system. In geographic coordinates it is a division: divide the target ground distance by 111,320 metres per degree of latitude, and for longitude divide additionally by the cosine of the working latitude, or simply use the latitude figure and accept slightly finer longitude resolution away from the equator.

#3. Snap on a shared absolute grid

Snap every geometry in the dataset against the same absolute grid origin, never relative to each geometry’s own bounding box. Shared boundaries survive because two coincident vertices round to the same grid point; they fracture the moment two layers use different origins or grid sizes.

#4. Re-validate topology

Snapping can collapse a ring whose extent is smaller than one grid cell, can turn a nearly-degenerate triangle into a line, and can introduce self-intersections where two vertices merge. All three are detectable. Validate before and after, compare the counts, and fail the job rather than shipping a dataset with new invalid geometries.

#5. Report size and error together

The output of this step is two numbers, not one: bytes saved and maximum positional displacement. A pipeline that reports only the first will eventually quantize a survey layer to metres and nobody will notice until a boundary dispute.

#Production Implementation

The function below quantizes a GeoDataFrame on a shared absolute grid, validates topology on both sides of the change, and returns the accuracy report the step is obliged to produce. It slots into the batch conversion pipeline immediately before the GeoParquet write.

# Requires: shapely>=2.0, geopandas>=1.0, numpy>=1.26  (Python 3.10+)
from __future__ import annotations

from dataclasses import dataclass

import geopandas as gpd
import numpy as np
from shapely import set_precision, get_coordinates
from shapely.errors import GEOSException

# One degree of latitude is ~111,320 m everywhere; longitude shrinks with
# cos(latitude), so using the latitude figure gives a grid at least as fine
# as requested in both axes.
METRES_PER_DEGREE = 111_320.0


@dataclass(frozen=True)
class PrecisionReport:
    grid: float
    features: int
    max_shift: float
    mean_shift: float
    invalid_before: int
    invalid_after: int
    collapsed: int


def grid_for_accuracy(accuracy_m: float, *, geographic: bool) -> float:
    """Translate a ground accuracy in metres into a grid in CRS units."""
    if accuracy_m <= 0:
        raise ValueError("accuracy must be positive")
    return accuracy_m / METRES_PER_DEGREE if geographic else accuracy_m


def quantize(
    frame: gpd.GeoDataFrame,
    accuracy_m: float,
    *,
    fail_on_new_invalid: bool = True,
) -> tuple[gpd.GeoDataFrame, PrecisionReport]:
    """Snap every ordinate to a shared absolute grid, then prove it was safe.

    Returns the quantized frame and a report carrying both the byte-side win
    (implied by the grid) and the accuracy cost (max/mean displacement).
    """
    if frame.crs is None:
        raise ValueError("a CRS is required before quantizing — grids are in CRS units")
    if frame.empty:
        raise ValueError("nothing to quantize")

    geographic = bool(frame.crs.is_geographic)
    grid = grid_for_accuracy(accuracy_m, geographic=geographic)

    invalid_before = int((~frame.geometry.is_valid).sum())
    before = frame.geometry.to_numpy()

    try:
        # set_precision snaps against an absolute grid anchored at the CRS
        # origin — identical for every geometry, so shared edges stay shared.
        snapped = set_precision(before, grid_size=grid)
    except GEOSException as exc:
        raise RuntimeError(f"precision reduction failed on grid {grid}: {exc}") from exc

    out = frame.copy()
    out.geometry = snapped

    collapsed = int(sum(1 for g in snapped if g is None or g.is_empty))
    invalid_after = int((~out.geometry.is_valid).sum())
    if fail_on_new_invalid and invalid_after > invalid_before:
        raise RuntimeError(
            f"grid {grid} introduced {invalid_after - invalid_before} invalid "
            f"geometries — use a finer grid or repair upstream"
        )

    shifts = []
    for old, new in zip(before, snapped):
        if old is None or new is None or new.is_empty:
            continue
        a, b = get_coordinates(old), get_coordinates(new)
        if a.shape != b.shape:      # a vertex was removed by the snap
            shifts.append(grid)
            continue
        shifts.append(float(np.abs(a - b).max()))

    scale = METRES_PER_DEGREE if geographic else 1.0
    arr = np.asarray(shifts) * scale if shifts else np.zeros(1)

    return out, PrecisionReport(
        grid=grid,
        features=len(frame),
        max_shift=float(arr.max()),
        mean_shift=float(arr.mean()),
        invalid_before=invalid_before,
        invalid_after=invalid_after,
        collapsed=collapsed,
    )

For the narrowing lever, the safe form is a scaled integer rather than a float32 — it makes the precision explicit, it is exactly representable, and Parquet’s delta encoding acts on it directly:

# Requires: numpy>=1.26, pyarrow>=16 — coordinates as scaled int32 offsets
import numpy as np


def to_scaled_offsets(coords: np.ndarray, origin: tuple[float, float], scale: int):
    """Encode ordinates as int32 offsets from a dataset origin.

    scale = 100 stores centimetres in a metre-based CRS; the int32 range then
    covers +/- 21,474 km from the origin, which is the whole planet.
    """
    if coords.ndim != 2 or coords.shape[1] != 2:
        raise ValueError("expected an (n, 2) coordinate array")
    shifted = (coords - np.asarray(origin, dtype=np.float64)) * scale
    rounded = np.rint(shifted)
    if np.abs(rounded).max() > np.iinfo(np.int32).max:
        raise ValueError("origin/scale combination overflows int32 — re-centre the origin")
    return rounded.astype(np.int32)

#Reference Matrix

Measured on a 4.2 million-feature national parcel layer in EPSG:27700 (metres), written as GeoParquet with ZSTD level 3 and 128 MB row groups. Sizes are the geometry column only.

The float32 row is the trap. It looks competitive with the 10 cm grid on size while delivering worse and — crucially — non-uniform accuracy: the error depends on the magnitude of the coordinate, so it is small near the projection origin and grows with distance from it. A dataset that validates fine in one region fails silently in another. The int32 row achieves better compression with bounded, uniform, documented error, which is why it is the right narrowing choice when narrowing is warranted at all.

#Failure Modes and Gotchas

Quantizing after re-projection instead of before. Re-projection is a non-linear transform, so a grid that was 10 cm in the source CRS becomes an irregular, latitude-dependent spacing in the target. Snap in the CRS you intend to store, as the final step before writing, and never between two projection operations.

Different grids for layers that must align. Parcels snapped to 1 cm and the roads that bound them snapped to 1 m will not share edges, and every downstream overlay will produce slivers. Grid size is a dataset-family decision, not a per-file one; record it in the metadata the way the metadata preservation workflow records CRS.

Collapsing features smaller than the grid. A 60 cm utility access chamber quantized to a 1 m grid becomes a point or vanishes. The implementation above counts collapses precisely so this is caught in the report rather than discovered in a field survey. Check the count, do not just log it.

Assuming the compressor will do it for you. ZSTD at level 19 cannot compress random bits any better than level 3 can — that is what random means. Teams that respond to a large geometry column by raising the compression level spend a great deal of CPU for a few percent, when quantization was available for a fraction of the cost. Fix the entropy first, then tune the compression level.

Treating precision as reversible. It is not. Once the low bits are gone they are gone, so keep an unmodified archival copy of survey-grade source data and quantize only the delivery copies. Storage for one master is cheap; a re-survey is not.

#Frequently Asked Questions

How many decimal places of latitude and longitude do I actually need?

At the equator, one degree of longitude is about 111,320 metres, so the sixth decimal place is roughly 11 centimetres and the seventh is about 1 centimetre. Cadastral and survey data justifies seven; utility and asset inventories are well served by six; administrative boundaries, land cover, and most analytical layers are fully served by five, which is about 1.1 metres. Digitised-from-imagery datasets rarely have real accuracy past five, so storing eight or nine places records digitisation noise at full cost.

Why does rounding coordinates make a compressed file smaller when the values are the same width?

A float64 always occupies eight bytes, so rounding does not shrink the raw column. It shrinks the compressed column, because compression exploits repetition and predictability. Full-precision coordinates have effectively random low-order mantissa bits, which are incompressible; snapping to a grid makes those bits constant across the whole column, so the compressor encodes them once. The saving lands entirely in the entropy coding stage, which is why it shows up only after ZSTD or Snappy runs.

Is float32 safe for storing geographic coordinates?

Not for global data in degrees. A float32 carries about 7 significant decimal digits, and a longitude near 100 degrees consumes three of them before the decimal point, leaving roughly 1 metre of resolution at best and worse near the extremes. It is usable for projected coordinates in a local CRS with a small numeric range, or for data whose accuracy is genuinely tens of metres. For anything that will be re-projected, joined, or used for measurement, keep float64 and control size through quantization instead.

Does coordinate quantization break shared boundaries between polygons?

Only if you snap each polygon independently with different grid origins, which you should never do. When every geometry in the dataset is snapped to the same absolute grid, two polygons that shared a vertex still share it afterwards, because both round to the same grid point. Slivers appear when the grid is coarser than the smallest real feature, when different layers use different grids, or when snapping happens after a re-projection rather than before.

#Related

• Coordinate Precision Reduction for Smaller GeoParquet — the quantization step end to end, with the validation harness
• WKB vs Native Arrow Geometry Encoding — the container choice that decides whether delta encoding can act at all
• Simplifying Geometries Before Compression — removing vertices rather than digits, and when each is right
• ZSTD Compression Levels for Geospatial Data — the stage that runs after encoding, and why it cannot fix entropy
• Dictionary Encoding for Categorical GIS Attributes — the same argument applied to the attribute columns

← Back to Compression, Chunking & Spatial Indexing