Reading GTX Geoid Model Files in Python

A .gtx file is about the simplest binary grid format in geodesy — a forty-byte header and a block of floats — which makes it a good place to get the details exactly right, because there is nothing to hide behind. This guide, part of geoid and vertical datum transformations, parses one, interpolates it, and covers the two conventions that catch every first implementation: the longitude range and the row order.

The Format #

A GTX file begins with a forty-byte, big-endian header of six values: the latitude and longitude of the south-west corner as doubles in degrees, the latitude and longitude increments as doubles in degrees, and the number of rows and columns as 32-bit integers. After it comes rows × columns 32-bit floats, in metres, ordered south to north by row and west to east within a row.

Figure — a forty-byte header and a block of floats, with two conventions to get right.

Two properties of that description are where implementations go wrong. The file is big-endian regardless of the machine reading it, so a native-order read produces values of order 10³⁸ — obviously wrong, and therefore harmless. And the longitude origin is often expressed in the 0–360 range rather than −180–180, so a file covering the conterminous United States may declare a south-west longitude of 230 rather than −130. That one is not obviously wrong: it produces a clean out-of-extent rejection for every query, which reads like a bug in the query rather than in the reader.

Complete Runnable Implementation #

from __future__ import annotations

import struct
from dataclasses import dataclass
from pathlib import Path

import numpy as np

HEADER = struct.Struct(">4d2i")          # big-endian: 4 doubles, 2 ints
NODATA = -88.8888                        # the sentinel used by several models


@dataclass(frozen=True)
class GtxGrid:
    """A parsed GTX geoid model."""

    lat0: float
    lon0: float
    dlat: float
    dlon: float
    values: np.ndarray                   # (rows, cols) float64, metres

    @classmethod
    def read(cls, path: Path) -> "GtxGrid":
        raw = path.read_bytes()
        if len(raw) < HEADER.size:
            raise ValueError(f"{path} is too short to contain a GTX header")
        lat0, lon0, dlat, dlon, rows, cols = HEADER.unpack_from(raw, 0)
        if rows <= 0 or cols <= 0:
            raise ValueError(f"{path}: implausible grid dimensions {rows}x{cols}")
        expected = HEADER.size + rows * cols * 4
        if len(raw) < expected:
            raise ValueError(
                f"{path}: header declares {rows}x{cols} nodes ({expected} bytes) "
                f"but the file is {len(raw)} bytes"
            )
        flat = np.frombuffer(raw, dtype=">f4", count=rows * cols, offset=HEADER.size)
        values = flat.reshape(rows, cols).astype(np.float64)
        values[np.isclose(values, NODATA, atol=1e-3)] = np.nan
        # Normalise a 0..360 origin to -180..180 so callers use one convention.
        if lon0 > 180.0:
            lon0 -= 360.0
        return cls(float(lat0), float(lon0), float(dlat), float(dlon), values)

    @property
    def extent(self) -> tuple[float, float, float, float]:
        rows, cols = self.values.shape
        return (self.lat0, self.lat0 + (rows - 1) * self.dlat,
                self.lon0, self.lon0 + (cols - 1) * self.dlon)

    def undulation(self, lat: float, lon: float) -> float:
        """Bilinear geoid undulation in metres. Raises outside the extent."""
        rows, cols = self.values.shape
        fi = (lat - self.lat0) / self.dlat
        fj = (lon - self.lon0) / self.dlon
        if not (0.0 <= fi <= rows - 1 and 0.0 <= fj <= cols - 1):
            s, n, w, e = self.extent
            raise ValueError(
                f"({lat:.6f}, {lon:.6f}) is outside the model extent "
                f"[{s:.3f}, {n:.3f}] x [{w:.3f}, {e:.3f}]"
            )
        i = min(int(np.floor(fi)), rows - 2)
        j = min(int(np.floor(fj)), cols - 2)
        u, v = fi - i, fj - j
        block = self.values[i:i + 2, j:j + 2]
        if not np.all(np.isfinite(block)):
            raise ValueError(
                f"({lat:.6f}, {lon:.6f}) sits on a cell with an unmodelled node"
            )
        w4 = np.array([(1 - u) * (1 - v), (1 - u) * v, u * (1 - v), u * v])
        return float(w4 @ block.reshape(4))

Parameter Reference #

Worked Example #

from pathlib import Path

grid = GtxGrid.read(Path("g2018u0.gtx"))
print(grid.extent)
# (24.0, 58.0, -130.0, -60.0)

n = grid.undulation(45.500000, -122.600000)
h_ellipsoidal = 61.230
print(f"N = {n:.3f} m, H = {h_ellipsoidal - n:.3f} m")
# N = -22.114 m, H = 83.344 m

