CoordTransform Cadastral & High-Precision Python

How to Parse NTv2 .gsb Files with Python

Parsing an NTv2 .gsb file means reading its master header, walking every subgrid, and loading the per-node latitude and longitude shift surfaces into memory as deterministic float64 arrays — exactly, with no byte drift — so a cadastral datum transformation reproduces monument coordinates to survey-grade tolerance. This is the full read operation that extends understanding NTv2 grid shift files in Python, the parent guide on the .gsb binary structure, and it underpins every grid-based shift defined under Core Transformation Fundamentals & Standards. Where extracting grid metadata from .gsb files programmatically stops at the headers, this page reads the shift block itself, which is where a single mis-aligned seek silently corrupts a deliverable.

Deterministic .gsb parsing flow with byte-order detection and per-subgrid validation A top-to-bottom flow. The parser opens the .gsb binary, then a decision checks whether there are enough bytes for the master header; if not, control routes to a fallback that returns an empty subgrid list. If yes, the parser detects byte order by the invariant NUM_OREC equals 11, then loops over each subgrid reading its header and shift arrays. A second decision asks whether rows times cols equals GS_COUNT and the shifts are plausible: if yes the subgrid is validated, if no it is flagged and the parse fails. Validated subgrids are returned. Open .gsb (binary) enough bytes for master header? no Route to fallback (empty subgrid list) yes Detect byte order NUM_OREC == 11 For each sub-grid: read header + shift arrays rows·cols == GS_COUNT? shifts plausible? yes no VALIDATED sub-grid FLAGGED / fail Return sub-grids

Figure — deterministic .gsb parsing with byte-order detection, geometry validation, and fallback routing.

What the Shift Block Actually Stores

After each 176-byte subgrid header (11 records of 16 bytes: an 8-byte ASCII label plus an 8-byte value) comes the shift block: GS_COUNT nodes, each a 16-byte record of four IEEE-754 single-precision floats — latitude shift, longitude shift, latitude accuracy, longitude accuracy. The bounds (S_LAT, N_LAT, E_LONG, W_LONG) and the shifts are in arc-seconds when GS_TYPE is SECONDS, and NTv2 follows the geodetic convention that longitude increases westward, so E_LONG < W_LONG numerically. Nodes are stored row by row from the lower-left corner: longitude varies fastest from the east bound toward the west bound, then latitude steps north.

Byte layout of a .gsb file and the row-major node order The top bar splits the file into three regions: a master header of 176 bytes (11 records), a subgrid header of 176 bytes (11 records), and a shift block of GS_COUNT nodes of 16 bytes each. A header record expands into an 8-byte ASCII label followed by an 8-byte value decoded as a float64 or int. A shift node expands into four 4-byte single-precision floats: latitude shift, longitude shift, latitude accuracy, and longitude accuracy. Decoding the node block as 8-byte doubles instead of 4-byte floats misaligns every value. The lower panel shows a 3 by 3 grid whose nodes are numbered 0 to 8 in storage order, running east to west across each row from the south-west corner, then stepping north to the next row. Master header 176 B · 11 records Subgrid header 176 B · 11 records Shift block GS_COUNT × 16 B nodes 8-byte ASCII label 8-byte value float64 / int one 16-byte header record lat shift f32 · 4 B lon shift f32 · 4 B lat acc f32 · 4 B lon acc f32 · 4 B one 16-byte shift node = 4 × float32 Node storage order — row-major from the SW corner longitude: east → west (varies fastest) · latitude: south → north (row steps) dashed = wrap to next row N 0 1 2 3 4 5 6 7 8 SW NE

Figure — the .gsb byte map: 176-byte headers of 8+8 records, then a shift block of 16-byte nodes (four float32), stored row-major from the south-west corner.

The grid geometry is fully determined by the header, and the derived node count must equal the stored GS_COUNT or the file is corrupt or the byte order was misread:

nϕ=ϕN−ϕSΔϕ+1,nλ=λW−λEΔλ+1,N=nϕ nλn_{\phi} = \frac{\phi_{N}-\phi_{S}}{\Delta\phi}+1, \qquad n_{\lambda} = \frac{\lambda_{W}-\lambda_{E}}{\Delta\lambda}+1, \qquad N = n_{\phi}\,n_{\lambda}

A shift is an angular correction, so converting it to ground distance needs the local ellipsoidal radii — the meridian radius M(ϕ)M(\phi) for the latitude component and the prime-vertical radius N(ϕ)N(\phi) for the longitude component, as derived in projection math fundamentals for cadastral surveys:

