Error Distribution Modeling in Python for Cadastral & High-Precision Coordinate Transformation

Error distribution modeling is the sub-task within Algorithmic Math & Geodetic Workflows that converts the residuals of a coordinate transformation into a defensible, propagated uncertainty budget. When survey observations are migrated across datum realizations, projection zones, or legacy grid systems, the residual deviations left behind are not computational noise to be discarded; they are quantifiable departures from ground truth that must be bounded, propagated through every transformation stage, and recorded for audit. This guide implements that workflow in deterministic Python: a weighted least-squares estimator with full covariance propagation, an error-ellipse decomposition, and the verification and logging that survey-grade work demands under ISO 19111 coordinate reference system metadata expectations.

Survey-grade pipelines must rigorously separate systematic bias from random noise before any distribution can be modeled. Systematic components include grid scale distortion, datum shift offsets, and epoch-dependent crustal motion; random components encompass instrumental variance, atmospheric refraction residuals, and control point measurement uncertainty. Deterministic behavior is enforced through strict np.float64 arithmetic, explicit propagation of covariance matrices, and automated routing to a lower-order model when a control network degrades. Localized networks are typically corrected with affine transformations for local grids, whose linear form yields analytically tractable error distributions; broader jurisdictional extents that exhibit non-uniform deformation instead call for polynomial shift algorithms for regional adjustments. Both families feed the same stochastic machinery described below.

The deterministic error-distribution pipeline: covariance enters at the top, the condition-number gate guards the inversion (routing degenerate systems to the affine fallback), and every path terminates at the tolerance gate that decides commit versus review.

The mathematical prerequisite: propagation of variance

Every coordinate this pipeline emits inherits its uncertainty from the observations that produced it, mapped through the transformation’s sensitivity. The governing relation is the general law of propagation of variance: for an output vector related to inputs by a (locally linear) map with Jacobian $\mathbf{J}$ , the output covariance is

\mathbf{\Sigma}_{x} = \mathbf{J}\,\mathbf{\Sigma}_{\ell}\,\mathbf{J}^{\mathsf{T}}

where $\mathbf{\Sigma}_{\ell}$ is the covariance of the observations. For a parametric transformation fitted by weighted least squares, the design matrix $\mathbf{A}$ plays the role of the Jacobian. With the weight matrix defined as the inverse observation covariance, $\mathbf{W} = \mathbf{\Sigma}_{\ell}^{-1}$ , the normal matrix and the estimated parameters are

\mathbf{N} = \mathbf{A}^{\mathsf{T}}\mathbf{W}\mathbf{A}, \qquad \hat{\mathbf{x}} = \mathbf{N}^{-1}\mathbf{A}^{\mathsf{T}}\mathbf{W}\,\boldsymbol{\ell}

The dispersion of those parameters follows directly from the propagation law, scaled by the a-posteriori variance factor $\hat{\sigma}_0^2$ that absorbs any mismatch between the assumed and the realized observation precision:

\mathbf{\Sigma}_{\hat{x}} = \hat{\sigma}_0^2\,\mathbf{N}^{-1}, \qquad \hat{\sigma}_0^2 = \frac{\mathbf{v}^{\mathsf{T}}\mathbf{W}\mathbf{v}}{n-u}

Here $\mathbf{v}$ is the residual vector, $n$ the number of observation equations, and $u$ the number of estimated parameters, so $n-u$ is the redundancy (degrees of freedom). A variance factor near $1.0$ confirms the stochastic model matches reality; values far from unity signal mis-weighted observations or unmodeled systematic shifts, exactly the diagnostic that least squares adjustment for control networks formalizes with a global chi-square test. Because $\hat{\sigma}_0^2$ and $\mathbf{N}^{-1}$ are the only quantities needed to reconstruct the full positional uncertainty of any derived point, they are the two outputs this page treats as non-negotiable.

Step-by-step implementation

Step 1 — Build the polynomial design matrix

