Least Squares Adjustment of Control Points: A Deterministic Python Affine Fit

Fitting the six parameters of a plane affine transformation for local grids from redundant control points is a weighted least squares problem that must resolve to a survey-grade tolerance — typically a residual band of 0.005 m for cadastral control — while staying deterministic and auditable under ISO 19111:2019 coordinate operation metadata. This page is a focused implementation within the broader Algorithmic Math & Geodetic Workflows reference: given matched source and target monuments, it estimates the affine coefficients, propagates a unit-weight variance, and gates the result against a statutory threshold before any coordinate is committed.

The affine model and its weighted normal equations #

A two-dimensional affine transformation maps each source coordinate $(x, y)$ to a target coordinate $(x', y')$ through six independent coefficients — four describing combined scale, rotation, and shear, and two describing translation:

Each control point contributes two observation equations, so $n$ points yield $2n$ rows in a design matrix $\mathbf{A}$ of width six. Stacking the unknowns as $\mathbf{p} = [a_{11}, a_{12}, t_x, a_{21}, a_{22}, t_y]^{\mathsf{T}}$ and the observed target values as $\boldsymbol{\ell}$, the residual model is $\mathbf{v} = \mathbf{A}\mathbf{p} - \boldsymbol{\ell}$. With a diagonal weight matrix $\mathbf{P} = \mathbf{\Sigma}_\ell^{-1}$ built from per-point precision, the weighted least squares solution comes from the normal equations

The quality of the fit is summarised by the a-posteriori unit-weight standard deviation $\hat{\sigma}_0 = \sqrt{\mathbf{v}^{\mathsf{T}}\mathbf{P}\mathbf{v}\,/\,(2n-6)}$, where $2n-6$ is the redundancy. Three non-collinear points give a zero-redundancy exact fit; the adjustment becomes meaningful only with four or more. Because this estimator is the linear ($\text{order}=1$) case of the same machinery used for polynomial shift algorithms for regional adjustments, its residuals and covariance feed directly into downstream error distribution modeling.

Complete runnable implementation #

The function below is self-contained: it accepts matched control points and optional per-point standard deviations, builds the design and weight matrices explicitly in float64, monitors the condition number before inverting, and routes to an SVD-based least squares solve when the normal matrix is ill-conditioned. Rounding is deferred — the returned coefficients keep full IEEE 754 precision so the caller controls deliverable formatting.

import numpy as np
from numpy.typing import NDArray
from numpy.linalg import LinAlgError, cond, lstsq, solve

CONDITION_THRESHOLD: float = 1e10   # collinear / degenerate geometry above this


def adjust_affine_control_points(
    source: NDArray[np.float64],
    target: NDArray[np.float64],
    point_sigmas: NDArray[np.float64] | None = None,
    tolerance_m: float = 0.005,        # cadastral residual band, metres
) -> dict[str, object]:
    """Estimate a 2D affine transform from control points by weighted least squares.

    Solves v = A p - l with weights P = Sigma_l^-1 for the six affine
    coefficients [a11, a12, tx, a21, a22, ty], satisfying the linear
    coordinate-operation model of ISO 19111 (clause 5.4, parametric transform).

    Args:
        source: (n, 2) array of source-frame coordinates, metres.
        target: (n, 2) array of matching target-frame coordinates, metres.
        point_sigmas: optional (n,) positional standard deviations, metres;
            defaults to unit weights (ordinary least squares).
        tolerance_m: maximum permitted absolute coordinate residual, metres.

    Returns:
        Dict with affine matrix, translation, residuals, sigma0, RMSE,
        condition number, the solver path taken, and a pass/fail flag.
    """
    src = np.asarray(source, dtype=np.float64)   # explicit double precision
    tgt = np.asarray(target, dtype=np.float64)
    if src.shape != tgt.shape or src.ndim != 2 or src.shape[1] != 2:
        raise ValueError("source and target must be congruent (n, 2) arrays.")
    n = src.shape[0]
    if n < 3:
        raise ValueError("At least three non-collinear control points required.")

    # 2n x 6 design matrix: two observation rows (x', y') per control point.
    A = np.zeros((2 * n, 6), dtype=np.float64)
    A[0::2, 0:3] = np.column_stack([src[:, 0], src[:, 1], np.ones(n)])  # x' row
    A[1::2, 3:6] = np.column_stack([src[:, 0], src[:, 1], np.ones(n)])  # y' row
    L = tgt.reshape(-1).astype(np.float64)        # interleaved [x'_0, y'_0, ...]

    # Diagonal weights P = 1 / sigma^2, one value shared by each point's two rows.
    if point_sigmas is None:
        w = np.ones(2 * n, dtype=np.float64)
    else:
        sig = np.asarray(point_sigmas, dtype=np.float64).reshape(-1)
        if sig.shape[0] != n or np.any(sig <= 0):
            raise ValueError("point_sigmas must hold n positive values.")
        w = np.repeat(1.0 / sig**2, 2)

    # Weighted normal equations N = A^T P A, U = A^T P L (P diagonal -> scale rows).
    Aw = A * w[:, None]
    N = A.T @ Aw
    U = A.T @ (w * L)

    c_num = float(cond(N))
    solver = "cholesky_normal"
    try:
        if c_num > CONDITION_THRESHOLD:
            raise LinAlgError(f"cond(N) = {c_num:.2e} exceeds {CONDITION_THRESHOLD:.0e}.")
        p = solve(N, U)                           # primary path
    except LinAlgError:
        solver = "svd_lstsq_fallback"             # degenerate geometry
        sqrt_w = np.sqrt(w)[:, None]
        p = lstsq(A * sqrt_w, L * sqrt_w[:, 0], rcond=None)[0]

    # Residuals in the target frame and survey-grade quality measures.
    v = (A @ p) - L                               # 2n residual vector, metres
    dof = 2 * n - 6
    sigma0 = float(np.sqrt((v * w) @ v / dof)) if dof > 0 else float("nan")
    point_residual = np.sqrt(v[0::2] ** 2 + v[1::2] ** 2)   # per-point magnitude
    rmse = float(np.sqrt(np.mean(point_residual**2)))
    max_residual = float(np.max(np.abs(v)))

    return {
        "affine_matrix": p[[0, 1, 3, 4]].reshape(2, 2),   # [[a11, a12], [a21, a22]]
        "translation": p[[2, 5]],                          # [tx, ty], metres
        "residuals": v.reshape(n, 2),
        "sigma0": sigma0,
        "rmse_m": rmse,
        "max_residual_m": max_residual,
        "condition_number": c_num,
        "degrees_of_freedom": dof,
        "solver_path": solver,
        "within_tolerance": bool(max_residual <= tolerance_m),
    }

