Datum Transformations & Projection Accuracy in Spatial SQLite

Reprojection looks like arithmetic and behaves like a lookup. Convert a national-grid coordinate to WGS 84 on a workstation and again inside a field application, and the two answers can differ by two metres — not because either is wrong, but because the two machines had different transformation data available and each silently picked the best option it could find. Nothing warns you. The coordinates are plausible, the layer renders, and the discrepancy only surfaces when someone stands where the map says the boundary is.

This guide is part of the Core Architecture & Format Standards for Spatial SQLite section. It explains where that variability comes from, how to make a transformation deterministic, and how to verify a reprojected layer against something other than hope.

Prerequisites #

• PROJ 8.0+ — the transformation-pipeline model and the projinfo tool both assume it
• pyproj 3.4+ — exposes the transformer selection and accuracy metadata this guide relies on
• A working proj.db reachable through PROJ_LIB or the library’s compiled-in default
• Two or more known control points whose coordinates you have in both reference systems
• Familiarity with the registry tables covered in SpatiaLite Metadata Tables Explained

Concept & Specification Reference #

Three distinct things get called “the projection”, and separating them is most of the battle.

A projection error is loud: the numbers change magnitude, the layer lands in the ocean, someone notices within a minute. A datum error is quiet, because the wrong answer is only a few metres from the right one — comfortably inside the error people attribute to GPS.

The critical point is that a transformation is not implied by the two reference systems. Between a national datum and WGS 84 there are usually several published routes: a three-parameter shift that is simple and approximate, a seven-parameter one that is better, and a grid-based one that is best and requires a data file most installations do not have. PROJ picks the most accurate route whose data it can actually find, which means the route depends on the machine.

Step-by-Step Implementation #

1. Ask which route PROJ is going to take #

Before transforming anything, find out what the library intends to do. pyproj exposes the candidate operations and their published accuracy, which turns an invisible decision into a visible one.

# List the transformation routes PROJ knows between two reference systems
from pyproj.transformer import TransformerGroup

tg = TransformerGroup("EPSG:27700", "EPSG:4326")

for t in tg.transformers:
    print(f"{t.description}\n    accuracy: {t.accuracy} m")

# Routes PROJ knows about but cannot run, because a grid is missing
for missing in tg.unavailable_operations:
    print("UNAVAILABLE:", missing.name)
    for grid in missing.grids:
        print("    needs grid:", grid.short_name, "->", grid.url)

unavailable_operations is the field that matters. A machine that lists a grid-based route there is a machine that will silently fall back to something less accurate, and the URL it prints is exactly what needs installing.

2. Pin the transformation rather than the endpoints #

Once you know which route you want, name it. pyproj accepts an explicit pipeline, which removes the selection step entirely and makes the result reproducible regardless of what else the machine has installed.

# Deterministic transformation: state the route, not just the endpoints
from pyproj import Transformer

# Selection by authority code — the operation is now fixed, not chosen
transformer = Transformer.from_pipeline(
    "+proj=pipeline "
    "+step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 "
    "+x_0=400000 +y_0=-100000 +ellps=airy "
    "+step +proj=hgridshift +grids=uk_os_OSTN15_NTv2_OSGBtoETRS.tif "
    "+step +proj=unitconvert +xy_in=rad +xy_out=deg"
)

lon, lat = transformer.transform(651409.903, 313177.270)
print(f"{lon:.9f}, {lat:.9f}")

The pipeline is verbose, and that is the point: every step is stated, so two machines running it either produce identical output or fail loudly because a grid is absent. A Transformer.from_crs call produces neither guarantee.

3. Record the transformation alongside the data #

A container that has been reprojected should say how. The registry records the reference system, not the route taken to reach it, so the route belongs in your own metadata — a small table in the container itself is the version least likely to be separated from the data.

-- A provenance table travels with the container
CREATE TABLE IF NOT EXISTS _crs_provenance (
    layer_name   TEXT NOT NULL,
    source_srs   INTEGER NOT NULL,
    target_srs   INTEGER NOT NULL,
    operation    TEXT NOT NULL,      -- the authority code or pipeline string
    accuracy_m   REAL,
    proj_version TEXT NOT NULL,
    applied_at   TEXT NOT NULL,
    PRIMARY KEY (layer_name, applied_at)
);

Recording the PROJ version matters as much as the operation. Published transformations are revised, and a pipeline that produced one answer under PROJ 8 can produce a slightly different one under PROJ 9 — a difference that is impossible to explain later without knowing which version ran.

4. Bundle the grids you depend on #

A grid-based transformation is only reproducible if the grid ships with the application. PROJ resolves grid files through its data directory, so bundling means placing the files where the library looks and confirming that it found them.

# Confirm the grid the pipeline needs is actually resolvable
import pyproj
from pathlib import Path

print("PROJ data dirs:", pyproj.datadir.get_data_dir())

grid = Path(pyproj.datadir.get_data_dir()) / "uk_os_OSTN15_NTv2_OSGBtoETRS.tif"
if not grid.exists():
    raise RuntimeError(
        f"grid missing at {grid} — transformations will fall back silently"
    )

Raising here is deliberate. The failure this check exists to prevent is not an exception — it is a successful run producing coordinates two metres out, and the only way to convert that into something noticeable is to refuse to start.

5. Reproject the layer, then re-register it #

Reprojecting geometry and updating the registry are two operations, and doing only the first produces a layer whose coordinates and declared reference system disagree — the same class of error as calling set_crs where to_crs was meant.

# Reproject a layer in place and update its registered SRS together
import sqlite3
from shapely import from_wkb, to_wkb
from shapely.ops import transform

conn = sqlite3.connect("survey.gpkg", isolation_level=None)
conn.enable_load_extension(True)
conn.load_extension("mod_spatialite")

