CoordTransform Cadastral & High-Precision Python

Extracting Grid Metadata from .gsb Files Programmatically

Extracting grid metadata from a .gsb file means reading the NTv2 master header and the first subgrid header to recover the datum pair, ellipsoid axes, spatial extent, and node spacing — bit-for-bit — before a single coordinate is shifted. This is the validation entry point for understanding NTv2 grid shift files in Python and, more broadly, for any audit-ready pipeline built on Core Transformation Fundamentals & Standards. The operation must satisfy ISO 19111:2019 §C.5, which requires every grid-based coordinate operation to expose its source and target reference frames and its domain of validity. A metadata extractor that drifts by even one byte offset silently mislabels the datum or clips the subgrid, so the tolerance here is exact equality, not an approximation.

Extracting `.gsb` grid metadata, step by step The extractor reads the first 352 bytes (the 176-byte master header plus the 176-byte first subgrid header), detects byte order by checking that NUM_OREC decodes to 11 — raising if neither endianness works — decodes the master header records at offset record times 16 plus 8, jumps to byte 176 to read the subgrid extent and spacing, then gates on the derived node count matching GS_COUNT before emitting a frozen, typed GsbMetadata object; a mismatch or unreadable header raises instead. Read first 352 bytes master header 176 B + first subgrid 176 B NUM_OREC == 11? try big-endian then little yes neither raise ValueError not valid NTv2 Decode master header (records 0–10) value field at record × 16 + 8 counts · GS_TYPE · SYSTEM_F/T · ellipsoid axes Jump to byte 176 — first subgrid header S_LAT N_LAT E_LONG W_LONG LAT_INC LONG_INC GS_COUNT n_φ · n_λ == GS_COUNT ? no raise — corrupt or byte order misread yes Emit GsbMetadata frozen · strictly typed

How the .gsb Header Is Laid Out

The NTv2 format does not store metadata as a flat array of numbers. Both the master header and every subgrid header are 11 records of 16 bytes, and each record is split into an 8-byte ASCII label followed by an 8-byte value. The master header’s records carry NUM_OREC, NUM_SREC, NUM_FILE, GS_TYPE, VERSION, SYSTEM_F, SYSTEM_T, and the four ellipsoid axes MAJOR_F/MINOR_F/MAJOR_T/MINOR_T. The spatial extent and node spacing do not live in the master header — they live in the subgrid header that begins at byte 176, in the records S_LAT, N_LAT, E_LONG, W_LONG, LAT_INC, LONG_INC, and GS_COUNT. Count fields are 4-byte integers padded to the 8-byte value slot; axes and bounds are IEEE-754 doubles. Two conventions trip up almost every first parser: NTv2 expresses latitude and longitude bounds in arc-seconds (when GS_TYPE is SECONDS), and longitude is positive west, so E_LONG is the numerically smaller western-positive bound and W_LONG the larger.

The two 176-byte NTv2 headers, record by record The master header occupies bytes 0 to 175 as eleven 16-byte records carrying NUM_OREC, NUM_SREC, NUM_FILE, GS_TYPE, VERSION, SYSTEM_F, SYSTEM_T and the four ellipsoid axes MAJOR_F, MINOR_F, MAJOR_T and MINOR_T. The first subgrid header begins at byte 176 with SUB_NAME, PARENT, CREATED, UPDATED and then the geometry records S_LAT, N_LAT, E_LONG, W_LONG, LAT_INC, LONG_INC and GS_COUNT. Each record is an 8-byte ASCII label followed by an 8-byte value, so 11 records make up a 176-byte header. Master header — bytes 0–175 NUM_OREC NUM_SREC NUM_FILE GS_TYPE VERSION SYSTEM_F SYSTEM_T MAJOR_F MINOR_F MAJOR_T MINOR_T First subgrid header — bytes 176–351 SUB_NAME PARENT CREATED UPDATED S_LAT N_LAT E_LONG W_LONG LAT_INC LONG_INC GS_COUNT 8-byte ASCII label e.g. "S_LAT " 8-byte value IEEE-754 double / 4-byte int Each record = 16 bytes; 11 records = one 176-byte header.

