When to Use Manual Grid Interpolation over pyproj

Manual grid interpolation earns its keep only in a handful of situations, and knowing them keeps you from reinventing an engine that already works — a judgement that belongs to Choosing a Transformation API: PROJ, pyproj, or Manual Grids within Batch Transformation & Automation for Cadastral Coordinate Pipelines. For any standard, EPSG-registered NTv2 grid, pyproj applies the correct bilinear shift at native speed and you should let it. Reading the grid and interpolating yourself is justified only when you need something the engine deliberately hides: the individual node values behind a shift, the per-node accuracy estimates, support for a grid PROJ has never heard of, or byte-for-byte deterministic control over the kernel. This page marks those boundaries and gives a self-contained interpolator that respects them.

The value is never a different coordinate. A correct manual reader applies the same bilinear kernel PROJ does, so on a registered grid it agrees to the sub-millimetre. What changes is visibility. Where pyproj returns only the shifted point, a manual reader hands you the four corner nodes, their accuracy columns, and the interpolation weights — the raw material for a defensible per-point uncertainty rather than a blanket figure. The mechanics of opening the .gsb and loading its node arrays are covered in understanding NTv2 grid shift files in Python; this page assumes that reader exists and focuses on when to reach for it.

Figure — the four corner nodes and their accuracies, which manual interpolation exposes and pyproj hides.

The Narrow Cases That Justify It #

Four situations make reading the grid yourself the right call. Per-node accuracy inspection — a QA pass that must confirm the four nodes bracketing a monument are themselves well-determined, not merely that the interpolated point looks plausible. Custom or unregistered grids — a survey authority’s provisional or internal .gsb that carries no EPSG operation code, so Transformer.from_crs cannot find it and only a file-path reader will do. Embedding uncertainty from the accuracy columns — attaching a defensible per-point standard deviation propagated from the node accuracies rather than asserting a blanket tolerance. Deterministic bilinear control — pinning the exact kernel and rounding so the result is byte-reproducible independent of the installed PROJ build.

The accuracy propagation is the case with real mathematical content. Treating the four node accuracies $\sigma_{ij}$ as independent, the bilinear weights $w_{ij}$ propagate them into the interpolated shift’s variance as

That per-point $\sigma_s$ is exactly what an engine-level API gives you no way to compute, because it never surfaces the node accuracies. If none of these four cases applies, the correct choice is pyproj, and the PROJ pipeline strings versus the pyproj Transformer API comparison decides which of its two forms to use.

Complete Runnable Implementation #

The interpolator below is self-contained given the loaded node arrays: the latitude and longitude shift surfaces plus their matching accuracy surfaces, all in arc-seconds. It rejects the -999.0 null sentinel, refuses to extrapolate past the grid edge, and returns both the shifted coordinate and the propagated uncertainty.

from __future__ import annotations
from dataclasses import dataclass

import numpy as np

NULL_FLAG = -999.0        # NTv2 sentinel: node is unmodelled, never interpolate it


@dataclass(frozen=True)
class ShiftResult:
    lon_deg: float
    lat_deg: float
    sigma_lat_sec: float      # propagated 1-sigma of the latitude shift
    sigma_lon_sec: float      # propagated 1-sigma of the longitude shift


