Transforming GeoPackage Layers Without Losing Precision

A GeoPackage stores geometry as float64 and its CRS as WKT2, which makes it the one mainstream container that can hold a cadastral coordinate and its identity intact — and it is still possible to lose both by transforming carelessly. This guide, part of geospatial file formats and CRS metadata, transforms a layer in place-preserving fashion: the geometry through a pinned operation at full precision, the CRS entry replaced rather than mutated, and the epoch and operation recorded where a reader will find them.

Where Precision Actually Escapes #

Three places, none of them the arithmetic.

Figure — three leaks, none of them in the arithmetic.

The read. A reader that materialises geometry through a text representation — well-known text, or a coordinate string — truncates at whatever precision that representation uses. Reading well-known binary keeps float64 all the way.

The write-back. Writing transformed coordinates into a layer whose geometry column was created for a different CRS leaves the spatial reference entry describing the old one, so the file says one thing and contains another. The CRS entry has to be written as part of the same operation.

The attribute copy. Cadastral layers often carry easting and northing as attributes as well as geometry, for convenience. Those columns are frequently REAL with an application-level format applied, and they are the ones that quietly end up at two decimal places while the geometry beside them is exact.

Complete Runnable Implementation #

from __future__ import annotations

import sqlite3
from dataclasses import dataclass
from pathlib import Path

import numpy as np
from pyproj import CRS, Transformer


@dataclass(frozen=True)
class LayerTransform:
    """Transform one GeoPackage layer, preserving precision and identity."""

    gpkg: Path
    layer: str
    source_crs: CRS
    target_crs: CRS
    pipeline: str | None = None       # pin the operation where it matters
    epoch: float | None = None

    def transformer(self) -> Transformer:
        if self.pipeline:
            return Transformer.from_pipeline(self.pipeline)
        return Transformer.from_crs(self.source_crs, self.target_crs, always_xy=True)

    def register_target_crs(self, conn: sqlite3.Connection) -> int:
        """Insert the target CRS as its own srs_id; never mutate an existing row.

        Editing the definition of an srs_id already referenced by other layers
        silently relabels them. A new row keeps the two layers honest even when
        both live in the same file.
        """
        auth = self.target_crs.to_authority()
        if auth is None:
            raise ValueError("target CRS has no authority code; register it explicitly")
        org, code = auth
        srs_id = int(code)
        row = conn.execute(
            "SELECT definition FROM gpkg_spatial_ref_sys WHERE srs_id = ?", (srs_id,)
        ).fetchone()
        wkt2 = self.target_crs.to_wkt(version="WKT2_2019")
        if row is None:
            conn.execute(
                "INSERT INTO gpkg_spatial_ref_sys "
                "(srs_name, srs_id, organization, organization_coordsys_id, "
                " definition, description) VALUES (?, ?, ?, ?, ?, ?)",
                (self.target_crs.name, srs_id, org, srs_id, wkt2,
                 f"written by the transformation pipeline; epoch {self.epoch}"),
            )
        elif row[0] != wkt2:
            raise ValueError(
                f"srs_id {srs_id} already exists with a different definition; "
                f"resolve this rather than overwriting it"
            )
        return srs_id

    def transform_points(self, xy: np.ndarray) -> np.ndarray:
        """Transform an (n, 2) float64 array. Non-finite rows are masked, not dropped."""
        if xy.dtype != np.float64 or xy.ndim != 2 or xy.shape[1] != 2:
            raise ValueError("expected an (n, 2) float64 array")
        good = np.isfinite(xy).all(axis=1)
        out = np.full_like(xy, np.nan)
        if good.any():
            e, n = self.transformer().transform(xy[good, 0], xy[good, 1])
            out[good, 0] = e
            out[good, 1] = n
        if not np.isfinite(out[good]).all():
            raise ValueError("transformation produced non-finite output for finite input")
        return out

Parameter Reference #

Worked Example #

import numpy as np
from pathlib import Path
from pyproj import CRS

job = LayerTransform(
    gpkg=Path("parcels.gpkg"),
    layer="boundaries",
    source_crs=CRS.from_epsg(4267),          # NAD27
    target_crs=CRS.from_epsg(6318),          # NAD83(2011)
    pipeline="+proj=pipeline +step +proj=hgridshift +grids=us_noaa_conus.tif",
    epoch=2010.00,
)
xy = np.array([[-122.600000, 45.500000],
               [-122.601000, 45.501000]], dtype=np.float64)
print(np.round(job.transform_points(xy), 8))
# [[-122.60117800   45.50005900]
#  [-122.60217800   45.50105900]]