Once the bounds are decoded, the grid geometry is fully determined. The node counts along each axis follow directly from the extent and spacing, and their product must equal the stored GS_COUNT:

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 non-integer nϕn_\phi or nλn_\lambda, or a product that disagrees with GS_COUNT, means the header is corrupt or the byte order was misread — both must fail extraction rather than propagate into the interpolation stage.

Complete Runnable Implementation

The parser below reads the master header and the first subgrid header, auto-detects byte order from the NUM_OREC == 11 invariant (NTv2 files ship in both big- and little-endian), and returns a frozen, strictly typed metadata object. The pure parse_gsb_metadata function works on an in-memory buffer so it is independently testable; extract_gsb_metadata is the thin file-reading wrapper.

from __future__ import annotations

import struct
from dataclasses import dataclass

# NTv2 master and subgrid headers are each 11 records of 16 bytes
# (8-byte ASCII label + 8-byte value). GS_TYPE "SECONDS" => bounds in
# arc-seconds, longitude positive WEST. ISO 19111:2019 §C.5 requires a
# grid-based operation to expose source/target frames and its extent.
RECORD_SIZE = 16
HEADER_RECORDS = 11
MASTER_HEADER_SIZE = RECORD_SIZE * HEADER_RECORDS    # 176 bytes
SUBGRID_HEADER_SIZE = RECORD_SIZE * HEADER_RECORDS   # 176 bytes
SECONDS_PER_DEGREE = 3600.0


@dataclass(frozen=True)
class GsbMetadata:
    byte_order: str          # ">" big-endian or "<" little-endian
    overview_records: int    # NUM_OREC (always 11)
    subgrid_records: int     # NUM_SREC
    file_count: int          # NUM_FILE (number of subgrids)
    shift_unit: str          # GS_TYPE, e.g. "SECONDS"
    version: str             # VERSION
    source_datum: str        # SYSTEM_F
    target_datum: str        # SYSTEM_T
    major_f: float           # source ellipsoid semi-major axis (m)
    minor_f: float           # source ellipsoid semi-minor axis (m)
    major_t: float           # target ellipsoid semi-major axis (m)
    minor_t: float           # target ellipsoid semi-minor axis (m)
    sub_name: str            # first subgrid SUB_NAME
    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)
    gs_count: int            # node count of the first subgrid


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 _double(buf: bytes, order: str, offset: int) -> float:
    """Unpack an 8-byte IEEE-754 double from a record value field."""
    return struct.unpack_from(order + "d", buf, offset)[0]


def _int(buf: bytes, order: str, offset: int) -> int:
    """NTv2 stores record counts as a 4-byte int padded to 8 bytes."""
    return struct.unpack_from(order + "i", buf, offset)[0]


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 _int(buf, order, _value_offset(0)) == HEADER_RECORDS:
            return order
    raise ValueError("NUM_OREC != 11 under either byte order; not valid NTv2.")


def parse_gsb_metadata(buf: bytes) -> GsbMetadata:
    """Parse the master header and first subgrid header from a .gsb buffer."""
    if len(buf) < MASTER_HEADER_SIZE + SUBGRID_HEADER_SIZE:
        raise ValueError("Buffer too short for a master + subgrid header.")
    order = _detect_byte_order(buf)

    # --- master header: records 0..10 ---
    num_orec = _int(buf, order, _value_offset(0))
    num_srec = _int(buf, order, _value_offset(1))
    num_file = _int(buf, order, _value_offset(2))
    gs_type = _ascii(buf, _value_offset(3))
    version = _ascii(buf, _value_offset(4))
    system_f = _ascii(buf, _value_offset(5))
    system_t = _ascii(buf, _value_offset(6))
    major_f = _double(buf, order, _value_offset(7))
    minor_f = _double(buf, order, _value_offset(8))
    major_t = _double(buf, order, _value_offset(9))
    minor_t = _double(buf, order, _value_offset(10))

    # --- first subgrid header: records 0..10, starting at byte 176 ---
    b = MASTER_HEADER_SIZE
    sub_name = _ascii(buf, _value_offset(0, b))
    s_lat = _double(buf, order, _value_offset(4, b))
    n_lat = _double(buf, order, _value_offset(5, b))
    e_long = _double(buf, order, _value_offset(6, b))
    w_long = _double(buf, order, _value_offset(7, b))
    lat_inc = _double(buf, order, _value_offset(8, b))
    long_inc = _double(buf, order, _value_offset(9, b))
    gs_count = _int(buf, order, _value_offset(10, b))

    return GsbMetadata(
        byte_order=order, overview_records=num_orec, subgrid_records=num_srec,
        file_count=num_file, shift_unit=gs_type, version=version,
        source_datum=system_f, target_datum=system_t,
        major_f=major_f, minor_f=minor_f, major_t=major_t, minor_t=minor_t,
        sub_name=sub_name, s_lat=s_lat, n_lat=n_lat, e_long=e_long,
        w_long=w_long, lat_inc=lat_inc, long_inc=long_inc, gs_count=gs_count,
    )


