CoordTransform Cadastral & High-Precision Python

Debugging Polynomial Shift Residuals in GIS

Debugging polynomial shift residuals in GIS is the targeted operation of decomposing the vector difference between observed and model-predicted control coordinates, then routing the result against a survey-grade tolerance — typically a maximum radial residual of 0.02 m for Class I cadastral control under ISO 19111 coordinate-operation accuracy reporting. This page extends the geodetic conversion math for ellipsoid-to-Cartesian pipelines, where a misaligned forward conversion silently masquerades as polynomial drift; here we isolate the residual itself as the primary diagnostic and treat its spatial structure, not its mean magnitude, as the signal worth reading.

Polynomial shift residual debugging decision flow Flowchart: build the design matrix from modelled coordinates, then test its condition number. A value above 1e12 routes to an affine order-1 fallback and re-solves; otherwise compute the per-point radial residual for every control mark. A tolerance gate then splits three ways: maximum residual within tolerance passes and logs an audit record; up to twice tolerance is a soft fail that flags and re-measures monuments; beyond twice tolerance is a hard fail that aborts to a least-squares adjustment. Build design matrix A cond(A) ≤ 1e12 ? Affine fallback re-solve, order 1 Per-point radial residual Rᵢ = √(ΔEᵢ² + ΔNᵢ²) max(R) vs tol PASS Log audit record + proceed SOFT FAIL Re-measure monuments, clean & re-solve HARD FAIL Abort → least-squares adjustment > 1e12 ≤ 1e12 ≤ tol ≤ 2·tol > 2·tol

Concept: What a Shift Residual Actually Measures

A polynomial shift model approximates non-rigid deformation across a control network by fitting low-order surfaces to the East and North components independently. For a second-order model, each predicted coordinate is a linear combination of the design-row basis [1, E, N, E2, EN, N2][1,\,E,\,N,\,E^2,\,EN,\,N^2] evaluated at the modeled coordinate. The residual at control point ii is the leftover that the surface could not absorb:

ΔEi=Eiobs−∑kckE pk ⁣(Eimod,Nimod),Ri=ΔEi2+ΔNi2 \Delta E_i = E_i^{\text{obs}} - \sum_{k} c^{E}_{k}\,p_k\!\left(E_i^{\text{mod}}, N_i^{\text{mod}}\right), \qquad R_i = \sqrt{\Delta E_i^{2} + \Delta N_i^{2}}

The radial magnitude RiR_i is the audit metric; the aggregate is the root-mean-square error RMSE=1n∑iRi2\text{RMSE} = \sqrt{\frac{1}{n}\sum_i R_i^2}. The critical insight for debugging is that residuals are diagnostic before they are aggregated. A network can pass an RMSE check while hiding a single displaced monument, and it can fail RMSE purely because the chosen polynomial order is too low to absorb real curvature. Both cases demand the per-point vector field, not the scalar summary — which is why this protocol never averages before it inspects, and why control geometry must satisfy the rank requirement of the order: a second-order fit needs at least six non-collinear points or the coefficient matrix becomes rank-deficient and the residuals are meaningless.

Residual anomalies also frequently originate upstream of the solver. Before trusting any number below, confirm every coordinate sits in one consistent projected CRS with explicit linear units — mixing US survey feet with international metres scales residuals by 3.28084 and invalidates the tolerance outright. Vertical and epoch consistency belongs to validating datum alignment with control points; resolve those first so the residual reflects model error rather than registration error.

Complete Runnable Implementation

The following self-contained function enforces explicit float64 precision, validates matrix conditioning, computes per-component residuals, and routes to a lower-order affine fallback when the design matrix is ill-conditioned. It is runnable as-is on Python 3.10+ with NumPy.

import numpy as np
from typing import Any, Dict, Tuple
import warnings