Pinning the pipeline is the difference between a layer whose provenance can be stated and one whose operation was chosen by whatever PROJ release ran, which is the argument made at length in PROJ pipeline strings vs the pyproj Transformer API.

Figure — the order that keeps geometry and CRS consistent at every moment.

Validation Check #

def assert_geometry_and_srs_agree(conn: sqlite3.Connection, layer: str,
                                  expected_srs_id: int) -> None:
    """The layer's declared srs_id must match what the geometry now contains."""
    row = conn.execute(
        "SELECT srs_id FROM gpkg_geometry_columns WHERE table_name = ?", (layer,)
    ).fetchone()
    assert row is not None, f"layer {layer} has no geometry column entry"
    assert row[0] == expected_srs_id, (
        f"layer {layer} declares srs_id {row[0]} but was written in "
        f"{expected_srs_id}; the file now contradicts itself"
    )

A file that contradicts itself is worse than one with no CRS at all: the reader has no reason to doubt it, and the coordinates are wrong by whatever the two systems differ by.

Common Mistakes #

Mutating an existing gpkg_spatial_ref_sys row. Convenient, and it relabels every other layer referencing that srs_id without touching their geometry. Insert a new row, and if the id is taken with a different definition, stop rather than reconcile silently.

Transforming through well-known text. Round-tripping geometry through a text representation to get at the coordinates truncates them at the representation’s precision. Read and write well-known binary, keep float64 throughout, and only format to text at the point of export.

Attribute columns left stale. A layer carrying easting and northing attributes alongside its geometry now has two coordinate systems in one row: transformed geometry, untransformed attributes. Either update them in the same pass or drop them — a stale attribute is more dangerous than an absent one, because it looks authoritative.

What to Check After the Transformation #

A transformed layer looks the same as an untransformed one, so the confidence has to come from checks rather than from inspection. Four are worth running before the file is handed on.

Figure — four post-transformation checks, and the specific failure each catches.

Round-trip a sample. Transform a few hundred points back to the source CRS and compare against the originals; closure to a fraction of a millimetre confirms the operation ran in the direction intended. A systematic offset means the pipeline was applied inverted, which is easy to do and invisible in the output.

Compare against control. Where the area contains published control, transform those positions and compare — this is the only check that says the coordinates are right rather than merely self-consistent, and it is what the deliverable’s residual statistics are built from anyway.

Confirm the extent moved as expected. A datum shift moves a layer by decimetres; a projection change moves it by hundreds of kilometres. Comparing the bounding box before and after against the expected magnitude catches a wrong operation instantly and costs nothing.

Verify the CRS entry and the geometry agree. The assertion given earlier does this, and it is the check that catches the specific failure mode of an in-place edit: geometry updated, spatial reference row left describing the old system.

Running all four takes seconds on a sample and turns “the transformation completed without errors” into “the transformation produced what it was supposed to produce”, which are very different claims.

Frequently Asked Questions #

Should I transform in place or write a new file? #

Write a new file. In-place transformation makes the previous state unrecoverable, and it leaves a window in which the geometry has been updated and the CRS entry has not. A new file also gives the deliverable a natural identity for the manifest and the audit hash.

Does a GeoPackage carry the coordinate epoch? #

Through its WKT2 definition it can, using the coordinate-metadata construct, and support varies by reader. Because it varies, record the epoch in the sidecar as well — duplication is cheap and the failure mode of an epoch that only exists in a field nothing reads is expensive.

How do I handle a layer with mixed geometry types? #

Transform the coordinates, not the geometry model: extract every coordinate pair, transform them as one array, and write them back into the same structure. The transformation is per point and indifferent to whether the point belongs to a line or a polygon ring, and going through a per-feature loop instead is where both time and precision get lost.

What about the spatial index? #

Rebuild it after transforming. A GeoPackage R-tree index stores bounding boxes in the old coordinate system, so leaving it in place produces a file whose spatial queries return the wrong features — correct geometry, wrong index, and no error anywhere.

Can several layers in one file use different CRSs? #

They can — the format allows a per-layer spatial reference — and it is worth avoiding in a deliverable. A recipient who opens one layer and infers the file’s CRS from it will be wrong about the others, and geometry operations across layers will silently mix systems. Where mixed CRSs are unavoidable, say so in the manifest and name the CRS per layer.

Related #

• Geospatial file formats and CRS metadata — the parent topic
• Reading and writing .prj files with Python
• Vectorizing coordinate transforms with NumPy and Dask
• PROJ pipeline strings vs the pyproj Transformer API
• Validating coordinate precision to millimeter standards