def extract_gsb_metadata(path: str) -> GsbMetadata:
    """Read the headers of a .gsb file and return its grid metadata."""
    with open(path, "rb") as fh:
        buf = fh.read(MASTER_HEADER_SIZE + SUBGRID_HEADER_SIZE)
    return parse_gsb_metadata(buf)

Inline Parameter Reference

Every field the extractor returns, with its unit and the range a valid continental grid stays within. Datum and ellipsoid fields come from the master header; extent and spacing from the first subgrid header.

Field Type Units Valid range Significance
byte_order str — ">" or "<" Endianness recovered from NUM_OREC == 11
shift_unit str — SECONDS / MINUTES / DEGREES GS_TYPE; fixes how bounds are scaled
source_datum / target_datum str — 8-char label SYSTEM_F / SYSTEM_T, resolve to EPSG codes
major_t / minor_t float metres 6.35e6–6.40e6 Target ellipsoid axes; checked against EPSG
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
gs_count int nodes nϕ⋅nλn_\phi \cdot n_\lambda Length of the subgrid shift block

Minimal Worked Example

Rather than ship a binary fixture, the example packs a valid NAD27→NAD83 header in memory (Clarke 1866 source, GRS 1980 target, a quarter-degree continental grid) and parses it straight back — so the block runs end to end with no external .gsb file.

import struct


def _record(label: str, value: bytes) -> bytes:
    return label.ljust(8).encode("ascii")[:8] + value


def build_demo_header(order: str = "<") -> bytes:
    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]

    master = b"".join([
        _record("NUM_OREC", pad_int(11)),
        _record("NUM_SREC", pad_int(1)),
        _record("NUM_FILE", pad_int(1)),
        _record("GS_TYPE ", txt("SECONDS")),
        _record("VERSION ", txt("NTv2.0")),
        _record("SYSTEM_F", txt("NAD27")),
        _record("SYSTEM_T", txt("NAD83")),
        _record("MAJOR_F ", dbl(6378206.4)),         # Clarke 1866
        _record("MINOR_F ", dbl(6356583.8)),
        _record("MAJOR_T ", dbl(6378137.0)),         # GRS 1980
        _record("MINOR_T ", dbl(6356752.314140)),
    ])

    s_lat, n_lat = 90000.0, 180000.0       # 25.0 N .. 50.0 N  (arc-seconds)
    e_long, w_long = 241200.0, 450000.0    # 67.0 W .. 125.0 W (+ve west)
    lat_inc = long_inc = 900.0             # 0.25 deg node spacing
    rows = int((n_lat - s_lat) / lat_inc) + 1
    cols = int((w_long - e_long) / long_inc) + 1

    subgrid = b"".join([
        _record("SUB_NAME", txt("NA0")),
        _record("PARENT  ", txt("NONE")),
        _record("CREATED ", txt("20240101")),
        _record("UPDATED ", txt("20240101")),
        _record("S_LAT   ", dbl(s_lat)),
        _record("N_LAT   ", dbl(n_lat)),
        _record("E_LONG  ", dbl(e_long)),
        _record("W_LONG  ", dbl(w_long)),
        _record("LAT_INC ", dbl(lat_inc)),
        _record("LONG_INC", dbl(long_inc)),
        _record("GS_COUNT", pad_int(rows * cols)),
    ])
    return master + subgrid


