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:

Fit and hold-out error against polynomial degree Two traces against polynomial degree from one to six. The fit RMSE falls monotonically from 0.082 to 0.016 metres, because a higher-degree polynomial can always fit the points it was given more closely. The hold-out RMSE falls to a minimum of 0.034 metres at degree three and then rises steeply to 0.121 metres at degree six. Degree three is the choice; anything beyond it is fitting the noise in the control set. 0.000 0.050 0.100 1 2 3 4 5 6 polynomial degree m fit RMSE (m) hold-out RMSE (m)

Figure — fit error always falls with degree; hold-out error is the one that tells the truth.

Δ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.

Term count and control requirement by polynomial degree Four degrees with their term count per coordinate, the minimum number of control points and a practical minimum. Degree one has three terms and needs three points, ten in practice. Degree two has six terms and needs six, twenty in practice. Degree three has ten terms and needs ten, thirty-five in practice. Degree four has fifteen terms, needs fifteen and fifty in practice, at which point the honest question is whether a grid would serve better than a polynomial. Terms Minimum In practice degree 1 3 3 10 degree 2 6 6 20 degree 3 10 10 35 degree 4 15 15 50 — use a grid

Figure — term count grows quadratically with degree, and so does the control set you need.

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.

Frequently Asked Questions

Is there a rule of thumb for the degree?

Start at degree two and only go higher if hold-out error genuinely improves; in practice degree three is the ceiling for most regional work. The reason is not aesthetic: each degree adds terms quadratically, and the control needed to constrain them grows with it, so a degree-four surface fitted to twenty points is describing the noise in those twenty points in considerable detail.

What happens outside the control extent?

A polynomial diverges, and it does so fast — that is the defining behaviour of the model class. A degree-three surface that fits to two centimetres inside its control hull can be metres out fifty kilometres beyond it, with no warning in any statistic computed from the fit. Clamp the model to its hull and raise rather than extrapolate.

Should I use cross-validation or a hold-out set?

K-fold cross-validation makes better use of a small control set and gives a more stable estimate, which matters because control points are expensive. A single hold-out set is simpler and is fine when you have plenty of control. Either way the point is the same: the degree must be chosen on points the fit has not seen.