The design matrix encodes how each transformation parameter influences each observation. Keeping its construction in a single, order-parameterised function means the affine case ( $\text{order}=1$ ) and the regional polynomial case share one code path, so a fallback never silently changes the model’s geometry.

import numpy as np
from numpy.typing import NDArray


def construct_design(coords: NDArray[np.float64], order: int) -> NDArray[np.float64]:
    """Assemble the design matrix of monomials x**i * y**j with i + j <= order.

    Includes the constant term (i = j = 0) so the translation component of the
    transformation is estimated, satisfying the affine model required by
    ISO 19111 coordinate operation parameterisation.
    """
    coords = np.asarray(coords, dtype=np.float64)  # explicit double precision
    if coords.ndim != 2 or coords.shape[1] != 2:
        raise ValueError("coords must be an (n, 2) array of planar E, N values.")

    terms: list[NDArray[np.float64]] = []
    for i in range(order + 1):
        for j in range(order + 1 - i):
            terms.append((coords[:, 0] ** i * coords[:, 1] ** j).reshape(-1, 1))
    return np.hstack(terms)

Step 2 — Estimate parameters and propagate covariance

This is the core estimator. It enforces double precision, monitors the condition number of the normal matrix before inversion, and routes to an affine fit when a higher-order system is degenerate, recording which model actually produced the answer. Every transformed coordinate that leaves this function carries the parameter covariance needed to reconstruct its own error ellipse.

import warnings
from typing import Any
import numpy as np
from numpy.typing import NDArray
from numpy.linalg import LinAlgError, cond, inv


def construct_design(coords: NDArray[np.float64], order: int) -> NDArray[np.float64]:
    coords = np.asarray(coords, dtype=np.float64)
    terms: list[NDArray[np.float64]] = []
    for i in range(order + 1):
        for j in range(order + 1 - i):
            terms.append((coords[:, 0] ** i * coords[:, 1] ** j).reshape(-1, 1))
    return np.hstack(terms)


def propagate_transformation_error(
    source_pts: NDArray[np.float64],
    target_pts: NDArray[np.float64],
    obs_cov: NDArray[np.float64],
    poly_order: int = 1,
    max_condition: float = 1e12,
    fallback_enabled: bool = True,
) -> dict[str, Any]:
    """Fit a transformation and propagate observation covariance to its parameters.

    Implements the propagation law Sigma_x = sigma0^2 * N^-1 with explicit
    float64 arithmetic and deterministic fallback routing, emitting an
    ISO 19111-style uncertainty record for downstream audit.
    """
    X_src = np.asarray(source_pts, dtype=np.float64)
    X_tgt = np.asarray(target_pts, dtype=np.float64)
    C_obs = np.asarray(obs_cov, dtype=np.float64)

    if X_src.shape[0] != X_tgt.shape[0]:
        raise ValueError("Control point arrays must be congruent.")
    if X_src.shape[0] < 3:
        raise ValueError("Minimum three control points required for rank sufficiency.")

    if C_obs.shape != (2, 2):
        raise ValueError("obs_cov must be the homogeneous 2x2 per-point covariance.")
    # Homogeneous network: collapse the per-point covariance to a scalar point
    # weight w = 1 / mean-axis-variance, then W = w * I_n is the n x n weight
    # matrix W = Sigma_l^-1 used in the weighted normal equations.
    point_var = float(np.trace(C_obs) / 2.0)
    n_pts = X_src.shape[0]
    W = np.eye(n_pts, dtype=np.float64) / point_var

    def _fit(order: int) -> tuple[Any, ...]:
        A = construct_design(X_src, order)
        N = A.T @ W @ A                        # weighted normal equations
        U = A.T @ W @ X_tgt
        c_num = float(cond(N))
        if c_num > max_condition:
            raise LinAlgError(f"Condition number {c_num:.2e} exceeds tolerance.")
        params = inv(N) @ U
        residuals = X_tgt - A @ params         # n x 2 residual field
        dof = (n_pts - A.shape[1]) * 2         # 2 axes observed per control point
        # a-posteriori variance factor sigma0^2 = v^T W v / (n - u), dimensionless
        sigma0_sq = float(np.sum((W @ residuals) * residuals) / dof) if dof > 0 else 1.0
        return params, residuals, N, c_num, dof, sigma0_sq

    model_used = "polynomial"
    try:
        params, residuals, N, c_num, dof, sigma0_sq = _fit(poly_order)
    except LinAlgError:
        if not fallback_enabled:
            raise RuntimeError("Matrix inversion failed and fallback disabled.")
        warnings.warn("Higher-order system degenerate. Routing to affine fallback.", RuntimeWarning)
        model_used = "affine_fallback"
        params, residuals, N, _c, dof, sigma0_sq = _fit(1)
        c_num = float("inf")

    C_params = inv(N) * sigma0_sq      # Sigma_x_hat = sigma0^2 * N^-1

    return {
        "transformation_params": params,
        "parameter_covariance": C_params,
        "post_fit_variance_factor": sigma0_sq,
        "condition_number": c_num,
        "degrees_of_freedom": dof,
        "model_used": model_used,
    }

