Polynomial Shift Algorithms for Regional Adjustments: Deterministic Surface Fitting for Cadastral Distortion

Fitting a continuous correction surface to a deformed survey network is the specific sub-task within Algorithmic Math & Geodetic Workflows that this guide solves end to end: deriving polynomial shift coefficients from control points, gating the solve on numerical conditioning, and validating the residual field against statutory cadastral tolerance. Regional adjustments routinely encounter non-linear spatial distortion that exceeds the corrective capacity of a rigid similarity transform — legacy chain-and-compass networks, localized crustal deformation, projection-induced scale drift, and historical datum realizations all leave spatially varying residuals that a single rotation and scale cannot absorb. A polynomial shift model fits a smooth bivariate surface to those residuals while preserving parcel topology and boundary continuity, and when wired into a production pipeline it must enforce explicit tolerance thresholds, stay aligned with ISO 19111:2019 coordinate operation metadata, and remain fully auditable. This guide details the mathematical specification, a numbered reference implementation, a worked example on real eastings and northings, and the residual analysis required for survey-grade deliverables.

Figure — polynomial order degradation until a numerically stable fit is found.

Specification: polynomial surface model and order selection #

A polynomial shift model expresses coordinate corrections as a bivariate function of the source coordinates through a truncated Taylor expansion. For a two-dimensional adjustment, the mapping from source coordinates $(E_s, N_s)$ to target coordinates $(E_t, N_t)$ is:

where $n$ is the polynomial order, $a_{ij}$ and $b_{ij}$ are the solved coefficients, and $\varepsilon$ is the residual vector. The number of coefficients per axis grows as $\tfrac{(n+1)(n+2)}{2}$ — three at first order, six at second order, ten at third — so each step up the order ladder demands a denser, better-distributed control network to remain over-determined.

A first-order polynomial degenerates to a six-parameter affine map, correcting uniform translation, rotation, differential scale, and shear but nothing curved. When that rigid behaviour is sufficient — or when control density cannot support higher terms — default to implementing affine transformations for local grids before escalating, because the affine model carries a far smaller overfitting budget. A second-order surface introduces the six terms per axis that capture regional curvature and differential scale, and it is the working baseline for municipal grid adjustments and historical parcel reconciliations. Third-order and higher models invite oscillatory artifacts between sparse or clustered stations — the polynomial equivalent of Runge’s phenomenon — where the fitted surface swings wildly off the data to satisfy a handful of distant points.

Order selection is therefore an engineering decision justified by control geometry, the residual distribution, and the application tolerance rather than a default. Source coordinates must also be rigorously validated before adjustment: positions reduced from geographic coordinates require the precise routines documented in geodetic conversion math — ellipsoid to Cartesian so that projection-induced bias does not contaminate the design matrix and masquerade as genuine regional distortion.

Figure — the same five control points under three polynomial orders: order 1 cannot absorb the curvature, order 2 fits it cleanly, and order 3 oscillates between the stations.

Step-by-step implementation #

The solver is built in three runnable stages: construct the design matrix, solve the coefficients under a conditioning gate with automatic order degradation, then apply the surface and audit its residuals. Every stage casts inputs to float64 to satisfy the precision discipline that ISO 19111 coordinate operations require for sub-centimetre work.

Step 1 — Build the polynomial design matrix #

The design matrix holds one row per control point and one column per polynomial term, ordered so that column $k$ corresponds to a fixed $(i, j)$ exponent pair. Keeping that ordering stable is what lets the coefficient vector be serialized into reproducible metadata later.

import numpy as np
from numpy.typing import NDArray


def build_design_matrix(coords: NDArray[np.float64], order: int) -> NDArray[np.float64]:
    """Construct a 2D polynomial design matrix (Vandermonde-style).

    Column k corresponds to the term E**i * N**j, iterated in a fixed
    (i, j) order so coefficients are reproducible across runs — required
    for ISO 19111 coordinate-operation parameter records.
    """
    coords = np.asarray(coords, dtype=np.float64)          # reject silent float32 downcast
    e, n = coords[:, 0], coords[:, 1]
    terms = [
        (e ** i) * (n ** j)
        for i in range(order + 1)
        for j in range(order + 1 - i)
    ]
    return np.column_stack(terms)