The negative undulation is the point worth pausing on: over much of North America the geoid sits below the ellipsoid, so the orthometric height is larger than the ellipsoidal one. A sign convention taken from a European example, where the undulation is positive, produces a height wrong by twice the undulation — about 44 metres here, which is at least unmistakable.

Figure — three misreads of the format, and the height error each produces.

Validation Check #

def assert_reasonable_undulations(grid: GtxGrid) -> None:
    """Global geoid undulations lie roughly within +/-110 m of the ellipsoid."""
    finite = grid.values[np.isfinite(grid.values)]
    assert finite.size, "the model contains no valid nodes"
    assert np.abs(finite).max() < 120.0, (
        f"maximum |N| is {np.abs(finite).max():.1f} m, which is outside the "
        f"physical range — the file was probably read in the wrong byte order"
    )

A physical plausibility check is the fastest way to catch an endianness error, and it costs one pass over the array at load time.

Common Mistakes #

Reading in native byte order. GTX is big-endian by specification. On a little-endian machine a native read produces values around 10³⁸, which the plausibility assertion catches immediately — the reason to have the assertion is that without it the failure surfaces as an absurd height much later.

Assuming north-to-south row order. GTX stores the southernmost row first, the opposite of many image formats. Flipping it produces undulations mirrored about the centre latitude, which are wrong by up to tens of metres and look entirely plausible in the middle of the grid.

Forgetting the 0–360 longitude origin. A file declaring lon0 = 230.0 covers the same area as one declaring −130.0, and a reader that does not normalise rejects every query in the model’s own coverage. Normalising once at load, as above, keeps the convention question out of every call site.

Sanity-Checking a Model at Load Time #

Four checks at load time cost a single pass over the array and catch every common misread of the format.

Figure — four load-time checks, each catching a different misread.

Magnitude. Global geoid undulations lie roughly between −110 and +90 metres. Anything outside that range is an endianness or a scaling error, not geodesy.

Extent plausibility. The declared extent should cover a region that makes sense for the model’s name and file size. A model claiming a longitude span of 300 degrees is one whose origin was not normalised out of the 0–360 convention.

Node count against file size. The header’s row and column counts imply an exact byte count; a mismatch means either a truncated download or a header read at the wrong offset. Checking it turns a confusing interpolation failure into a clear file error.

Coverage. The fraction of nodes carrying the no-data sentinel says how much of the declared extent is actually modelled. A model that is largely no-data over your working area is the wrong model, and finding that out at load is much cheaper than finding it out per query.

def load_report(grid: GtxGrid) -> dict[str, float]:
    finite = np.isfinite(grid.values)
    return {
        "rows": float(grid.values.shape[0]),
        "cols": float(grid.values.shape[1]),
        "coverage": float(finite.mean()),
        "min_m": float(np.nanmin(grid.values)),
        "max_m": float(np.nanmax(grid.values)),
    }

Logging that dictionary once per model load gives every later question about a height an answer that is already on the record.

Frequently Asked Questions #

Is GTX still current? #

It is widely deployed and increasingly superseded by GeoTIFF-based grids, which carry their own CRS and metadata rather than relying on a forty-byte header. A cadastral pipeline reads legacy formats for decades, so supporting GTX remains worthwhile; new work should prefer the self-describing format where the model is published in both.

Should the undulation be interpolated bilinearly or bicubically? #

Follow the model’s own guidance where it gives any: geoid surfaces are smooth, and several national models specify a cubic kernel, where the accuracy gain is real. Where nothing is specified, bilinear is the safe default and its error at typical one-arc-minute spacing is a few millimetres — the comparison is in bilinear vs bicubic interpolation for grid shifts.

What uncertainty should I attach to the interpolated value? #

The model’s published accuracy, which for a recent national geoid is typically one to three centimetres, plus the interpolation term. That total is usually the largest single contribution to a GNSS-derived orthometric height — larger than the GNSS observation — and omitting it makes the height look better than it is.

Can I use a GTX file outside its stated extent? #

No. Extrapolating a geoid model past its published extent produces a number with no defensible uncertainty, and the reader above raises rather than clamping for exactly that reason. Route the query to a global model with its own, larger uncertainty and record that the fallback was used.

Should the parsed grid be cached between runs? #

Yes, and cache the parsed array rather than re-reading the file: the parse is cheap but the byte-order conversion over a large national model is not free, and a pipeline that reads it per point will spend most of its time there. Key the cache on the file checksum rather than the path, so a re-issued model invalidates it automatically.

Related #

• Geoid and vertical datum transformations — the parent topic
• Applying geoid undulation for orthometric heights
• Detecting null shift sentinels in grid files
• How to parse NTv2 GSB files with Python
• Interpolation methods for grid shift surfaces