Step 3 — Decompose covariance into a positional error ellipse

A 2×2 covariance submatrix is not human-readable; its eigen-decomposition is. The eigenvalues give the squared semi-axis lengths and the dominant eigenvector gives the orientation. Scaling by the chi-square quantile $k = \chi^2_{2,\,0.95} = 5.991$ converts the standard ellipse into a 95% confidence region in the plane:

a = \sqrt{k\,\lambda_{1}}, \qquad b = \sqrt{k\,\lambda_{2}}, \qquad \theta = \operatorname{atan2}(u_{1y}, u_{1x})

Using the symmetric solver numpy.linalg.eigh guarantees real eigenvalues and orthonormal eigenvectors even when accumulated floating-point drift would push a general solver into complex territory. Rendering these regions for review is covered in visualizing geodetic error ellipses with Matplotlib.

import numpy as np
from numpy.typing import NDArray


def covariance_to_error_ellipse(
    cov_2x2: NDArray[np.float64],
    chi2_quantile: float = 5.991,  # chi-square, 2 dof, 95% confidence
) -> dict[str, float]:
    """Convert a 2x2 positional covariance into a confidence error ellipse.

    Returns semi-major (a) and semi-minor (b) axes in the covariance's length
    units and the orientation theta in radians, measured from the +E axis.
    """
    cov = np.asarray(cov_2x2, dtype=np.float64)
    if cov.shape != (2, 2):
        raise ValueError("Expected a 2x2 positional covariance submatrix.")
    if not np.allclose(cov, cov.T, atol=1e-12):
        raise ValueError("Covariance must be symmetric.")

    eigvals, eigvecs = np.linalg.eigh(cov)      # ascending eigenvalues
    if np.any(eigvals < -1e-12):
        raise ValueError("Covariance is not positive semi-definite.")
    eigvals = np.clip(eigvals, 0.0, None)       # damp tiny negative drift

    order = np.argsort(eigvals)[::-1]           # major axis first
    lam_major, lam_minor = eigvals[order]
    major_vec = eigvecs[:, order[0]]

    return {
        "semi_major": float(np.sqrt(chi2_quantile * lam_major)),
        "semi_minor": float(np.sqrt(chi2_quantile * lam_minor)),
        "theta_rad": float(np.arctan2(major_vec[1], major_vec[0])),
    }

Parameter and return-value reference

The estimator’s contract is summarised below. Units follow the projected working frame (metres for planar coordinates, metres² for the variance terms of the observation covariance).

