Polynomial Shift Algorithms for Regional Adjustments

Regional coordinate adjustments in cadastral and high-precision surveying routinely encounter non-linear spatial distortions that exceed the corrective capacity of rigid similarity transforms. Legacy survey networks, localized crustal deformation, projection-induced scale drift, and historical datum realizations require deterministic polynomial shift models capable of absorbing spatially varying residuals while preserving parcel topology and boundary continuity. When integrated into production coordinate transformation pipelines, these algorithms must enforce explicit tolerance thresholds, maintain strict alignment with ISO 19111 coordinate reference system metadata, and remain fully auditable against EPSG registry conventions. The mathematical and operational foundations for these workflows reside within the broader discipline of Algorithmic Math & Geodetic Workflows, where least-squares adjustment, condition equation formulation, and error propagation are treated as engineered processes rather than heuristic approximations.

flowchart TD
    A["Source ↔ target control points"] --> B["Build design matrix<br/>order n (default 2)"]
    B --> C["SVD least-squares solve"]
    C --> D{"cond ≤ threshold<br/>and full rank?"}
    D -->|"yes"| E(["Coefficients + RMSE"])
    D -->|"no"| F["Reduce order: n → n−1"]
    F --> G{"n ≥ 1?"}
    G -->|"yes"| B
    G -->|"no"| H(["Fail — no stable model"])

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

Mathematical Formulation and Order Selection #

Polynomial shift models express coordinate corrections as a function of source coordinates through a truncated Taylor series expansion. For a two-dimensional cadastral adjustment, the transformation from source coordinates $(E_s, N_s)$ to target coordinates $(E_t, N_t)$ is defined as:

where $n$ denotes the polynomial order, $a_{ij}$ and $b_{ij}$ are the solved coefficients, and $\varepsilon$ represents the residual vector. First-order polynomials degenerate to affine transformations, which correct uniform scale, rotation, and translation but cannot capture regional curvature or differential scale. Second-order models introduce six coefficients per axis and serve as the industry baseline for municipal grid adjustments, historical parcel reconciliations, and localized datum shifts. Third-order and higher models introduce significant overfitting risk when control networks are sparse or geometrically clustered, frequently producing oscillatory artifacts between control stations.

Order selection must be justified by control geometry, residual distribution, and application tolerance. When rigid-body behavior is sufficient or when control density is insufficient to support higher-order terms, practitioners should default to Implementing Affine Transformations for Local Grids before escalating to polynomial formulations. Input coordinates must be rigorously validated prior to adjustment; source data derived from geographic coordinates require precise conversion routines as documented in Geodetic Conversion Math: Ellipsoid to Cartesian to prevent projection-induced bias from contaminating the design matrix.

Deterministic Implementation Architecture #

Production-grade polynomial solvers must prioritize numerical stability over algorithmic complexity. The adjustment is solved via weighted least squares using the design matrix $A$, observation vector $L$, and parameter vector $X$:

The normal equations $A^T P A X = A^T P L$ yield the coefficient vector, but direct inversion of $A^T P A$ is numerically unstable for ill-conditioned control networks. Production pipelines should employ QR decomposition or SVD-based least-squares solvers, enforce explicit floating-point precision boundaries, and implement automatic fallback routing when matrix rank deficiency or excessive condition numbers are detected. When control geometry fails to support polynomial terms, the pipeline must gracefully degrade to rigid parameter estimation, following the methodology outlined in Calculating Helmert transformation parameters in NumPy.

import numpy as np
from typing import Tuple, Dict, Union, Optional
import warnings

