CoordTransform Cadastral & High-Precision Python

Understanding NTv2 Grid Shift Files in Python

Parsing an NTv2 .gsb grid is the binary-I/O sub-task that the rest of Core Transformation Fundamentals & Standards depends on: every NTv2-based datum shift you run — and every audit you have to defend — begins by reading this fixed-layout file correctly, in the right byte order, with the shift arrays promoted to a precision that will not leak millimetres before interpolation even starts. This page narrows the standards framework to one concrete operation: opening a National Transformation version 2 grid shift file, decoding its master and sub-grid headers, and loading the latitude/longitude shift surfaces into memory as deterministic, validatable NumPy arrays. Unlike a parametric Helmert shift that carries seven constants, NTv2 stores an empirically derived correction at every node of a dense geographic grid, which is precisely why it can model the non-uniform distortion between a legacy survey network and a modern geocentric frame — and precisely why a single mis-aligned seek silently corrupts a cadastral deliverable. The chosen method that determines whether you reach for an NTv2 grid at all is covered in NADCON vs NTv2: choosing the right datum shift; here we assume NTv2 has been selected and focus on reading it faithfully.

The NTv2 .gsb grid hierarchy A single master overview header (NUM_OREC, NUM_SREC, NUM_FILE and the source and target datums) sits above a sequence of sub-grids numbered 0 through n. Each sub-grid carries its own extent, node spacing and shift arrays. Sub-grids nest: sub-grid 1 is a child that refines sub-grid 0, and sub-grid n refines further. Sub-grid 0 expands to show that every node is a 16-byte record of four float32 values: latitude shift, longitude shift, latitude accuracy and longitude accuracy, expressed in arc-seconds with longitude positive west. Master header (176 bytes) NUM_OREC · NUM_SREC · NUM_FILE · source/target datums Sub-grid 0 extent · spacing · shifts Sub-grid 1 nested child grid Sub-grid n finest refinement nests refines Per node — 16 bytes latitude shift · longitude shift latitude accuracy · longitude accuracy 4 × float32 · arc-seconds · lon +W

Figure — the NTv2 .gsb hierarchy: one master header over nested sub-grids of shift nodes.

Binary Specification: the .gsb Layout

The .gsb format is a record-oriented binary structure. Both the master (overview) header and every sub-grid header are exactly 11 records of 16 bytes, where each record is an 8-byte ASCII keyword followed by an 8-byte value — giving a fixed 176-byte header block. After each sub-grid header come its shift nodes, each a 16-byte record of four IEEE 754 single-precision floats: latitude shift, longitude shift, latitude accuracy, longitude accuracy.

Three properties of the specification dominate any correct reader:

  • Byte order is not fixed. The original Canadian grids are big-endian, while many grids redistributed through open toolchains are little-endian. The reliable discriminator is NUM_OREC, the first numeric value, which must decode to 11. Whichever endianness yields 11 is the file’s byte order; guessing wrong turns every double into garbage.
  • Angular units are arc-seconds, longitude positive west. Sub-grid extents (S_LAT, N_LAT, E_LONG, W_LONG) and the node shifts themselves are expressed in seconds of arc, and NTv2 follows the geodetic convention that longitude increases westward. Converting a shift to metres on the ground therefore requires the local meridian and prime-vertical radii — see the arc-second-to-metre reduction in projection math fundamentals for cadastral surveys.
  • -999.0 is the null sentinel. A node carrying the null flag is unmodelled; interpolating across it instead of rejecting the query is a classic source of out-of-extent error.

The sub-grid header field offsets, measured from the start of the 176-byte block, are fixed by the specification:

