CoordTransform Cadastral & High-Precision Python

Choosing Polynomial Degree for Regional Fits

Selecting the degree of a polynomial coordinate-shift surface is the decision that quietly sets the accuracy ceiling of every regional adjustment, and it belongs squarely under polynomial shift algorithms for regional adjustments within the wider Algorithmic Math & Geodetic Workflows reference. Push the degree too low and the surface cannot follow the regional distortion; push it too high and it starts modelling the noise in your control points, so the choice is never “as high as the data allows” but “as high as held-out data can justify”.

Bias-variance trade-off against polynomial degree Two curves plotted against increasing polynomial degree on the horizontal axis and error on the vertical axis. The training RMSE curve falls monotonically toward zero as degree rises. The validation RMSE curve falls to a minimum at a moderate degree and then rises sharply as overfitting sets in. A dashed vertical line marks the optimum degree at the validation minimum. RMSE degree training validation optimum

Figure — training error keeps falling; validation error marks the usable degree.

The Trade-off You Are Actually Balancing

A regional shift surface expresses the residual displacement Δx\Delta x between two frames as a bivariate polynomial in the plan coordinates u=E−E0u = E - E_0 and v=N−N0v = N - N_0, centred and scaled about the control extent to keep the design matrix well-conditioned:

Δx(u,v)=∑i=0d∑j=0d−iaij uivj\Delta x(u, v) = \sum_{i=0}^{d}\sum_{j=0}^{d-i} a_{ij}\, u^{i} v^{j}

Raising the degree dd adds coefficients and always lowers the training residual, because more free parameters can bend the surface through more control points. What it cannot lower is the error on coordinates the fit never saw. The honest scoring metric is the cross-validated root-mean-square error over a held-out set of mm points,

RMSEval=1m∑k=1m(Δx^k−Δxk)2\text{RMSE}_{\text{val}} = \sqrt{\frac{1}{m}\sum_{k=1}^{m}\left(\Delta \hat{x}_k - \Delta x_k\right)^2}

which first falls as the surface captures genuine regional trend, then rises once the extra terms start fitting measurement noise. A second, independent penalty grows with degree: the condition number of the normal-equation matrix climbs steeply, so the coefficient solution becomes numerically fragile long before the mathematics forbids it. Watching both validation RMSE and conditioning is the same defensive posture used when debugging polynomial shift residuals in GIS, where an over-fitted surface hides its error inside the control set.

Complete Runnable Implementation

The function fits every degree from one up to a ceiling, scores each on a held-out validation split, records the condition number, and returns the degree with the lowest validation RMSE — refusing any degree whose conditioning has already exploded. Coordinates are centred and scaled before the Vandermonde terms are built so the comparison across degrees is fair.

from __future__ import annotations

from dataclasses import dataclass

import numpy as np


@dataclass(frozen=True)
class DegreeScore:
    degree: int
    train_rmse: float
    val_rmse: float
    condition_number: float
    n_terms: int


@dataclass(frozen=True)
class DegreeSelection:
    best_degree: int
    scores: list[DegreeScore]


def _design_matrix(u: np.ndarray, v: np.ndarray, degree: int) -> np.ndarray:
    """Bivariate polynomial terms u^i v^j with i + j <= degree."""
    cols = [u ** i * v ** j
            for i in range(degree + 1)
            for j in range(degree + 1 - i)]
    return np.column_stack(cols)


