Combining Grid and Parameter Uncertainty in Python

A grid-based shift and a parametric one fail in different ways, and a pipeline that can use either has to combine their uncertainties without pretending they are the same kind of number. This guide, part of uncertainty propagation through transformation chains, assembles the four terms a real chain contributes — grid accuracy nodes, interpolation error, parameter covariance and observation noise — decides which of them are independent, and produces a single covariance a deliverable can quote.

Which Terms Are Independent #

Quadrature addition assumes independence, so the first job is to decide what is independent of what. Four terms, four answers:

Figure — over what distance each term stays correlated, and what that implies.

• Observation noise is independent between points and between stages. It adds cleanly.
• Grid accuracy — the per-node accuracy values NTv2 records carry — describes how well the grid models the true distortion at that node. Between points in different cells it is largely independent; between points in the same cell it is essentially identical.
• Interpolation error is spatially correlated in exactly the same way, and for the same reason: it is a property of the cell, not of the point.
• Parameter uncertainty, where a parametric fallback ran, is common to every point transformed with that parameter set.

The pattern is that everything except observation noise is correlated over some spatial scale. That does not change the arithmetic for a single point — all four terms simply add in quadrature — but it changes what may be said about aggregates, and it is the reason the implementation below keeps the terms separate rather than summing them immediately.

Complete Runnable Implementation #

from __future__ import annotations

from dataclasses import dataclass, field

import numpy as np


@dataclass(frozen=True)
class Budget:
    """The uncertainty terms of one transformed point, kept separate.

    All values are one-sigma in metres. `independent` terms shrink over an
    aggregate; `correlated` terms do not, which is why they are not summed here.
    """

    observation: float
    interpolation: float
    grid_accuracy: float
    parameters: float = 0.0
    method: str = "grid"

    def __post_init__(self) -> None:
        for name in ("observation", "interpolation", "grid_accuracy", "parameters"):
            v = getattr(self, name)
            if v < 0.0 or not np.isfinite(v):
                raise ValueError(f"{name} must be a finite, non-negative sigma")
        if self.method not in ("grid", "parametric"):
            raise ValueError("method must be 'grid' or 'parametric'")

    @property
    def independent(self) -> float:
        """Terms that differ point to point — the relative-accuracy part."""
        return float(np.hypot(self.observation, self.interpolation))

    @property
    def correlated(self) -> float:
        """Terms shared across points — the systematic part."""
        return float(np.hypot(self.grid_accuracy, self.parameters))

    @property
    def total(self) -> float:
        """One-sigma total for the absolute position of this point."""
        return float(np.hypot(self.independent, self.correlated))

    def scaled(self, factor: float) -> float:
        """Total at a stated confidence factor (1.96 in 1-D, 2.4477 in 2-D at 95%)."""
        return self.total * factor


def aggregate(budgets: list[Budget]) -> dict[str, float]:
    """Dataset-level statement from per-point budgets.

    The independent part shrinks as 1/sqrt(n) over the set; the correlated part
    does not shrink at all. Reporting a dataset uncertainty that ignores this
    understates it by roughly sqrt(n).
    """
    if not budgets:
        raise ValueError("no budgets to aggregate")
    n = len(budgets)
    ind = np.sqrt(np.mean([b.independent ** 2 for b in budgets]) / n)
    cor = float(np.max([b.correlated for b in budgets]))
    return {
        "n": float(n),
        "independent_mean_m": float(ind),
        "correlated_m": cor,
        "dataset_sigma_m": float(np.hypot(ind, cor)),
        "worst_point_sigma_m": float(max(b.total for b in budgets)),
    }

Parameter Reference #

Worked Example #

Two points from the same batch — one served by the grid, one by the parametric fallback because it fell outside the grid extent:

Figure — a mixed batch has two populations, and one number describes neither.

import numpy as np

inside = Budget(observation=0.012, interpolation=0.004,
                grid_accuracy=0.010, method="grid")
outside = Budget(observation=0.012, interpolation=0.0,
                 grid_accuracy=0.0, parameters=0.055, method="parametric")

for name, b in (("grid", inside), ("fallback", outside)):
    print(f"{name:9s} total {b.total * 1000:5.1f} mm  "
          f"95% 2-D {b.scaled(2.4477) * 1000:5.1f} mm")