Annotated 176-byte NTv2 sub-grid header The sub-grid header is 11 records of 16 bytes, totalling 176 bytes. Each record is an 8-byte ASCII keyword in the left half followed by its 8-byte value in the right half, so every value begins 8 bytes after its keyword. Reading top to bottom: SUB_NAME value at byte offset 8, PARENT at 24, then CREATED and UPDATED metadata, S_LAT at 72, N_LAT at 88, E_LONG at 104, W_LONG at 120, LAT_INC at 136, LONG_INC at 152, and GS_COUNT at 168. Offsets are measured from the start of the 176-byte block. Sub-grid header — 11 records × 16 bytes = 176 bytes 8-byte keyword 8-byte value 0 SUB_NAME sub-grid name @8 16 PARENT parent / NONE @24 32 CREATED date metadata 48 UPDATED date metadata 64 S_LAT south limit (sec) @72 80 N_LAT north limit (sec) @88 96 E_LONG east limit (sec, +W) @104 112 W_LONG west limit (sec, +W) @120 128 LAT_INC lat spacing (sec) @136 144 LONG_INC lon spacing (sec) @152 160 GS_COUNT rows × cols (int32) @168 Each value begins 8 bytes after its keyword · offsets from the start of the block
Record Keyword Value offset Type Meaning
0 SUB_NAME 8 char[8] sub-grid name
1 PARENT 24 char[8] parent name (NONE at top level)
4 S_LAT 72 float64 south latitude limit (arc-sec)
5 N_LAT 88 float64 north latitude limit (arc-sec)
6 E_LONG 104 float64 east longitude limit (arc-sec, +W)
7 W_LONG 120 float64 west longitude limit (arc-sec, +W)
8 LAT_INC 136 float64 latitude node spacing (arc-sec)
9 LONG_INC 152 float64 longitude node spacing (arc-sec)
10 GS_COUNT 168 int32 node count = rows × cols

The node grid dimensions follow directly from the extents and spacing, with one extra node because the grid is inclusive of both edges:

rows=N_LAT−S_LATLAT_INC+1,cols=W_LONG−E_LONGLONG_INC+1\text{rows} = \frac{N\_LAT - S\_LAT}{LAT\_INC} + 1, \qquad \text{cols} = \frac{W\_LONG - E\_LONG}{LONG\_INC} + 1

Step-by-Step Implementation

The reader below uses only struct, numpy, and the standard library — no opaque wrapper — so that every precision decision is explicit and auditable. Build it up in four steps; each step is a runnable, type-hinted unit and each carries an inline comment naming the specification rule it satisfies.

Step 1 — Detect byte order deterministically 

import struct
import logging
from pathlib import Path
from dataclasses import dataclass

import numpy as np

logger = logging.getLogger("ntv2")

HEADER_RECORDS = 11            # NTv2 spec: 11 records per header block
RECORD_SIZE = 16              # 8-byte keyword + 8-byte value
HEADER_SIZE = HEADER_RECORDS * RECORD_SIZE   # 176-byte fixed header
NODE_SIZE = 16               # 4 x float32: lat/lon shift + lat/lon accuracy
NULL_FLAG = -999.0           # sentinel marking an unmodelled node


def detect_byte_order(raw: bytes) -> str:
    """Return '>' (big) or '<' (little) by testing NUM_OREC, which the
    spec fixes at 11. Whichever endianness decodes to 11 is the file's."""
    if raw[:8] != b"NUM_OREC":
        raise ValueError(f"Not an NTv2 file: leading keyword {raw[:8]!r}")
    if struct.unpack_from(">i", raw, 8)[0] == HEADER_RECORDS:  # big-endian probe
        return ">"
    if struct.unpack_from("<i", raw, 8)[0] == HEADER_RECORDS:  # little-endian probe
        return "<"
    raise ValueError("NUM_OREC decodes to neither 11 BE nor LE; file is corrupt")

Step 2 — Decode the master header and count sub-grids 

@dataclass(frozen=True)
class MasterHeader:
    byte_order: str
    num_overview_records: int   # NUM_OREC, always 11
    num_subgrid_records: int    # NUM_SREC, records per sub-grid header
    num_subgrids: int           # NUM_FILE, count of sub-grids to read
    gs_type: str                # GS_TYPE, e.g. 'SECONDS'