def choose_polynomial_degree(
    east: np.ndarray,
    north: np.ndarray,
    shift: np.ndarray,
    max_degree: int = 4,
    val_fraction: float = 0.3,
    max_condition: float = 1e10,
    seed: int = 12345,
) -> DegreeSelection:
    """Pick the shift-surface degree that minimises held-out validation RMSE.

    Args:
        east, north: control-point plan coordinates (metres).
        shift: the coordinate shift to model at each control point (metres).
        max_degree: highest polynomial degree to trial.
        val_fraction: fraction of points reserved for validation.
        max_condition: degrees whose design matrix exceeds this conditioning
            are rejected regardless of their validation RMSE.
        seed: fixed RNG seed so the train/validation split is reproducible.

    Returns:
        A DegreeSelection with the winning degree and per-degree scores.
    """
    e = np.asarray(east, dtype=np.float64)
    n = np.asarray(north, dtype=np.float64)
    s = np.asarray(shift, dtype=np.float64)
    if not (e.shape == n.shape == s.shape):
        raise ValueError("east, north and shift must share one length")

    # Centre and scale to the control extent (numerical conditioning, not geodesy).
    u = (e - e.mean()) / (float(np.ptp(e)) or 1.0)
    v = (n - n.mean()) / (float(np.ptp(n)) or 1.0)

    rng = np.random.default_rng(seed)
    idx = rng.permutation(u.size)
    n_val = max(1, int(round(val_fraction * u.size)))
    val_idx, train_idx = idx[:n_val], idx[n_val:]

    scores: list[DegreeScore] = []
    best_degree, best_val = 1, np.inf
    for degree in range(1, max_degree + 1):
        a_train = _design_matrix(u[train_idx], v[train_idx], degree)
        if a_train.shape[0] < a_train.shape[1]:
            break  # more coefficients than training points: not identifiable
        coef, *_ = np.linalg.lstsq(a_train, s[train_idx], rcond=None)
        cond = float(np.linalg.cond(a_train))

        train_pred = a_train @ coef
        val_pred = _design_matrix(u[val_idx], v[val_idx], degree) @ coef
        train_rmse = float(np.sqrt(np.mean((train_pred - s[train_idx]) ** 2)))
        val_rmse = float(np.sqrt(np.mean((val_pred - s[val_idx]) ** 2)))
        scores.append(DegreeScore(degree, train_rmse, val_rmse, cond, coef.size))

        if cond <= max_condition and val_rmse < best_val:
            best_val, best_degree = val_rmse, degree

    return DegreeSelection(best_degree=best_degree, scores=scores)

Parameter Reference

Name Type Units Valid range Cadastral significance
east, north np.ndarray metres control extent plan coordinates of the control points
shift np.ndarray metres finite the displacement being modelled
max_degree int — 1…≈5 ceiling on trialled degrees
val_fraction float — 0<x<1 share reserved for validation
max_condition float — ~1e8…1e12 conditioning veto on unstable fits
best_degree int — 1…max degree with lowest validation RMSE
condition_number float — ≥ 1 numerical stability of that degree’s normal equations

Worked Example

A synthetic quadratic distortion with light noise recovers degree two, not the highest degree offered — the extra terms buy nothing on held-out points.

import numpy as np

rng = np.random.default_rng(7)
e = rng.uniform(0, 1000, 40)
n = rng.uniform(0, 1000, 40)
true = 0.02 + 3e-5 * e - 4e-5 * n + 1e-7 * e * n     # quadratic-ish trend
obs = true + rng.normal(0, 0.003, e.size)            # 3 mm control noise
sel = choose_polynomial_degree(e, n, obs, max_degree=4, seed=7)
print(sel.best_degree)
for sc in sel.scores:
    print(sc.degree, round(sc.train_rmse, 4), round(sc.val_rmse, 4))
# 2
# 1 0.0096 0.0066
# 2 0.002 0.0036
# 3 0.0016 0.0038
# 4 0.0012 0.0043

Validation Check

Confirm the winner beats the linear baseline on validation while its training error is not driven to zero — a training RMSE near machine precision is a symptom, not a success.

best = next(s for s in sel.scores if s.degree == sel.best_degree)
linear = next(s for s in sel.scores if s.degree == 1)
assert best.val_rmse <= linear.val_rmse, "chosen degree lost to the linear baseline"
assert best.train_rmse > 1e-6, "training residual collapsed — surface is memorising noise"
assert best.condition_number < 1e10, "selected degree is numerically unstable"

Common Mistakes

Raising the degree until the training residual is zero
With enough terms a polynomial can pass through every control point exactly, driving the training RMSE to zero — and the surface then oscillates wildly between points. Zero training residual is the signature of overfitting, not accuracy. Score on data the fit never touched.
Fitting and scoring on the same points
Evaluating the surface on the very points used to solve its coefficients always flatters higher degrees, because they have more freedom to hug those points. Reserve a held-out split (or use k-fold cross-validation) and compare degrees on that untouched set so the reported error reflects real predictive skill.
Ignoring the condition number
The normal-equation matrix conditioning grows sharply with degree, especially on un-centred metre-scale coordinates. A high condition number means the coefficients are dominated by round-off long before the validation curve says stop. Centre and scale the coordinates, and veto any degree whose conditioning has exploded.
Extrapolating past the control extent
A polynomial fit is only trustworthy inside the convex hull of the control it was solved from. Outside that footprint the higher-order terms diverge rapidly, so a degree that scores well on interior validation can still be badly wrong at the edge of a parcel. Constrain use to the control extent, or fall back to a global parametric transformation there.