Δn=δϕ3600⋅π180 M(ϕ),Δe=δλ3600⋅π180 N(ϕ)cos⁡ϕ\Delta n = \frac{\delta\phi}{3600}\cdot\frac{\pi}{180}\,M(\phi), \qquad \Delta e = \frac{\delta\lambda}{3600}\cdot\frac{\pi}{180}\,N(\phi)\cos\phi

One more rule dominates a correct reader: -999.0 is the null-shift sentinel. A node carrying it is unmodelled, and interpolating across it instead of rejecting the query is a classic out-of-extent defect — so the parser preserves it as NaN rather than averaging it away.

Complete Runnable Implementation

The parser below reads a .gsb buffer, auto-detects byte order from the NUM_OREC == 11 invariant (NTv2 ships in both big- and little-endian), walks every subgrid, validates rows · cols == GS_COUNT, and returns frozen, strictly typed subgrid objects with the shift surfaces as float64 NumPy arrays. The pure parse_ntv2_shifts works on an in-memory buffer so it is independently testable; read_ntv2_gsb is the thin file-reading wrapper.

from __future__ import annotations

import struct
from dataclasses import dataclass

import numpy as np

# NTv2 master and subgrid headers are each 11 records of 16 bytes
# (8-byte ASCII label + 8-byte value). After a subgrid header come GS_COUNT
# nodes; each node is four float32: lat shift, lon shift, lat acc, lon acc.
# GS_TYPE "SECONDS" => bounds/shifts in arc-seconds, longitude positive WEST.
RECORD_SIZE = 16
HEADER_RECORDS = 11
HEADER_SIZE = RECORD_SIZE * HEADER_RECORDS    # 176 bytes
NODE_SIZE = 16                                # 4 x float32
NULL_SHIFT = -999.0                           # unmodelled-node sentinel


@dataclass(frozen=True)
class NTv2Subgrid:
    name: str
    parent: str
    s_lat: float             # southern bound (arc-seconds, +ve north)
    n_lat: float             # northern bound (arc-seconds)
    e_long: float            # eastern bound  (arc-seconds, +ve WEST)
    w_long: float            # western bound  (arc-seconds, +ve WEST)
    lat_inc: float           # latitude node spacing  (arc-seconds)
    long_inc: float          # longitude node spacing (arc-seconds)
    rows: int                # n_phi
    cols: int                # n_lambda
    lat_shift: np.ndarray    # float64 arc-seconds, shape (rows, cols)
    lon_shift: np.ndarray    # float64 arc-seconds, shape (rows, cols)


def _value_offset(record_index: int, base: int = 0) -> int:
    """Byte offset of a record's value field (8-byte label precedes it)."""
    return base + record_index * RECORD_SIZE + 8


def _ascii(buf: bytes, offset: int) -> str:
    """Decode an 8-byte ASCII value field, trimming spaces and NULs."""
    return buf[offset:offset + 8].decode("ascii", "replace").rstrip(" \x00")


def _detect_byte_order(buf: bytes) -> str:
    """NUM_OREC is always 11; whichever order yields 11 is the file's order."""
    for order in (">", "<"):
        if struct.unpack_from(order + "i", buf, _value_offset(0))[0] == HEADER_RECORDS:
            return order
    raise ValueError("NUM_OREC != 11 under either byte order; not valid NTv2.")