rows = conn.execute("SELECT fid, geom FROM parcels").fetchall()
reprojected = [
    (to_wkb(transform(transformer.transform, from_wkb(blob))), fid)
    for fid, blob in rows if blob is not None
]

conn.execute("BEGIN IMMEDIATE")
conn.executemany(
    "UPDATE parcels SET geom = GeomFromWKB(?, 4326) WHERE fid = ?", reprojected
)
conn.execute("UPDATE gpkg_contents SET srs_id = 4326 WHERE table_name = 'parcels'")
conn.execute(
    "UPDATE gpkg_geometry_columns SET srs_id = 4326 WHERE table_name = 'parcels'"
)
conn.execute("COMMIT")

Validation & Verification #

Control points are the only real verification. Take coordinates you know in both systems — a survey monument, a benchmark, a corner whose position is published — transform them, and measure the residual.

# Verify the pipeline against known control points
CONTROLS = [
    # (easting, northing, expected_lon, expected_lat)
    (651409.903, 313177.270, 1.716073972, 52.658007833),
    (438700.000, 114800.000, -0.939621000, 50.868331000),
]

for e, n, elon, elat in CONTROLS:
    lon, lat = transformer.transform(e, n)
    # ~111 km per degree of latitude; adequate for a residual sanity check
    dy = (lat - elat) * 111_320
    dx = (lon - elon) * 111_320 * abs(lat) ** 0
    print(f"residual: {dx:+.3f} m east, {dy:+.3f} m north")
    assert abs(dx) < 0.2 and abs(dy) < 0.2, "transformation route is not the expected one"

A residual of a few centimetres means the grid-based route ran. A residual of one to two metres means it did not, and the assertion is what turns that from a fact nobody noticed into a failed build.

Cross-check the registry too, because a correct transformation written into a mis-registered layer is still wrong:

-- GeoPackage: the two registry rows must agree with each other
SELECT c.table_name, c.srs_id AS contents_srs, g.srs_id AS geometry_srs
FROM gpkg_contents c
JOIN gpkg_geometry_columns g ON g.table_name = c.table_name
WHERE c.data_type = 'features' AND c.srs_id <> g.srs_id;

Common Failure Modes & Fixes #

Coordinates land in the Gulf of Guinea #

Diagnosis: Longitude and latitude were swapped, or a projected coordinate was interpreted as degrees. Both produce values near zero, and zero-zero is off the coast of West Africa. Fix: Check axis order. EPSG:4326 is formally latitude-first, and different libraries disagree about whether to honour that; pyproj respects it unless the transformer is created with always_xy=True, which is almost always what code assuming longitude-first expects.

The same file reprojects differently on two machines #

Diagnosis: Different transformation routes, because one machine has a grid the other does not. Fix: Pin the pipeline explicitly and assert the grid is present at startup. Comparing TransformerGroup(...).unavailable_operations on both machines identifies the difference in one command.

ST_Transform returns NULL for every row #

Diagnosis: The target SRID has no definition in the container’s reference-system registry, or PROJ cannot open proj.db. SpatiaLite returns NULL rather than raising in both cases. Fix: Insert the missing definition as described in Managing Spatial Reference Systems in SQLite, and assert proj_version() returns a value during startup.

Everything shifts by a consistent amount after a library upgrade #

Diagnosis: The published transformation was revised, or the new PROJ version prefers a different route. Fix: This is not a bug to work around — it is usually a correction. Compare the two routes with projinfo, decide which one your published data should use, and pin it. Recording the PROJ version in the provenance table is what makes this diagnosable at all.

Areas computed after reprojection are wrong #

Diagnosis: Area was computed in a geographic reference system, where the units are degrees and the result is meaningless. Fix: Compute area in a projected system appropriate to the extent, or use a geodesic area function. Reprojecting to EPSG:4326 for storage and interchange is normal; computing measurements there is not.

Performance Notes #

Transformation cost is dominated by two things, and neither is the arithmetic. The first is transformer construction: resolving a route involves database lookups and, for grid-based routes, opening and indexing a grid file. Construct the transformer once and reuse it — a loop that calls Transformer.from_crs per feature can spend the overwhelming majority of its time on setup.

The second is vertex count. pyproj transforms arrays far faster than it transforms scalars, so passing whole coordinate arrays rather than iterating point by point makes a large difference on polygon layers. shapely.ops.transform is convenient and works per geometry; where a layer is large enough for the difference to matter, extracting coordinates into arrays, transforming in bulk, and rebuilding geometry is measurably faster.

Finally, reprojection is a full rewrite of every geometry, which means it invalidates every bounding box in the spatial index. Budget the index rebuild as part of the reprojection rather than as an afterthought — on a large layer it commonly costs more than the transformation did.

Child Pages #

Pages in this section go deeper on individual reference-system tasks:

• How to Add a Custom CRS to gpkg_spatial_ref_sys — inserting a definition a container does not ship with
• Choosing Between EPSG:4326 and a Projected CRS — storage, measurement and interchange pull in different directions
• Datum Shifts and Why Coordinates Move by Metres — the mechanism behind the quiet error
• How to Bundle PROJ Grid Files with an Application — shipping the data a deterministic pipeline needs
• Verifying a Reprojection with Known Control Points — the residual check as a build gate

Frequently Asked Questions #

Related #

• Core Architecture & Format Standards for Spatial SQLite — parent section: format contracts, registries and file structure
• SpatiaLite Metadata Tables Explained — where reference-system definitions live and how layers bind to them
• Managing Spatial Reference Systems in SQLite — inserting and binding definitions correctly
• Geometry Validity & Topology Repair — the check a reprojected layer owes itself