def read_master_header(raw: bytes) -> MasterHeader:
    bo = detect_byte_order(raw)
    # Numeric values sit 8 bytes after each 8-byte keyword (NTv2 record layout).
    num_orec = struct.unpack_from(bo + "i", raw, 8)[0]    # NUM_OREC
    num_srec = struct.unpack_from(bo + "i", raw, 24)[0]   # NUM_SREC
    num_file = struct.unpack_from(bo + "i", raw, 40)[0]   # NUM_FILE = sub-grid count
    gs_type = raw[48:56].decode("ascii", "replace").strip()  # GS_TYPE units flag
    if gs_type and gs_type != "SECONDS":
        # Spec permits other units, but cadastral grids are SECONDS in practice.
        logger.warning("Unexpected GS_TYPE %r; arc-second reduction assumes SECONDS", gs_type)
    return MasterHeader(bo, num_orec, num_srec, num_file, gs_type)

Step 3 — Decode each sub-grid header 

@dataclass(frozen=True)
class SubgridMetadata:
    name: str
    parent: str
    lat_min: float       # S_LAT, arc-seconds
    lat_max: float       # N_LAT, arc-seconds
    lon_min: float       # E_LONG, arc-seconds, positive west
    lon_max: float       # W_LONG, arc-seconds, positive west
    lat_spacing: float   # LAT_INC, arc-seconds
    lon_spacing: float   # LONG_INC, arc-seconds
    rows: int
    cols: int
    gs_count: int        # GS_COUNT = rows * cols


def read_subgrid_header(block: bytes, bo: str) -> SubgridMetadata:
    name = block[8:16].decode("ascii", "replace").strip()    # SUB_NAME value
    parent = block[24:32].decode("ascii", "replace").strip()  # PARENT value ('NONE' at top)
    lat_min = struct.unpack_from(bo + "d", block, 72)[0]      # S_LAT
    lat_max = struct.unpack_from(bo + "d", block, 88)[0]      # N_LAT
    lon_min = struct.unpack_from(bo + "d", block, 104)[0]     # E_LONG (+W)
    lon_max = struct.unpack_from(bo + "d", block, 120)[0]     # W_LONG (+W)
    lat_inc = struct.unpack_from(bo + "d", block, 136)[0]     # LAT_INC
    lon_inc = struct.unpack_from(bo + "d", block, 152)[0]     # LONG_INC
    gs_count = struct.unpack_from(bo + "i", block, 168)[0]    # GS_COUNT
    # Both edges are nodes, hence the +1 (NTv2 grid is inclusive of its bounds).
    rows = round((lat_max - lat_min) / lat_inc) + 1
    cols = round((lon_max - lon_min) / lon_inc) + 1
    if rows * cols != gs_count:
        raise ValueError(
            f"{name}: derived node count {rows * cols} != GS_COUNT {gs_count}"
        )
    return SubgridMetadata(
        name=name, parent=parent, lat_min=lat_min, lat_max=lat_max,
        lon_min=lon_min, lon_max=lon_max, lat_spacing=lat_inc,
        lon_spacing=lon_inc, rows=rows, cols=cols, gs_count=gs_count,
    )

Step 4 — Load the shift arrays at float64 