Name	Direction	Type	Units	Valid range	Cadastral significance
`source_pts`	in	`NDArray[float64]` `(n, 2)`	metres	n ≥ 3	Control points in the source frame; rank source of the fit
`target_pts`	in	`NDArray[float64]` `(n, 2)`	metres	congruent with source	Known coordinates in the target datum/projection
`obs_cov`	in	`NDArray[float64]` `(2, 2)`	metres²	symmetric, positive-definite	Stochastic model of observation precision
`poly_order`	in	`int`	—	1–3	1 = affine; 2–3 = regional surface fit
`max_condition`	in	`float`	—	1e6 – 1e14	Upper bound on `cond(N)` before fallback fires
`transformation_params`	out	`NDArray[float64]`	mixed	—	Estimated coefficients per axis
`parameter_covariance`	out	`NDArray[float64]`	parameter²	—	Drives every downstream error ellipse
`post_fit_variance_factor`	out	`float`	—	ideally ≈ 1.0	Goodness-of-fit of the stochastic model
`condition_number`	out	`float`	—	≤ `max_condition`	Numerical reliability of the inversion
`model_used`	out	`str`	—	`polynomial` / `affine_fallback`	Audit flag recording which model produced the result

Worked example with real-world coordinates

Consider a small cadastral network being re-fitted from a legacy local grid onto a national projection — for example transforming working coordinates toward the British National Grid (EPSG:27700) plane. Four control monuments are observed with an isotropic precision of $\sigma = 0.02\,\text{m}$ , giving an observation covariance of $\operatorname{diag}(4\times10^{-4}, 4\times10^{-4})\,\text{m}^2$ . The script below fits the affine model, reports the variance factor, and decomposes the covariance of the first estimated coefficient into an error ellipse.

import numpy as np


def construct_design(coords, order):
    coords = np.asarray(coords, dtype=np.float64)
    terms = []
    for i in range(order + 1):
        for j in range(order + 1 - i):
            terms.append((coords[:, 0] ** i * coords[:, 1] ** j).reshape(-1, 1))
    return np.hstack(terms)


# Four control monuments: source local grid -> target EPSG:27700 metres.
src = np.array([[0.0, 0.0], [100.0, 0.0], [0.0, 100.0], [100.0, 100.0]], dtype=np.float64)
tgt = np.array(
    [[529100.00, 181200.00],
     [529200.01, 181199.98],
     [529100.02, 181300.01],
     [529199.99, 181299.97]],
    dtype=np.float64,
)
obs_cov = np.diag([4e-4, 4e-4]).astype(np.float64)   # sigma = 0.02 m per axis

point_var = float(np.trace(obs_cov) / 2.0)         # 4e-4 m^2 per axis
W = np.eye(src.shape[0]) / point_var               # n x n point weight matrix
A = construct_design(src, 1)
N = A.T @ W @ A
params = np.linalg.inv(N) @ (A.T @ W @ tgt)
residuals = tgt - A @ params
dof = (src.shape[0] - A.shape[1]) * 2
sigma0_sq = float(np.sum((W @ residuals) * residuals) / dof)

print(f"degrees of freedom : {dof}")
print(f"variance factor    : {sigma0_sq:.3f}")
print(f"max |residual| (m) : {np.max(np.abs(residuals)):.4f}")
# Expected (approx): dof = 2, variance factor ~ 0.6, max residual ~ 0.01 m

The peak residual of roughly $0.01\,\text{m}$ sits inside the declared observation precision, and a dimensionless variance factor of order $0.6$ — below unity because only two degrees of freedom remain — confirms the affine model is consistent with the declared $\sigma = 0.02\,\text{m}$ stochastic model rather than straining against it. If an independent check monument were withheld and re-predicted, its transformed position should fall inside the 95% ellipse produced by Step 3 — the standard acceptance test before a coordinate is committed.

Verification and residual analysis

Convergence and a plausible variance factor are necessary but not sufficient. Production routing demands an explicit tolerance gate and a structured, serialisable record so that the decision to accept or reject a coordinate is reproducible during a cadastral audit. The snippet below computes per-point planar residuals, the network RMSE, and emits a JSON-ready log entry keyed to the survey-grade tolerance. The same threshold discipline is generalised in tuning transformation thresholds for survey-grade work.

import json
import logging
from typing import Any
import numpy as np
from numpy.typing import NDArray

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger("error_distribution")