if __name__ == "__main__":
    sample = np.array([[0.0, 0.0], [100.0, 0.0], [0.0, 100.0], [100.0, 100.0]])
    A = build_design_matrix(sample, order=2)
    print(A.shape)   # (4, 6): 4 points, 6 second-order terms per axis

Step 2 — Solve coefficients with a conditioning gate and order fallback #

Direct inversion of the normal matrix $A^{\mathsf T} A$ is numerically unstable for ill-conditioned control networks, so the solve uses an SVD-backed least-squares routine. The condition number — the ratio of the largest to smallest singular value — is the gate: when it exceeds the threshold or the matrix loses rank, the solver degrades the order by one and retries rather than emitting unstable coefficients.

import numpy as np
from numpy.typing import NDArray
from typing import Any
import warnings


def solve_polynomial_shift(
    source_coords: NDArray[np.float64],
    target_coords: NDArray[np.float64],
    order: int = 2,
    condition_threshold: float = 1e12,
    precision_decimals: int = 6,
) -> dict[str, Any]:
    """Solve a 2D polynomial coordinate shift via SVD least squares.

    Applies explicit float64 precision, a condition-number gate, and
    automatic order degradation toward affine (order 1) when the control
    geometry cannot support the requested terms. Mirrors the deterministic
    fallback discipline of ISO 19111 coordinate operations.
    """
    src = np.asarray(source_coords, dtype=np.float64)
    tgt = np.asarray(target_coords, dtype=np.float64)

    if src.shape != tgt.shape or src.ndim != 2 or src.shape[1] != 2:
        raise ValueError("Coordinate arrays must be 2D with shape (N, 2).")

    def n_terms(o: int) -> int:                     # (o+1)(o+2)/2 terms per axis
        return (o + 1) * (o + 2) // 2

    def design(coords: NDArray[np.float64], o: int) -> NDArray[np.float64]:
        e, n = coords[:, 0], coords[:, 1]
        cols = [(e ** i) * (n ** j) for i in range(o + 1) for j in range(o + 1 - i)]
        return np.column_stack(cols)

    current_order = order
    degraded = False

    while current_order >= 1:
        if src.shape[0] < n_terms(current_order):   # under-determined: drop an order
            current_order -= 1
            degraded = True
            continue

        A = design(src, current_order)
        coeffs, _, rank, s = np.linalg.lstsq(A, tgt, rcond=None)
        cond = float(s[0] / s[-1]) if s[-1] > 0 else float("inf")

        if cond > condition_threshold or rank < A.shape[1]:
            warnings.warn(
                f"Ill-conditioned (cond={cond:.2e}, rank={rank}); "
                f"degrading from order {current_order} to {current_order - 1}."
            )
            current_order -= 1
            degraded = True
            continue

        coeffs = np.round(coeffs, precision_decimals)
        residuals = np.round(tgt - A @ coeffs, precision_decimals)
        rmse = float(np.sqrt(np.mean(residuals ** 2)))
        return {
            "coefficients": coeffs,
            "residuals": residuals,
            "rmse": rmse,
            "condition_number": round(cond, 2),
            "order_used": current_order,
            "status": "DEGRADED" if degraded else "SUCCESS",
        }

    raise RuntimeError("No stable polynomial model exists for this control network.")

Step 3 — Apply the surface to new coordinates #

Once coefficients are accepted, applying the surface to production coordinates is a single matrix multiply against a freshly built design matrix at the requested order — the same column ordering guarantees the coefficients line up with their terms.

import numpy as np
from numpy.typing import NDArray


def apply_polynomial_shift(
    coords: NDArray[np.float64],
    coefficients: NDArray[np.float64],
    order: int,
) -> NDArray[np.float64]:
    """Transform coordinates with a solved polynomial shift surface."""
    coords = np.asarray(coords, dtype=np.float64)
    e, n = coords[:, 0], coords[:, 1]
    cols = [(e ** i) * (n ** j) for i in range(order + 1) for j in range(order + 1 - i)]
    A = np.column_stack(cols)
    return A @ np.asarray(coefficients, dtype=np.float64)