print(aggregate([inside] * 900 + [outside] * 100))
# grid      total  16.5 mm  95% 2-D  40.4 mm
# fallback  total  56.3 mm  95% 2-D 137.8 mm
# {'n': 1000.0, 'independent_mean_m': 0.00040..., 'correlated_m': 0.055,
#  'dataset_sigma_m': 0.05500..., 'worst_point_sigma_m': 0.0563...}

The aggregate line is the one worth reading twice. Averaging a thousand points drives the independent part down to under half a millimetre — and the dataset uncertainty barely moves from 55 mm, because the correlated term from the hundred fallback points does not average away. A pipeline that reported the mean of the per-point totals would have claimed about 20 mm for the dataset, understating it by nearly a factor of three.

Validation Check #

def assert_terms_present(b: Budget) -> None:
    """A budget with a missing term is the failure this guide exists to prevent."""
    if b.method == "grid":
        assert b.grid_accuracy > 0.0, "grid shift with no grid accuracy term"
        assert b.interpolation > 0.0, "grid shift with no interpolation term"
    else:
        assert b.parameters > 0.0, "parametric shift with no parameter term"
    assert b.observation > 0.0, "no observation uncertainty — was it defaulted to zero?"

Every one of those assertions fires on a zero that was never filled in, which is the realistic failure mode: an uncertainty pipeline usually breaks by omission rather than by miscalculation.

Common Mistakes #

Grid accuracy confused with interpolation error. They are different quantities from different sources. Grid accuracy is how well the modelled surface matches reality at a node, published with the grid; interpolation error is how well the value between nodes is reconstructed, measured as in measuring interpolation error against control points. Using one for both leaves a real term out of the budget.

A mixed dataset reported with one uncertainty. When some points were served by the grid and some by a fallback, their uncertainties differ by a factor of three or more. A single figure for the whole deliverable either flatters the fallback points or unfairly penalises the grid ones — and the fix is not a compromise number but a per-point figure plus a summary of both populations.

The correlated part averaged away. Dividing every term by the square root of the point count is the reflex, and it is wrong for anything shared between points. The aggregate function above shrinks only the independent part, which is the whole reason the terms are kept apart.

Recording the Budget Alongside the Coordinates #

A budget that exists only inside the pipeline is not evidence. The four terms, their sources and the method that produced each point should be written next to the coordinates — four floats and a short string per row — so a reviewer can see not just how uncertain a point is but why. That per-row record is also what makes a mixed dataset explainable: the hundred fallback points in the example above are not a defect, they are a documented consequence of the grid extent, and the record is the difference between those two readings of the same file.

Figure — five values per row that make an uncertainty statement reviewable.

Frequently Asked Questions #

Where do NTv2 per-node accuracy values come from? #

They are produced by the agency that computed the grid, usually from the residuals of the adjustment that generated it, and they vary across the surface — a well-surveyed urban area has smaller values than a sparsely controlled one. Reading them requires parsing the third and fourth floats of each node record, which is one of the few genuine reasons to write a manual grid reader rather than use the library, as when to use manual grid interpolation over pyproj discusses.

What if the grid publishes no accuracy values? #

Fall back to the operation’s declared accuracy, which is a single figure for the whole grid, and record that the per-node values were unavailable. That is coarser but honest. What should not happen is the term being dropped because there was no obvious number to put in it.

How do I decide the spatial scale over which a term is correlated? #

For interpolation error, the cell size — points in the same cell share nearly the same error, points several cells apart do not. For grid accuracy, the scale of the control network that produced the grid, which is usually tens of kilometres. For parameter uncertainty, the whole area the parameter set covers. Exact treatment requires a spatial covariance model; treating a term as fully correlated within its scale and independent beyond it is the conservative simplification and is what most specifications assume.

Does this change if the pipeline runs the operations in a different order? #

The total does not — quadrature addition is order-independent, and the coordinate stages contribute no uncertainty of their own. What changes is which stage a term is attributed to, and that matters only for the readability of the record. Keep the attribution stable across runs so two reports can be compared line by line.

Related #

• Uncertainty propagation through transformation chains — the parent topic
• Propagating covariance through a Helmert transformation
• Measuring interpolation error against control points
• Fallback routing strategies for missing grid files
• Tuning transformation thresholds for survey-grade work