meta = parse_gsb_metadata(build_demo_header(order="<"))
print(meta.source_datum, "->", meta.target_datum)
print("unit:", meta.shift_unit, "| byte order:", meta.byte_order)
print(f"lat {meta.s_lat / 3600:.3f}..{meta.n_lat / 3600:.3f} deg, "
      f"W-long {meta.e_long / 3600:.3f}..{meta.w_long / 3600:.3f} deg")
print("nodes:", meta.gs_count)
# NAD27 -> NAD83
# unit: SECONDS | byte order: <
# lat 25.000..50.000 deg, W-long 67.000..125.000 deg
# nodes: 23533

Validation Check

Extraction is only trustworthy once the geometry it reports is internally consistent: the target ellipsoid must match its EPSG-registered axis, the spacing must divide the extent into whole cells, and the derived node count must equal the stored GS_COUNT. The assertion below fails loudly on any of these, the same gate a survey-grade pipeline runs before ingesting the shift block.

import math

GRS80_SEMI_MAJOR_M = 6378137.0
EPSG_AXIS_TOL_M = 1e-4    # 0.1 mm on the semi-major axis


def assert_metadata_consistent(meta: GsbMetadata) -> None:
    """Survey-grade gate: ellipsoid, spacing, and node count must agree."""
    assert math.isclose(meta.major_t, GRS80_SEMI_MAJOR_M, abs_tol=EPSG_AXIS_TOL_M), \
        f"target semi-major {meta.major_t} != GRS80 {GRS80_SEMI_MAJOR_M}"
    lat_cells = (meta.n_lat - meta.s_lat) / meta.lat_inc
    lon_cells = (meta.w_long - meta.e_long) / meta.long_inc
    assert math.isclose(lat_cells, round(lat_cells), abs_tol=1e-6), "lat spacing"
    assert math.isclose(lon_cells, round(lon_cells), abs_tol=1e-6), "lon spacing"
    rows, cols = round(lat_cells) + 1, round(lon_cells) + 1
    assert meta.gs_count == rows * cols, \
        f"GS_COUNT {meta.gs_count} != rows*cols {rows * cols}"


assert_metadata_consistent(meta)   # passes for a well-formed grid

Common Mistakes

Packing every field as a bare double and ignoring the 8-byte labels
Each NTv2 record is an 8-byte ASCII label followed by an 8-byte value, and the spatial bounds live in the subgrid header at byte 176 — not in the master header. Reading the header as a flat run of doubles from offset 0 mislabels the datum, reads ellipsoid axes where a label string sits, and invents a non-existent extent. Always index the value field as record * 16 + 8 and read S_LAT/N_LAT/E_LONG/W_LONG from the subgrid block.
Hard-coding big-endian byte order
NTv2 files are distributed in both byte orders — Canadian originals differ from many derived national grids. A fixed >d flag will decode NUM_OREC as a garbage integer and the axes as denormals. Detect the order from the NUM_OREC == 11 invariant before reading anything else, and carry the detected order through every struct.unpack_from call.
Forgetting that longitude is positive west and bounds are in arc-seconds
When GS_TYPE is SECONDS, every bound is in arc-seconds and longitude increases westward, so E_LONG < W_LONG numerically. Treating the values as signed east-positive degrees flips the sign and rescales by 3600, displacing the grid by tens of degrees. Divide by 3600 for degrees and negate the longitude only when handing positions to an east-positive consumer.

Frequently Asked Questions

Do I need to read every subgrid header to validate a file?
For a metadata sanity check the master header plus the first subgrid header is enough to confirm the datum pair, byte order, and extent. For a full integrity pass, walk all NUM_FILE subgrids — each has its own 176-byte header and GS_COUNT shift block — and confirm the PARENT labels form a consistent nesting tree.
How do I map SYSTEM_F and SYSTEM_T to EPSG codes?
The 8-character labels (for example NAD27, NAD83) are free text, not codes. Maintain an explicit lookup from label to EPSG datum and cross-check the reported ellipsoid axes against the registered values, so a mislabelled file cannot masquerade as the wrong datum.
What should extraction do when GS_COUNT disagrees with the extent?
Treat it as a hard failure. A mismatch means either the byte order was misread or the file is truncated, and continuing would index past the shift block. Raise before any interpolation, then route to a known-good grid as described in the fallback strategy.