def audit_residuals(
    residuals: NDArray[np.float64],
    tolerance_m: float = 0.05,   # statutory planar tolerance
) -> dict[str, Any]:
    """Compare planar residuals against tolerance and emit an audit record."""
    res = np.asarray(residuals, dtype=np.float64)
    point_mag = np.sqrt(np.sum(res ** 2, axis=1))         # per-point magnitude
    rmse = float(np.sqrt(np.mean(point_mag ** 2)))
    worst_idx = int(np.argmax(point_mag))

    record = {
        "rmse_m": round(rmse, 6),
        "max_residual_m": round(float(point_mag[worst_idx]), 6),
        "worst_point_index": worst_idx,
        "tolerance_m": tolerance_m,
        "within_tolerance": bool(point_mag.max() <= tolerance_m),
    }
    logger.info("residual audit: %s", json.dumps(record))
    return record


# Minimal demonstration with a residual array (E, N) per control point.
demo_res = np.array([[0.011, -0.018], [0.020, 0.009], [-0.014, 0.021]], dtype=np.float64)
report = audit_residuals(demo_res, tolerance_m=0.05)
assert report["within_tolerance"], "Residuals exceed survey-grade tolerance."

The assert is the gate: a residual beyond the statutory band stops the pipeline rather than silently writing a non-compliant coordinate to the cadastral store. Upstream, geometric soundness of the inputs can be verified first with a Python script for geodetic inverse problem solving, which checks baseline distances and azimuths before the stochastic model ever runs.

Troubleshooting and gotchas

Singular or near-singular normal matrix. Tightly clustered or collinear control points inflate cond(N) and corrupt the inverse. The estimator’s max_condition gate catches this and routes to the affine model; if even the affine system is rank-deficient, add geometrically independent control rather than lowering the threshold. Always inspect the returned model_used flag — a silent demotion to affine_fallback changes the meaning of the parameter covariance.

Single- vs double-precision promotion. Passing float32 arrays (common from some I/O paths) silently degrades the variance factor and can flip small eigenvalues negative during decomposition. Every entry point here calls np.asarray(..., dtype=np.float64); never remove those casts.

Negative eigenvalues in the ellipse step. A covariance that is positive-definite in theory can show eigenvalues of order $-10^{-15}$ after repeated inversions. covariance_to_error_ellipse clips drift within tolerance but raises on genuinely indefinite input — do not widen the atol to suppress the error, since a truly negative eigenvalue means the upstream fit is unreliable.

Mis-scaled observation covariance. If obs_cov is supplied in the wrong units (e.g. millimetres² instead of metres²), the fit still converges but the variance factor lands orders of magnitude from $1.0$ . Treat any $\hat{\sigma}_0^2$ outside roughly $0.1$ – $10$ as a stochastic-model error, not a measurement outlier.

Degrees of freedom at zero. With exactly the minimum control (three points for an affine fit) the redundancy is zero, the variance factor is undefined, and the code falls back to $\hat{\sigma}_0^2 = 1.0$ . The covariance is then formally optimistic; add redundant control before treating the ellipse as defensible.

Algorithmic Math & Geodetic Workflows — parent guide to the deterministic transformation pipeline
Least squares adjustment for control networks — the global chi-square test and Gauss-Markov model behind the variance factor
Polynomial shift algorithms for regional adjustments — higher-order surface fitting for non-uniform distortion
Visualizing geodetic error ellipses with Matplotlib — rendering the confidence regions produced here
Python script for geodetic inverse problem solving — upstream geometric validation of control inputs

Error Distribution Modeling in Python for Cadastral & High-Precision Coordinate Transformation

The mathematical prerequisite: propagation of variance #

Step-by-step implementation #

Step 1 — Build the polynomial design matrix #

Step 2 — Estimate parameters and propagate covariance #

Step 3 — Decompose covariance into a positional error ellipse #

Parameter and return-value reference #

Worked example with real-world coordinates #

Verification and residual analysis #

Troubleshooting and gotchas #

Related guides #