Parameter and return-value reference #

Units follow the working plane: metres for coordinates and residuals, dimensionless for the variance factor and condition number.

Minimal worked example #

A legacy local grid is tied to four monuments observed on the British National Grid (EPSG:27700) plane, each with a positional precision near 0.02 m. The fit recovers an offset of roughly $(529100, 181200)$ m with near-unit scale, and the redundant fourth point keeps the residuals honest.

import numpy as np

src = np.array([[0.0, 0.0], [100.0, 0.0], [0.0, 100.0], [100.0, 100.0]], dtype=np.float64)
tgt = np.array(
    [[529100.00, 181200.00],
     [529200.01, 181199.98],
     [529100.02, 181300.01],
     [529199.99, 181299.97]],
    dtype=np.float64,
)
sigmas = np.full(src.shape[0], 0.02, dtype=np.float64)   # 20 mm per monument

result = adjust_affine_control_points(src, tgt, point_sigmas=sigmas, tolerance_m=0.05)
print(f"translation (m) : {np.round(result['translation'], 3)}")
print(f"sigma0          : {result['sigma0']:.3f}")
print(f"max residual (m): {result['max_residual_m']:.4f}")
print(f"solver path     : {result['solver_path']}")
# Expected (approx): translation ~ [529100.0, 181200.0], max residual ~ 0.01 m

Validation check #

The pass/fail decision must be explicit so a non-compliant transformation cannot silently write coordinates to the cadastral store. The returned within_tolerance flag is the gate; assert on it before propagating any parameter, exactly as the threshold discipline in tuning transformation thresholds for survey-grade work prescribes.

assert result["within_tolerance"], (
    f"Affine fit rejected: max residual {result['max_residual_m']:.4f} m "
    f"exceeds survey-grade tolerance."
)

A second, independent acceptance test withholds one monument from the fit and re-predicts it; its transformed position should agree with the held-back coordinate inside the same tolerance, which is the cross-check formalised in validating datum alignment with control points.

Common mistakes #

Single- vs double-precision promotion. Control coordinates arriving as float32 from a shapefile or CSV reader silently degrade the normal matrix; on a national grid the easting magnitude (~5×10⁵) leaves only millimetre resolution in single precision. Every entry point here calls np.asarray(..., dtype=np.float64) — never strip those casts to “save memory”.

Treating the exact fit as adjusted. With exactly three control points the redundancy $2n-6$ is zero, sigma0 is undefined, and the residuals are identically zero by construction. That is interpolation, not adjustment: it reports no quality information at all. Add at least one redundant monument before trusting within_tolerance.

Forgetting the per-row weight repeat. Each control point owns two observation rows (its x′ and y′ equations), so a per-point variance must be repeated across both rows with np.repeat(..., 2). Applying one weight per point instead of per row mis-scales the normal equations and pushes sigma0 away from unity even when the geometry is sound.

Related guides #

• Implementing affine transformations for local grids — parent guide to affine parameter derivation and similarity fallback
• Algorithmic Math & Geodetic Workflows — the deterministic transformation pipeline this fit belongs to
• Least squares adjustment for control networks — the full Gauss-Markov model and chi-square test behind sigma0
• Error distribution modeling in Python — propagating these residuals into confidence ellipses
• Validating datum alignment with control points — the held-back-monument acceptance test