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,

One bootstrap replicate, and what must be inside the loop The control set is sampled with replacement to the same size. The entire estimator is then refitted on that sample — the whole point is that parameter estimation is inside the loop, because a bootstrap that only resamples the final statistic measures nothing about the fit. The statistic of interest is recorded and the loop repeats B times. The recorded replicates form the sampling distribution from which the percentile interval is read. Control set n monuments Resample with replacement same size n Refit the FULL estimator parameters, not just the statistic Record the statistic one replicate B replicates percentile interval after B repeat

Figure — the resampling loop: refit the whole estimator, not just the statistic.

θ^(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.

Stability of a bootstrap interval against the number of replicates The run-to-run standard deviation of the estimated 95 percent bound, in metres, against the number of bootstrap replicates, computed from the one-over-square-root-of-B law scaled to 0.0142 metres at one thousand replicates. Fifty replicates wobble by about 0.064 metres, which is larger than the tolerance being tested. Beyond a thousand replicates the gain is marginal, which is why a thousand is the usual default. 0.000 0.020 0.040 0.060 50 250 500 1000 2000 bootstrap replicates B m run-to-run sd of the bound (m)

Figure — how unstable the interval is at small B: the estimate wobbles as 1/sqrt(B).

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.

Frequently Asked Questions

How many resamples are enough?

One thousand is the usual working default for a percentile interval, and the reason is visible in the stability of the estimate: the run-to-run variation of the interval bound falls as one over the square root of the replicate count, so fifty replicates give a bound that moves by more than the tolerance you are testing, while going from one thousand to ten thousand buys a factor of three for ten times the compute. If the decision is marginal, more replicates will not rescue it — more control points will.

Should I resample the residuals or the control points?

Resample the control points, together with everything attached to them, and refit the whole estimator on each sample. Resampling residuals alone holds the fitted parameters fixed, which measures the spread of the residuals rather than the uncertainty of the transformation — a much smaller number, and not the one being asked for. The rule is that whatever you want an uncertainty for must be recomputed inside the loop.

What does it mean when the bootstrap interval is wildly asymmetric?

Usually that one or two control points dominate the fit. When a single monument appears in a resample the estimate lands in one place, and when it is omitted it lands somewhere quite different, producing a bimodal distribution that a percentile interval reports as asymmetry. That is real information: the transformation depends on one observation more than it should, and the fix is more control, not a different interval method.