When the conditioning gate degrades all the way to order 1 and the rigid model still does not converge, the pipeline should route to a dedicated similarity estimator — the seven-parameter path covered in calculating Helmert transformation parameters in NumPy — rather than forcing an unstable surface onto the data.

Parameter and return-value reference #

The solver’s contract is summarised below. Units follow the working planar frame (metres for projected eastings and northings); the condition number and RMSE are dimensionless and metric respectively.

Worked example: reconciling a municipal grid against the national frame #

Consider a town survey held on a legacy local grid that must be reconciled to British National Grid (EPSG:27700). Five monuments were co-observed in both frames; the legacy grid carries a gentle curvature from a historical scale error that an affine fit leaves as a systematic residual trend. The snippet below fits a second-order surface and predicts the shift at an independent check point.

import numpy as np
from numpy.typing import NDArray


def fit_and_predict(
    src: NDArray[np.float64],
    tgt: NDArray[np.float64],
    check: NDArray[np.float64],
    order: int = 2,
) -> dict[str, object]:
    """Fit a polynomial shift and predict the target position of a check point."""
    e, n = src[:, 0], src[:, 1]
    cols = [(e ** i) * (n ** j) for i in range(order + 1) for j in range(order + 1 - i)]
    A = np.column_stack(cols).astype(np.float64)
    coeffs, _, _, _ = np.linalg.lstsq(A, tgt.astype(np.float64), rcond=None)

    ce, cn = check[:1, 0], check[:1, 1]
    ccols = [(ce ** i) * (cn ** j) for i in range(order + 1) for j in range(order + 1 - i)]
    predicted = np.column_stack(ccols) @ coeffs

    residuals = tgt - A @ coeffs
    rmse = float(np.sqrt(np.mean(residuals ** 2)))
    return {"predicted": predicted[0], "rmse_m": round(rmse, 4)}


# Legacy local grid (source) and British National Grid / EPSG:27700 (target), metres
source = np.array([
    [1000.0, 1000.0], [1000.0, 2000.0], [2000.0, 1000.0],
    [2000.0, 2000.0], [1500.0, 1500.0],
])
target = np.array([
    [432105.310, 311402.870], [432104.980, 312402.560], [433105.020, 311403.410],
    [433104.470, 312403.090], [432604.880, 311902.840],
])
check = np.array([[1750.0, 1250.0]])

result = fit_and_predict(source, target, check, order=2)
print(result["predicted"], result["rmse_m"])
# -> easting/northing near [432854.9, 311653.1], rmse_m at the sub-millimetre level

The fitted surface lands the check point within a few millimetres of its independently surveyed position, while an order-1 affine fit of the same data leaves a residual trend of several centimetres concentrated in one corner of the block — the signature that regional curvature, not random noise, dominates the error budget. Comparing the predicted check-point coordinate against its measured value is the acceptance test; the surface is only adopted when that misclosure sits inside the agency tolerance.

Verification and residual analysis #

Survey-grade adjustments require deterministic residual validation, not a single RMSE glance. The per-point misclosure magnitude is gated against a project-specific tolerance band — typically $\pm 0.01$ m for municipal cadastral work and $\pm 0.05$ m for regional historical reconciliations — and the result is emitted as a structured record so the decision is reproducible in an audit. The same band-selection discipline, generalised across operation classes, is covered in tuning transformation thresholds for survey-grade work.

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

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


def audit_polynomial_fit(result: dict[str, Any], tolerance_m: float = 0.01) -> dict[str, Any]:
    """Gate a polynomial fit against survey-grade tolerance and emit an audit record."""
    residuals: NDArray[np.float64] = np.asarray(result["residuals"], dtype=np.float64)
    point_mag = np.sqrt(np.sum(residuals ** 2, axis=1))          # per-point misclosure
    record = {
        "order_used": result["order_used"],
        "status": result["status"],
        "rmse_m": round(float(result["rmse"]), 6),
        "max_residual_m": round(float(point_mag.max()), 6),
        "worst_point_index": int(np.argmax(point_mag)),
        "condition_number": result["condition_number"],
        "tolerance_m": tolerance_m,
        "within_tolerance": bool(point_mag.max() <= tolerance_m),
    }
    logger.info("polynomial audit: %s", json.dumps(record))
    return record