def parse_ntv2_shifts(buf: bytes) -> list[NTv2Subgrid]:
    """Parse every subgrid's shift surfaces from an in-memory .gsb buffer."""
    if len(buf) < 2 * HEADER_SIZE:
        raise ValueError("Buffer too short for a master + subgrid header.")
    order = _detect_byte_order(buf)
    num_srec = struct.unpack_from(order + "i", buf, _value_offset(1))[0]  # NUM_SREC

    subgrids: list[NTv2Subgrid] = []
    offset = HEADER_SIZE  # skip the 176-byte master header
    for _ in range(num_srec):
        if offset + HEADER_SIZE > len(buf):
            raise ValueError("Truncated before a declared subgrid header.")
        b = offset
        name = _ascii(buf, _value_offset(0, b))
        parent = _ascii(buf, _value_offset(1, b))
        s_lat = struct.unpack_from(order + "d", buf, _value_offset(4, b))[0]
        n_lat = struct.unpack_from(order + "d", buf, _value_offset(5, b))[0]
        e_long = struct.unpack_from(order + "d", buf, _value_offset(6, b))[0]
        w_long = struct.unpack_from(order + "d", buf, _value_offset(7, b))[0]
        lat_inc = struct.unpack_from(order + "d", buf, _value_offset(8, b))[0]
        long_inc = struct.unpack_from(order + "d", buf, _value_offset(9, b))[0]
        gs_count = struct.unpack_from(order + "i", buf, _value_offset(10, b))[0]

        # Geometry is fixed by the header; it must agree with GS_COUNT.
        rows = round((n_lat - s_lat) / lat_inc) + 1
        cols = round((w_long - e_long) / long_inc) + 1
        if rows * cols != gs_count:
            raise ValueError(
                f"{name}: rows*cols {rows * cols} != GS_COUNT {gs_count} "
                "(corrupt header or misread byte order)."
            )

        node_start = offset + HEADER_SIZE
        node_bytes = gs_count * NODE_SIZE
        if node_start + node_bytes > len(buf):
            raise ValueError(f"{name}: shift block truncated.")

        # Read all four float32 columns at the file's byte order, then promote
        # to float64 so downstream interpolation does not leak single-precision.
        nodes = np.frombuffer(
            buf[node_start:node_start + node_bytes], dtype=order + "f4"
        ).reshape(gs_count, 4).astype(np.float64)

        # Nodes run south->north (rows), east->west (cols, +ve west).
        lat_shift = nodes[:, 0].reshape(rows, cols)
        lon_shift = nodes[:, 1].reshape(rows, cols)

        # Preserve -999.0 as NaN so it cannot be silently interpolated across.
        lat_shift = np.where(lat_shift == NULL_SHIFT, np.nan, lat_shift)
        lon_shift = np.where(lon_shift == NULL_SHIFT, np.nan, lon_shift)

        subgrids.append(NTv2Subgrid(
            name=name, parent=parent, s_lat=s_lat, n_lat=n_lat,
            e_long=e_long, w_long=w_long, lat_inc=lat_inc, long_inc=long_inc,
            rows=rows, cols=cols, lat_shift=lat_shift, lon_shift=lon_shift,
        ))
        offset = node_start + node_bytes
    return subgrids


def read_ntv2_gsb(path: str) -> list[NTv2Subgrid]:
    """Read a .gsb file from disk and return its parsed subgrids."""
    with open(path, "rb") as fh:
        return parse_ntv2_shifts(fh.read())

Inline Parameter Reference

Every field a parsed subgrid exposes, with its unit and the range a valid continental grid stays within. Bounds and spacing come from the subgrid header; the shift surfaces come from the node block that follows it.

Field Type Units Valid range Significance
s_lat / n_lat float arc-seconds −324000-324000–324000324000 Southern/northern bound (+ve north)
e_long / w_long float arc-seconds 00–12960001296000 Eastern/western bound (+ve west)
lat_inc / long_inc float arc-seconds >0> 0 Node spacing; must divide the extent
rows / cols int nodes ≥2\ge 2 nϕn_\phi / nλn_\lambda; product equals GS_COUNT
lat_shift ndarray arc-seconds ≈±60\approx \pm 60 Latitude correction per node; NaN = unmodelled
lon_shift ndarray arc-seconds ≈±60\approx \pm 60 Longitude correction per node; NaN = unmodelled

Minimal Worked Example

Rather than ship a binary fixture, the example packs a valid NAD27→NAD83 master header, one subgrid header, and a tiny 3×3 shift block in memory, then parses it straight back — so the block runs end to end with no external .gsb file.

import struct


def build_demo_gsb(order: str = "<") -> bytes:
    rec = lambda label, val: label.ljust(8).encode("ascii")[:8] + val
    pad_int = lambda n: struct.pack(order + "i", n) + b"\x00\x00\x00\x00"
    dbl = lambda x: struct.pack(order + "d", x)
    txt = lambda t: t.ljust(8).encode("ascii")[:8]

    s_lat, lat_inc = 90000.0, 900.0      # 25.00 N, 0.25 deg spacing
    e_long, long_inc = 241200.0, 900.0   # 67.00 W, +ve west
    rows = cols = 3
    n_lat = s_lat + (rows - 1) * lat_inc
    w_long = e_long + (cols - 1) * long_inc

    master = b"".join([
        rec("NUM_OREC", pad_int(11)), rec("NUM_SREC", pad_int(1)),
        rec("NUM_FILE", pad_int(1)), rec("GS_TYPE ", txt("SECONDS")),
        rec("VERSION ", txt("NTv2.0")), rec("SYSTEM_F", txt("NAD27")),
        rec("SYSTEM_T", txt("NAD83")), rec("MAJOR_F ", dbl(6378206.4)),
        rec("MINOR_F ", dbl(6356583.8)), rec("MAJOR_T ", dbl(6378137.0)),
        rec("MINOR_T ", dbl(6356752.314140)),
    ])
    subgrid = b"".join([
        rec("SUB_NAME", txt("NA0")), rec("PARENT  ", txt("NONE")),
        rec("CREATED ", txt("20240101")), rec("UPDATED ", txt("20240101")),
        rec("S_LAT   ", dbl(s_lat)), rec("N_LAT   ", dbl(n_lat)),
        rec("E_LONG  ", dbl(e_long)), rec("W_LONG  ", dbl(w_long)),
        rec("LAT_INC ", dbl(lat_inc)), rec("LONG_INC", dbl(long_inc)),
        rec("GS_COUNT", pad_int(rows * cols)),
    ])
    nodes = bytearray()
    for i in range(rows * cols):            # lat_sh, lon_sh, lat_acc, lon_acc
        nodes += struct.pack(order + "4f", 1.50 + 0.01 * i, -4.20 - 0.01 * i, 0.05, 0.05)
    return master + subgrid + bytes(nodes)


