CoordTransform Cadastral & High-Precision Python

Bootstrap Resampling for Transformation Uncertainty

When the closed-form covariance of a transformation is untrustworthy — small control sets, non-Gaussian residuals, or a fit whose linearization is shaky — resampling the control-point residuals gives an empirical uncertainty that makes no distributional promises it cannot keep. This technique extends error distribution modeling in Python under the wider Algorithmic Math & Geodetic Workflows reference, replacing an assumed error law with the one your own residuals actually exhibit.

Residual bootstrap for transformation parameters The original residuals feed a resampling block that draws them with replacement using a seeded generator. Each resample is added back to the fitted control values and the transformation is re-fitted, producing one parameter estimate. After B replicates the accumulated estimates form an empirical distribution from which percentile confidence intervals are read. Residuals from fit Resample w/ replace seeded RNG · re-fit → θ*(b) Accumulate B replicates Percentile CI repeat B times

Figure — each replicate re-fits from resampled residuals; the spread is the uncertainty.

Why Resample the Residuals

A least-squares transformation returns parameters θ^\hat{\boldsymbol\theta} and, in principle, a covariance σ02(ATPA)−1\sigma_0^2 (\mathbf{A}^{\mathsf T}\mathbf{P}\mathbf{A})^{-1}. That covariance is only as honest as its assumptions: that the residuals are zero-mean, homoscedastic, and roughly normal, and that the model is locally linear at the solution. When a cadastral fit rests on a handful of monuments, or the residuals are visibly skewed, those assumptions fail quietly and the reported error ellipse understates reality. The residual bootstrap sidesteps the assumption entirely. Holding the fitted control values fixed, it draws a resample r∗\mathbf{r}^{*} from the observed residuals with replacement, forms synthetic observations y∗=y^+r∗\mathbf{y}^{*} = \hat{\mathbf{y}} + \mathbf{r}^{*}, and re-fits. Repeating this BB times yields an empirical distribution of the estimator,

θ^(b)∗=f ⁣(y^+r(b)∗),b=1,…,B\hat{\boldsymbol\theta}^{*}_{(b)} = f\!\left(\hat{\mathbf{y}} + \mathbf{r}^{*}_{(b)}\right), \qquad b = 1, \dots, B

whose spread estimates the sampling variability of θ^\hat{\boldsymbol\theta}. A two-sided 100(1−α)%100(1-\alpha)\% confidence interval for any parameter or transformed coordinate is then read straight off the order statistics — the α/2\alpha/2 and 1−α/21-\alpha/2 percentiles of the BB replicates. Because the intervals come from the residuals themselves, they inherit whatever non-normal shape those residuals have, complementing the analytic ellipses drawn when visualizing geodetic error ellipses with Matplotlib.

Complete Runnable Implementation

The routine below bootstraps a four-parameter (similarity) planar transformation, but it accepts any fitter matching its signature. It fixes the RNG seed for reproducibility, resamples residuals rather than points, and reports percentile intervals for every parameter.

from __future__ import annotations

from collections.abc import Callable
from dataclasses import dataclass

import numpy as np


@dataclass(frozen=True)
class BootstrapCI:
    point_estimate: np.ndarray  # theta_hat from the original fit
    lower: np.ndarray           # per-parameter lower percentile bound
    upper: np.ndarray           # per-parameter upper percentile bound
    replicates: np.ndarray      # B-by-k array of bootstrap estimates


def fit_similarity(src: np.ndarray, dst: np.ndarray) -> np.ndarray:
    """Least-squares 4-parameter (a, b, tx, ty) planar similarity fit.

    Maps src -> dst by [X; Y] = [[a, -b], [b, a]] @ [x; y] + [tx; ty].
    Returns theta = (a, b, tx, ty).
    """
    x, y = src[:, 0], src[:, 1]
    n = x.size
    top = np.column_stack([x, -y, np.ones(n), np.zeros(n)])
    bot = np.column_stack([y, x, np.zeros(n), np.ones(n)])
    design = np.vstack([top, bot])
    obs = np.concatenate([dst[:, 0], dst[:, 1]])
    theta, *_ = np.linalg.lstsq(design, obs, rcond=None)
    return theta