class NTv2Grid:
    """In-memory NTv2 grid: metadata plus float64 latitude/longitude shift
    surfaces (in arc-seconds) keyed by sub-grid name."""

    def __init__(self, gsb_path: str | Path) -> None:
        self.path = Path(gsb_path)
        self.subgrids: list[SubgridMetadata] = []
        self.lat_shift: dict[str, np.ndarray] = {}
        self.lon_shift: dict[str, np.ndarray] = {}
        self._load()

    def _load(self) -> None:
        raw = self.path.read_bytes()
        master = read_master_header(raw)
        bo = master.byte_order
        offset = HEADER_SIZE  # master header consumed; sub-grids follow
        for _ in range(master.num_subgrids):
            meta = read_subgrid_header(raw[offset:offset + HEADER_SIZE], bo)
            self.subgrids.append(meta)
            offset += HEADER_SIZE
            node_bytes = meta.gs_count * NODE_SIZE
            # 4 interleaved float32 per node: lat_shift, lon_shift, lat_acc, lon_acc.
            dtype = np.dtype(bo + "f4")
            nodes = np.frombuffer(
                raw, dtype=dtype, count=meta.gs_count * 4, offset=offset
            ).reshape(meta.gs_count, 4)
            # Promote to float64 BEFORE any interpolation so the on-disk single
            # precision never propagates into the residual budget.
            self.lat_shift[meta.name] = nodes[:, 0].astype(np.float64).reshape(meta.rows, meta.cols)
            self.lon_shift[meta.name] = nodes[:, 1].astype(np.float64).reshape(meta.rows, meta.cols)
            offset += node_bytes

    def covers(self, lat_sec: float, lon_sec_west: float) -> str | None:
        """Return the name of the smallest sub-grid covering the point, or None.
        Inputs are arc-seconds; longitude must already be positive-west."""
        match: str | None = None
        finest = float("inf")
        for m in self.subgrids:
            inside = (m.lat_min <= lat_sec <= m.lat_max
                      and m.lon_min <= lon_sec_west <= m.lon_max)
            if inside and m.lat_spacing < finest:  # prefer the densest covering child
                finest, match = m.lat_spacing, m.name
        return match

The covers method encodes the parent/child rule that gives NTv2 its accuracy: where a dense child sub-grid overlaps a coarser parent, the child’s localized correction surface must win. Selecting the parent where a child exists is a real and silent precision regression. The downstream interpolation that consumes these arrays — and the gated behaviour when covers returns None — is handled in fallback routing strategies for missing grid files.

Parameter and Return-Value Reference

Name Type Units Valid range Cadastral significance
gsb_path str | Path — existing .gsb file the authoritative grid; its version must be logged for certification
byte_order str — > or < wrong order corrupts every double; never assume
lat_min / lat_max float arc-seconds sub-grid extent defines the domain of validity; queries outside must be rejected
lon_min / lon_max float arc-seconds (+W) sub-grid extent longitude is positive-west — sign errors here invert the grid
lat_spacing / lon_spacing float arc-seconds > 0 drives interpolation kernel choice and residual magnitude
gs_count int nodes rows × cols mismatch against derived dimensions means a truncated/corrupt file
lat_shift / lon_shift np.ndarray arc-seconds, float64 finite or -999.0 the correction surface; -999.0 marks unmodelled nodes
covers() → str | None — sub-grid name None is the signal to route to a fallback, never to extrapolate

Worked Example: a Canadian NTv2_0 Grid

Consider the NAD27 → NAD83 transformation over eastern Canada, EPSG operation NAD27 to NAD83 (3) (EPSG:1313), which is realized by the NTv2_0.gsb grid. Take a point near Ottawa at latitude 45.4215° N, longitude 75.6972° W. NTv2 works in positive-west arc-seconds, so the query coordinates are lat = 45.4215 * 3600 = 163517.4″ and lon_west = 75.6972 * 3600 = 272509.92″.

grid = NTv2Grid("NTv2_0.gsb")           # the published Canadian national grid

lat_sec = 45.4215 * 3600.0
lon_west_sec = 75.6972 * 3600.0          # positive-west, per NTv2 convention

name = grid.covers(lat_sec, lon_west_sec)
print(name)                              # e.g. 'CAeast' — the covering sub-grid
print(grid.lat_shift[name].shape)        # (rows, cols) of the float64 surface
print(grid.subgrids[0].lat_spacing)      # node spacing in arc-seconds

For this region the nodal shifts are on the order of a fraction of an arc-second — roughly a few metres on the ground once reduced through the meridian radius — which is exactly the non-uniform correction a single parametric translation cannot reproduce. The numeric extraction of those node values and their per-node accuracy estimates is the dedicated routine in extracting grid metadata from .gsb files programmatically, and the full bilinear interpolation that turns the surrounding four nodes into a point shift is built in how to parse NTv2 .gsb files with Python.

Verification and Residual Analysis