def compute_polynomial_residuals(
    observed: np.ndarray,
    modeled: np.ndarray,
    order: int = 2,
    tolerance_m: float = 0.02,
    fallback_order: int = 1,
    max_condition_number: float = 1e12,
) -> Dict[str, Any]:
    """Compute polynomial shift residuals with deterministic fallback routing.

    Residuals are the (E, N) differences between observed control coordinates
    and the polynomial-predicted coordinates. Accuracy reporting follows the
    ISO 19111:2019 coordinate-operation accuracy model. East and North are
    solved independently in explicit IEEE 754 float64.
    """
    obs = np.asarray(observed, dtype=np.float64)
    mod = np.asarray(modeled, dtype=np.float64)

    if obs.shape != mod.shape or obs.ndim != 2 or obs.shape[1] != 2:
        raise ValueError("Input arrays must be Nx2 (E, N) matrices.")

    n_points = obs.shape[0]

    def _design(coords: np.ndarray, poly_order: int) -> np.ndarray:
        E, N = coords[:, 0], coords[:, 1]
        if poly_order == 1:  # affine basis (ISO 19111 polynomial order 1)
            return np.column_stack([np.ones(n_points), E, N])
        if poly_order == 2:  # quadratic basis: 1, E, N, E^2, EN, N^2
            return np.column_stack(
                [np.ones(n_points), E, N, E ** 2, E * N, N ** 2]
            )
        raise NotImplementedError("Only 1st and 2nd order shifts are supported.")

    def _solve(A: np.ndarray, b: np.ndarray) -> Tuple[np.ndarray, float]:
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            coeffs, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
        return coeffs, float(np.linalg.cond(A))

    A = _design(mod, order)
    try:
        coeffs_e, cond_e = _solve(A, obs[:, 0])
        coeffs_n, cond_n = _solve(A, obs[:, 1])
        if max(cond_e, cond_n) > max_condition_number:
            warnings.warn(
                f"Condition number {max(cond_e, cond_n):.2e} exceeds limit; "
                f"routing to order {fallback_order}."
            )
            return compute_polynomial_residuals(
                obs, mod, order=fallback_order, tolerance_m=tolerance_m
            )
    except np.linalg.LinAlgError:
        warnings.warn("Singular matrix; routing to order 1 (affine).")
        return compute_polynomial_residuals(
            obs, mod, order=1, tolerance_m=tolerance_m
        )

    delta_e = obs[:, 0] - A @ coeffs_e
    delta_n = obs[:, 1] - A @ coeffs_n
    radial = np.sqrt(delta_e ** 2 + delta_n ** 2)

    max_residual = float(np.max(radial))
    rmse = float(np.sqrt(np.mean(radial ** 2)))

    return {
        "coefficients": {"E": coeffs_e, "N": coeffs_n},
        "residuals": {"dE": delta_e, "dN": delta_n, "R": radial},
        "diagnostics": {
            "order": order,
            "rank": int(np.linalg.matrix_rank(A)),
            "condition_number": max(cond_e, cond_n),
            "max_residual_m": max_residual,
            "rmse_m": rmse,
            "passes_tolerance": max_residual <= tolerance_m,
        },
    }

The fallback recursion is the same defensive pattern used when implementing affine transformations for local grids: an order-1 model is the deterministic safe harbour when the higher-order coefficient matrix cannot be solved stably.

Parameter and return reference

Name Type Units Valid range / meaning
observed np.ndarray Nx2 metres Field-surveyed control (E, N); N ≥ 6 for order 2
modeled np.ndarray Nx2 metres Source-frame coordinates fed to the polynomial
order int — 1 (affine) or 2 (quadratic)
tolerance_m float metres Max radial residual gate; 0.02 Class I, 0.05 Class II
max_condition_number float — Triggers fallback above 1e12
diagnostics.max_residual_m float metres Worst per-point radial residual
diagnostics.rmse_m float metres Network root-mean-square residual
diagnostics.passes_tolerance bool — True when max_residual_m ≤ tolerance_m

Minimal Worked Example

Generate a network whose observed coordinates follow an exact second-order surface of the modeled coordinates, then recover it and read the diagnostics:

import numpy as np

mod = np.array(
    [[0.0, 0.0], [1000.0, 0.0], [0.0, 1000.0], [1000.0, 1000.0],
     [500.0, 250.0], [250.0, 750.0], [750.0, 500.0]],
    dtype=np.float64,
)
E, N = mod[:, 0], mod[:, 1]
# Known true quadratic shift surface (sub-metre coefficients):
obs_e = 12.000 + 1.0000020 * E + 0.0000010 * N + 1e-9 * E * N
obs_n = -8.000 + 0.9999980 * N + 0.0000015 * E + 1e-9 * E ** 2
obs = np.column_stack([obs_e, obs_n])

result = compute_polynomial_residuals(obs, mod, order=2, tolerance_m=0.02)
print(round(result["diagnostics"]["rmse_m"], 6),
      result["diagnostics"]["passes_tolerance"])