def bootstrap_transform_uncertainty(
    src: np.ndarray,
    dst: np.ndarray,
    fitter: Callable[[np.ndarray, np.ndarray], np.ndarray] = fit_similarity,
    n_boot: int = 2000,
    alpha: float = 0.05,
    seed: int = 20260714,
) -> BootstrapCI:
    """Percentile confidence intervals for transformation parameters.

    Resamples the fit residuals with replacement, re-fits n_boot times, and
    returns the alpha/2 and 1-alpha/2 percentile bounds per parameter.

    Args:
        src: n-by-2 source control coordinates (metres).
        dst: n-by-2 target control coordinates (metres).
        fitter: callable(src, dst) -> parameter vector; default similarity.
        n_boot: number of bootstrap replicates; >= 1000 for stable tails.
        alpha: two-sided significance (0.05 -> 95% interval).
        seed: fixed RNG seed for a reproducible, auditable result.

    Returns:
        A BootstrapCI with the point estimate, bounds, and all replicates.
    """
    src = np.asarray(src, dtype=np.float64)
    dst = np.asarray(dst, dtype=np.float64)
    if src.shape != dst.shape or src.shape[1] != 2:
        raise ValueError("src and dst must be matching n-by-2 arrays")
    n = src.shape[0]

    theta_hat = fitter(src, dst)
    fitted = _apply(theta_hat, src)          # hat{y}
    residuals = dst - fitted                 # per-point (dx, dy) residuals

    rng = np.random.default_rng(seed)
    reps = np.empty((n_boot, theta_hat.size), dtype=np.float64)
    for b in range(n_boot):
        pick = rng.integers(0, n, size=n)    # resample residual rows w/ replacement
        synthetic = fitted + residuals[pick]
        reps[b] = fitter(src, synthetic)

    lo = np.percentile(reps, 100.0 * (alpha / 2.0), axis=0)
    hi = np.percentile(reps, 100.0 * (1.0 - alpha / 2.0), axis=0)
    return BootstrapCI(theta_hat, lo, hi, reps)


def _apply(theta: np.ndarray, pts: np.ndarray) -> np.ndarray:
    a, b, tx, ty = theta
    x, y = pts[:, 0], pts[:, 1]
    return np.column_stack([a * x - b * y + tx, b * x + a * y + ty])

Parameter Reference

Name Type Units Valid range Cadastral significance
src, dst np.ndarray metres n≥3, shape (n,2) matched control coordinates in both frames
fitter Callable — returns θ the transformation model being characterised
n_boot int count ≥ 1000 replicate count; too few gives noisy tails
alpha float — 0<x<1 two-sided level; 0.05 → 95% interval
seed int — any fixes the RNG so the interval is reproducible
lower, upper np.ndarray param units — percentile bounds per parameter

Worked Example

A similarity fit over six control points, lightly perturbed, returns a tight interval around the scale/rotation block and the translations. Re-running with the same seed reproduces the bounds exactly.

import numpy as np

src = np.array([[0.0, 0.0], [100.0, 0.0], [100.0, 100.0],
                [0.0, 100.0], [50.0, 20.0], [20.0, 80.0]])
rng = np.random.default_rng(1)
true = _apply(np.array([1.00002, 0.0003, 12.5, -7.25]), src)
dst = true + rng.normal(0, 0.004, true.shape)     # 4 mm control noise
ci = bootstrap_transform_uncertainty(src, dst, n_boot=2000, seed=42)
print(np.round(ci.point_estimate, 5))
print(np.round(ci.lower, 5))
print(np.round(ci.upper, 5))
# [ 1.00005  0.00028 12.49882 -7.24973]
# [ 1.00002  0.00026 12.49676 -7.25205]
# [ 1.00007  0.00031 12.50105 -7.24759]

Validation Check

Assert the point estimate lies inside its own interval and that the interval is reproducible under a fixed seed before the numbers reach an uncertainty statement.

assert np.all(ci.lower <= ci.point_estimate), "point estimate below its lower bound"
assert np.all(ci.point_estimate <= ci.upper), "point estimate above its upper bound"
again = bootstrap_transform_uncertainty(src, dst, n_boot=2000, seed=42)
assert np.allclose(ci.lower, again.lower), "intervals not reproducible — check the seed"

Common Mistakes

Using too few resamples
Percentile bounds are order statistics of the replicate cloud, so the tails are noisy when n_boot is small — a few hundred draws can move a 95% bound by millimetres between runs. Use at least a thousand, and more when you need the 99% tail, then confirm the interval is stable by re-running with a different seed.
Resampling the points instead of the residuals
Resampling whole control points (a case bootstrap) changes the network geometry every replicate and can drop a critical monument entirely, inflating the interval for reasons that have nothing to do with observation error. The residual bootstrap holds the fitted values fixed and resamples only the residuals, so the geometry is preserved and the spread reflects measurement noise.
Not fixing the random seed
An unseeded generator yields different intervals on every run, which is indefensible in an audit trail. Seed the generator explicitly with np.random.default_rng(seed) and record the seed alongside the interval so the exact uncertainty figure can be regenerated.
Ignoring spatial correlation in the residuals
The plain residual bootstrap assumes the per-point residuals are exchangeable. If the residuals are spatially correlated — common when a regional distortion is under-modelled — resampling them independently understates the true uncertainty. Whiten or block-resample the residuals, or model the distortion first, before trusting the interval.