def solve_polynomial_shift(
    source_coords: np.ndarray,
    target_coords: np.ndarray,
    order: int = 2,
    condition_threshold: float = 1e12,
    precision_decimals: int = 6
) -> Dict[str, Union[np.ndarray, float, str]]:
    """
    Solves a 2D polynomial coordinate shift using least-squares with explicit
    precision control and automatic fallback routing to affine (order=1).

    Parameters
    ----------
    source_coords : np.ndarray
        Shape (N, 2) array of [E, N] source coordinates.
    target_coords : np.ndarray
        Shape (N, 2) array of [E, N] target coordinates.
    order : int
        Polynomial order (1 or 2). Higher orders require dense, well-distributed control.
    condition_threshold : float
        Maximum acceptable condition number before triggering fallback.
    precision_decimals : int
        Explicit decimal places for output coefficients and residuals.

    Returns
    -------
    dict
        Contains 'coefficients', 'residuals', 'condition_number', 'status', and 'order_used'.
    """
    # Enforce explicit float64 precision to prevent silent downcasting
    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).")

    if src.shape[0] < 3:
        raise ValueError("Minimum 3 control points required for polynomial adjustment.")

    def _build_design_matrix(coords: np.ndarray, poly_order: int) -> np.ndarray:
        """Constructs Vandermonde-style design matrix for 2D polynomial."""
        terms = []
        for i in range(poly_order + 1):
            for j in range(poly_order + 1 - i):
                terms.append((coords[:, 0] ** i) * (coords[:, 1] ** j))
        return np.column_stack(terms)

    current_order = order
    fallback_triggered = False

    while current_order >= 1:
        try:
            A = _build_design_matrix(src, current_order)
            # Solve via SVD-backed least squares for numerical stability
            coeffs, residuals, rank, s = np.linalg.lstsq(A, tgt, rcond=None)

            # Condition number estimation from singular values
            cond_num = float(s[0] / s[-1]) if s[-1] > 0 else float('inf')

            if cond_num > condition_threshold or rank < A.shape[1]:
                warnings.warn(
                    f"Matrix ill-conditioned (cond={cond_num:.2e}, rank={rank}). "
                    f"Falling back to order {current_order - 1}."
                )
                current_order -= 1
                fallback_triggered = True
                continue

            # Explicit precision truncation
            coeffs = np.round(coeffs, decimals=precision_decimals)
            residuals = np.round(tgt - A @ coeffs, decimals=precision_decimals)
            rmse = float(np.sqrt(np.mean(residuals**2)))

            return {
                "coefficients": coeffs,
                "residuals": residuals,
                "rmse": rmse,
                "condition_number": round(cond_num, 2),
                "status": "FALLBACK_TO_AFFINE" if fallback_triggered else "SUCCESS",
                "order_used": current_order
            }

        except np.linalg.LinAlgError as e:
            if current_order > 1:
                current_order -= 1
                fallback_triggered = True
                continue
            raise RuntimeError(f"Polynomial solver failed at all orders: {e}") from e

    raise RuntimeError("Failed to converge on a valid transformation model.")

Residual Auditing and Threshold Enforcement #

Survey-grade adjustments require deterministic residual validation. Post-solution, the pipeline must evaluate the root-mean-square error (RMSE) against project-specific tolerance bands, typically $\pm 0.01$ m for municipal cadastral work and $\pm 0.05$ m for regional historical reconciliations. Residual vectors must be analyzed for spatial autocorrelation; systematic clustering indicates unmodeled distortion or erroneous control points. When residuals exceed tolerance, the pipeline should flag the adjustment for manual review rather than applying forced corrections. Threshold enforcement must be configurable via metadata-driven parameters to support multi-jurisdictional deployments.

Pipeline Integration and Metadata Alignment #

Production coordinate transformation pipelines must embed ISO 19111 coordinate operation metadata directly into the adjustment output. Each solved polynomial shift requires explicit recording of:

• Source and target CRS identifiers (EPSG codes or URNs)
• Polynomial order and coefficient matrix
• Condition number and RMSE at solve time
• Fallback routing status and order degradation path
• Timestamp and pipeline version hash

This metadata ensures full auditability, supports reproducible geospatial workflows, and enables automated compliance validation against agency standards. By enforcing explicit precision boundaries, deterministic fallback routing, and strict residual thresholds, polynomial shift algorithms transition from theoretical approximations to auditable, pipeline-ready engineering components.