# -> 0.0 True   (residuals at machine precision; the surface is recovered exactly)

Because the observed coordinates were produced by a genuine quadratic surface, the least-squares solver reconstructs the coefficients and the radial residuals collapse to the float64 noise floor — max_residual_m on the order of 1e-10 m. Perturb a single monument (for example obs[4] += [0.08, 0.0]) and the least-squares fit smears that 0.08 m displacement across the surface, lifting max_residual_m just past the 0.02 m gate (about 0.021 m) so passes_tolerance flips to False: the classic signature of a disturbed control mark, which the per-point vector field localises far more sharply than the aggregate RMSE.

Validation Check

Gate the diagnostics with an explicit survey-grade assertion before any coordinate is written downstream:

diag = result["diagnostics"]
assert diag["rank"] >= 6, "Rank deficient: control geometry cannot support order 2."
assert diag["passes_tolerance"], (
    f"Max residual {diag['max_residual_m']:.4f} m exceeds tolerance."
)

For a network that genuinely sits at the threshold, route the three lanes deterministically: pass when max(R) ≤ tolerance_m; soft fail when tolerance_m < max(R) ≤ 2 × tolerance_m, where you flag and re-measure the offending monuments and re-solve on the cleaned set; and hard fail when max(R) > 2 × tolerance_m or the rank check fails, where you abort and hand the network to a rigorous least-squares adjustment for control networks. Whichever lane fires, persist an audit record carrying the order, condition number, RMSE, pass/fail status, and the source/target EPSG codes so the result is reproducible and legally defensible.

Reading the Residual Vector Field

Plot the per-point dE/dN vectors with graduated symbology before aggregating — the orientation pattern names the defect:

Four residual vector-field signatures and their causes Four panels of control-point residual arrows. Top left: every arrow points the same way, a uniform bias from a datum translation, fixed by removing a rigid shift before fitting. Top right: arrows circulate tangentially about a pivot, a rotation from azimuth or grid-convergence mismatch, fixed by checking the central meridian and scale. Bottom left: arrows radiate outward from a centroid, a global scale-factor error, fixed by correcting the linear scale. Bottom right: arrows reverse direction across the network, real higher-order distortion or epoch drift, absorbed by raising to a second-order surface. Uniform bias → datum translation Fix: remove rigid shift before fitting Rotation → azimuth / convergence Fix: check central meridian & scale Radial → global scale error Fix: correct linear scale factor Reversing curvature → distortion Fix: raise to order 2 / check epoch
  1. Uniform directional bias. A consistent vector orientation across the whole network is a false-origin offset or datum translation; remove it with a rigid shift before fitting, which is what the surrounding polynomial shift algorithms for regional adjustments assume as a precondition.
  2. Rotational pattern. Tangential vectors circulating a pivot indicate azimuth or grid-convergence mismatch; verify the projection central meridian and scale factor against the projection math fundamentals for cadastral surveys.
  3. Radial expansion or contraction. Vectors diverging from a centroid signal a global scale-factor discrepancy; correct the linear scale parameter rather than raising the polynomial order.
  4. Reversing curvature. Residuals that flip direction across the polygon are real higher-order distortion or epoch drift; a second-order surface usually absorbs it, and if clustered residuals persist, inspect monument metadata, which is the failure mode covered by error-distribution modeling in Python.

Common Mistakes

Aggregating to RMSE before inspecting the vector field
A single displaced monument can pass an RMSE check while corrupting every parcel near it, and an under-ordered model can fail RMSE despite a clean network. Always render the per-point `dE`/`dN` vectors first; the orientation pattern, not the scalar mean, tells you whether the defect is registration, geometry, or model order.
Letting float32 promotion erode sub-centimetre residuals
Reading control coordinates from a source that defaults to `float32` (many raster and shapefile loaders do) caps precision near 1e-1 m on metre-scale eastings, swamping a 0.02 m tolerance. Cast every input array to `np.float64` explicitly — as the implementation does on entry — and never rely on implicit NumPy promotion mid-pipeline.
Fitting a second-order surface to collinear or clustered control
Six points are necessary but not sufficient for an order-2 fit: if they are near-collinear, the design matrix is rank-deficient, the condition number explodes, and the residuals are numerically meaningless. Check `diagnostics.rank` and `condition_number`, and let the affine fallback engage rather than reporting a spuriously small residual from an unstable solve.