A parser is only trustworthy once its output is pinned against an independent observation. After applying the interpolated shift, compute the horizontal residual between the transformed coordinate and a surveyed control monument, compare it to the agency tolerance, and emit a structured record so the result is auditable. The residual is the simple planar separation:

r=(Nout−Nctrl)2+(Eout−Ectrl)2r = \sqrt{(N_{\text{out}} - N_{\text{ctrl}})^2 + (E_{\text{out}} - E_{\text{ctrl}})^2}
import json
import math


def verify_residual(
    out_en: tuple[float, float],
    control_en: tuple[float, float],
    grid_name: str,
    tolerance_m: float = 0.020,   # 20 mm horizontal, typical cadastral gate
) -> dict[str, object]:
    """Compare a transformed point against a control monument and emit an
    audit record. Raises if the residual exceeds the survey-grade tolerance."""
    dn = out_en[0] - control_en[0]
    de = out_en[1] - control_en[1]
    residual = math.hypot(dn, de)            # planar separation in metres
    record = {
        "grid": grid_name,
        "residual_m": round(residual, 4),
        "tolerance_m": tolerance_m,
        "passed": residual <= tolerance_m,
    }
    logger.info("ntv2_verify %s", json.dumps(record))
    if not record["passed"]:
        raise ValueError(f"Residual {residual:.4f} m exceeds tolerance {tolerance_m} m")
    return record


# Expected: residual well under 20 mm against the certified monument.
rec = verify_residual((5031234.561, 445678.902), (5031234.567, 445678.910), "CAeast")
assert rec["passed"]

The structured record — grid name, residual, tolerance, and pass flag — is the minimum payload an agency submission needs. Pinning that residual against a network of monuments rather than a single point is the broader workflow in validating datum alignment with control points, and the CRS axis-order and unit setup that makes the easting/northing comparison reproducible is established in setting up high-precision coordinate reference systems.

Troubleshooting and Gotchas

Every shift value looks like astronomical garbage
The byte order is wrong. Do not hard-code >; probe NUM_OREC and accept whichever endianness decodes to 11. Canadian-origin grids are big-endian while many redistributed grids are little-endian, and a wrong order makes every float64 nonsensical without raising an exception.
Coordinates near the prime meridian or 180° drift the wrong way
NTv2 longitude is positive west and expressed in arc-seconds. A point at 75.6972° W is +272509.92″, not −272509.92″. Mixing the positive-east GIS convention with the positive-west grid convention silently mirrors the correction.
A point inside the national grid still returns no shift
It likely landed on a node carrying the -999.0 null sentinel, or inside a sub-grid hole. Treat the sentinel as "unmodelled" and reject the query; interpolating across it fabricates a shift. Confirm the point falls inside a real sub-grid extent before trusting the value.
Results are a metre off only in overlap regions
The coarse parent sub-grid was used where a dense child exists. Select the smallest-spacing covering sub-grid, not the first match — that parent/child precedence is what gives NTv2 its localized accuracy.
pyproj can't find the grid even though the file is present
When you later hand the grid to PROJ/pyproj, the PROJ_DATA environment variable (formerly PROJ_LIB) must point at the directory containing the .gsb, or PROJ will silently fall back to a lower-accuracy operation. Set it explicitly in any batch pipeline rather than relying on the install default.

Frequently Asked Questions

Why promote the shifts to float64 if the file stores float32?
The on-disk single precision is fine for storage, but bilinear or bicubic interpolation accumulates rounding. Promoting to float64 before interpolation keeps the on-disk quantisation from propagating into the residual budget, which matters at the millimetre tolerances cadastral work is graded against.
Do I need to read the accuracy columns?
Not to compute a shift, but yes for an audit-ready pipeline. Columns three and four of each node hold the latitude and longitude accuracy estimates; logging them lets you attach a defensible uncertainty to each transformed coordinate rather than asserting a blanket figure.
Can one `.gsb` contain more than one transformation?
A single file can hold many sub-grids, but they all realize one source-to-target datum pair. Different datum pairs live in different files with their own EPSG operation codes; never assume a file applies beyond the pair its master header declares.