grids = parse_ntv2_shifts(build_demo_gsb("<"))
g = grids[0]
print(f"{g.name}: {g.rows}x{g.cols} nodes, parent={g.parent}")
print("SW lat shift (arc-sec):", round(float(g.lat_shift[0, 0]), 3))
print("NE lon shift (arc-sec):", round(float(g.lon_shift[-1, -1]), 3))
# NA0: 3x3 nodes, parent=NONE
# SW lat shift (arc-sec): 1.5
# NE lon shift (arc-sec): -4.28

Validation Check

Reading the block is only trustworthy once the values it reports are physically plausible: continental NTv2 corrections stay well under one arc-minute, and no NaN sentinel should slip into a region you intend to interpolate. The assertion below fails loudly on either condition — the same gate a survey-grade pipeline runs before handing the surfaces to an interpolator.

import numpy as np

MAX_PLAUSIBLE_SHIFT_ARCSEC = 60.0   # national grids stay under ~1 arc-minute


def assert_shifts_plausible(sg: NTv2Subgrid) -> None:
    """Survey-grade gate: finite shifts must stay within a plausible bound."""
    finite = np.concatenate([
        sg.lat_shift[np.isfinite(sg.lat_shift)].ravel(),
        sg.lon_shift[np.isfinite(sg.lon_shift)].ravel(),
    ])
    peak = float(np.abs(finite).max())
    assert peak <= MAX_PLAUSIBLE_SHIFT_ARCSEC, \
        f'{sg.name}: peak shift {peak:.3f}" exceeds plausible bound'


assert_shifts_plausible(g)   # passes for a well-formed grid

Common Mistakes

Reading the shift nodes as float64 instead of float32
The bounds and spacing in the header are IEEE-754 doubles, but every shift node is four single-precision floats. Decoding the node block with <d (8 bytes) instead of <f4 (4 bytes) halves the node count, misaligns every subsequent value, and produces denormal garbage. Read the shift block as f4 at the file's byte order, then promote to float64 for downstream arithmetic.
Hard-coding big-endian byte order
NTv2 files ship in both orders — Canadian originals differ from many derived national grids. A fixed > flag will decode NUM_OREC as a garbage integer and the bounds as denormals. Detect the order from the NUM_OREC == 11 invariant before reading anything, and carry the detected order through every struct.unpack_from and the numpy dtype.
Interpolating across the -999.0 null sentinel
A node value of -999.0 marks an unmodelled cell, not a 999-arc-second correction. Treating it as a real shift drags a bilinear interpolation toward a nonsense value near the grid edge. Convert the sentinel to NaN on read so any query that touches an unmodelled node fails the extent check instead of returning a plausible-looking but invalid result.

Frequently Asked Questions

How do I locate a single node's shift from a latitude and longitude?
Convert the query to arc-seconds (longitude positive west), then compute the row as (phi - s_lat) / lat_inc and the column as (w_long_query - e_long) / long_inc. The integer parts index the lower-left node of the enclosing cell; the fractional parts feed the bilinear weights. Reject the query if either index falls outside 0 .. rows-1 / 0 .. cols-1.
Should I load the whole shift block or memory-map it?
A continental grid is a few megabytes, so reading it whole is fine and keeps the parser pure and testable. For very large multi-subgrid national files, numpy.memmap over the node region avoids copying, but you still validate rows*cols == GS_COUNT against the header before trusting any slice.
What should the parser do when a subgrid is corrupt or truncated?
Raise immediately rather than returning a partial surface — a truncated block would index past the buffer and silently mislabel later subgrids. The caller then routes to a known-good grid or a parametric method, as described in the fallback strategy.