def manual_grid_shift(
    lon_deg: float,
    lat_deg: float,
    lat_shift_sec: np.ndarray,     # latitude-shift surface, arc-seconds, float64
    lon_shift_sec: np.ndarray,     # longitude-shift surface, arc-seconds, +W
    lat_acc_sec: np.ndarray,       # latitude accuracy surface, arc-seconds
    lon_acc_sec: np.ndarray,       # longitude accuracy surface, arc-seconds
    lat_min_sec: float,
    lon_min_sec_west: float,       # NTv2 longitude is positive west
    lat_inc_sec: float,
    lon_inc_sec: float,
) -> ShiftResult:
    """Bilinearly interpolate an NTv2 shift by hand, rejecting null nodes and
    out-of-bounds queries, and propagate the per-node accuracies into a 1-sigma
    uncertainty on each shift component.

    Raises:
        ValueError: the point is outside the grid or borders a null node.
    """
    lat_sec = lat_deg * 3600.0
    lon_west_sec = -lon_deg * 3600.0            # +east degrees -> +west arc-seconds
    fj = (lat_sec - lat_min_sec) / lat_inc_sec
    fi = (lon_west_sec - lon_min_sec_west) / lon_inc_sec
    j0, i0 = int(np.floor(fj)), int(np.floor(fi))

    rows, cols = lat_shift_sec.shape
    if not (0 <= j0 < rows - 1 and 0 <= i0 < cols - 1):
        raise ValueError("query point is outside the grid; do not extrapolate")

    ty, tx = fj - j0, fi - i0                    # fractional cell position in [0, 1)
    w = (
        (1 - tx) * (1 - ty),   # w00
        tx * (1 - ty),         # w10
        (1 - tx) * ty,         # w01
        tx * ty,               # w11
    )

    def corners(surface: np.ndarray) -> tuple[float, float, float, float]:
        return (
            float(surface[j0, i0]), float(surface[j0, i0 + 1]),
            float(surface[j0 + 1, i0]), float(surface[j0 + 1, i0 + 1]),
        )

    lat_nodes = corners(lat_shift_sec)
    lon_nodes = corners(lon_shift_sec)
    if NULL_FLAG in lat_nodes or NULL_FLAG in lon_nodes:
        raise ValueError("cell borders a -999.0 null node; query is unmodelled")

    d_lat = sum(wi * ni for wi, ni in zip(w, lat_nodes))
    d_lon = sum(wi * ni for wi, ni in zip(w, lon_nodes))

    # Variance propagation: independent node accuracies, weighted bilinearly.
    lat_acc = corners(lat_acc_sec)
    lon_acc = corners(lon_acc_sec)
    var_lat = sum((wi * ai) ** 2 for wi, ai in zip(w, lat_acc))
    var_lon = sum((wi * ai) ** 2 for wi, ai in zip(w, lon_acc))

    lat_out = (lat_sec + d_lat) / 3600.0
    lon_out = -(lon_west_sec + d_lon) / 3600.0   # back to +east degrees
    return ShiftResult(
        lon_deg=lon_out,
        lat_deg=lat_out,
        sigma_lat_sec=var_lat ** 0.5,
        sigma_lon_sec=var_lon ** 0.5,
    )

Parameter Reference #

Worked Example #

import numpy as np

# A tiny 2x2 demo cell: uniform shift with modest node accuracies.
lat_shift = np.array([[1.20, 1.24], [1.22, 1.26]], dtype=np.float64)
lon_shift = np.array([[-0.80, -0.78], [-0.82, -0.79]], dtype=np.float64)
lat_acc = np.full((2, 2), 0.010, dtype=np.float64)   # 10 milli-arc-sec
lon_acc = np.full((2, 2), 0.012, dtype=np.float64)

res = manual_grid_shift(
    lon_deg=-75.6972, lat_deg=45.4215,
    lat_shift_sec=lat_shift, lon_shift_sec=lon_shift,
    lat_acc_sec=lat_acc, lon_acc_sec=lon_acc,
    lat_min_sec=45.4215 * 3600.0, lon_min_sec_west=75.6972 * 3600.0 - 30.0,
    lat_inc_sec=30.0, lon_inc_sec=30.0,
)
print(round(res.sigma_lat_sec, 4), round(res.sigma_lon_sec, 4))
# 0.0056 0.0067  -- a per-point uncertainty pyproj cannot report

Validation Check #

The propagated uncertainty must never exceed the largest contributing node accuracy, because bilinear weights sum to one and are non-negative — a cheap invariant that catches weight-sign bugs.

max_node_acc = max(lat_acc.max(), lon_acc.max())
assert res.sigma_lat_sec <= max_node_acc + 1e-12, "propagated sigma exceeds node accuracy"
assert res.sigma_lon_sec <= max_node_acc + 1e-12, "weight error in propagation"

Common Mistakes #

Related #

• Choosing a Transformation API: PROJ, pyproj, or Manual Grids — the parent decision on when each approach fits.
• PROJ pipeline strings vs the pyproj Transformer API — the sibling choice once you have decided to let PROJ do the work.
• Understanding NTv2 grid shift files in Python — the reader that loads the node arrays this page interpolates.
• Batch Transformation & Automation for Cadastral Coordinate Pipelines — the parent reference on building auditable transformation batches.