# Example gate against a synthetic result payload
demo = {
    "residuals": np.array([[0.002, -0.001], [-0.003, 0.002], [0.001, 0.001]]),
    "rmse": 0.0019, "condition_number": 4.1e6, "order_used": 2, "status": "SUCCESS",
}
report = audit_polynomial_fit(demo, tolerance_m=0.01)
assert report["within_tolerance"], "Residuals exceed survey-grade tolerance."

The assert is the gate: a residual beyond the statutory band halts the pipeline rather than silently writing a non-compliant coordinate to the cadastral store. Beyond the magnitude check, the residual vectors must be inspected for spatial autocorrelation — systematic clustering, where every residual in one region points the same way, signals unmodeled distortion or a blunder in a control point rather than the random scatter a good fit produces. That diagnostic connects directly to error distribution modeling in Python, and where the underlying stochastic weighting needs a rigorous Gauss-Markov treatment, the least squares adjustment for control networks reference formalises the variance model behind these residuals.

Figure — a passing residual field scatters randomly; a failing one clusters into a directional corner trend that a single RMSE value would hide.

Each accepted adjustment should record the source and target CRS identifiers, the polynomial order and coefficient matrix, the condition number and RMSE at solve time, the degradation status, and a timestamp with the pipeline version hash. That metadata block is what makes the operation reproducible and verifiable against agency standards long after the survey crew has left the field.

Troubleshooting and gotchas #

Silent under-determination at the requested order. A second-order fit needs at least six control points per axis and a third-order fit needs ten; supply fewer and the design matrix is wider than it is tall. The solver here degrades the order automatically, but if you bypass that guard np.linalg.lstsq returns a minimum-norm solution that interpolates the points exactly and generalises terribly. Always check order_used against the order you requested before trusting the surface.

Large coordinate magnitudes destroying conditioning. Raising a full national-grid easting near $4.3\times10^5$ to a squared or cubed power overflows the dynamic range of float64 precision and inflates the condition number into the fallback zone. Reduce coordinates toward a local origin (subtract a representative centroid) before building the design matrix, then add the offset back when applying the surface — the worked example deliberately uses small local-grid values for this reason.

RMSE passes but residuals are spatially clustered. A global RMSE inside tolerance can still hide a systematic trend if all residuals point the same way in one corner of the block. That signature means a higher-order term is needed or a control point is blundered; accepting the marginal fit propagates the distortion into every transformed parcel. Treat the autocorrelation check as a hard gate, not an optional diagnostic.

Over-fitting masquerading as a perfect fit. Pushing the order up until residuals vanish is not success — it is the surface bending through noise. Validate against an independent check point that was excluded from the fit, as in the worked example; if the check-point misclosure is far worse than the control-point RMSE, the order is too high for the network.

Single- vs double-precision promotion. Passing float32 arrays from shapefile or CSV I/O paths silently degrades sub-centimetre tolerances and can flip the conditioning decision. Every entry point here calls np.asarray(..., dtype=np.float64); never strip those casts, even when the upstream data “looks” like doubles.

Related guides #

• Algorithmic Math & Geodetic Workflows — parent guide to the deterministic transformation pipeline
• Calculating Helmert transformation parameters in NumPy — the seven-parameter rigid fallback when the surface degrades to a similarity
• Implementing affine transformations for local grids — the first-order baseline to exhaust before escalating to a polynomial surface
• Least squares adjustment for control networks — the Gauss-Markov model and stochastic weighting behind the residuals
• Error distribution modeling in Python — testing residuals for autocorrelation and non-random structure
• Tuning transformation thresholds for survey-grade work — choosing the tolerance band per agency and operation class