Low GCP Residuals, but the Georeferenced Map Is Still Wrong
On this page
Problem statement
The georeferencer's table of control points says every residual is a metre or two, the RMS is excellent — and the map is visibly wrong: roads miss their modern successors at the edge of the sheet, or the whole map sits a constant distance from where it should. Residuals measure how well the transformation passes through the points it was fitted to. They say nothing about the places in between, and nothing about errors that every point shares.
This page shows the three ways a low residual hides a bad map, on the Ordnance Survey six-inch sheet of Oakham (1884–1891) with seven control points matched to OpenStreetMap, and the number to look at instead.
Quick answer
Fit with fewer points or a more flexible transformation and the residual falls while the real error rises. Compare the residual with the error at points the fit did not see:
print(f"all 7 points: residual RMS {rms(errors(fit(GCPS), GCPS)):4.1f} m, leave-one-out RMS {rms(loo(GCPS)):4.1f} m")
for four in LOW_FOUR:
rest = [g for g in GCPS if g not in four]
print(f"4 of 7 points: residual RMS {rms(errors(fit(four), four)):4.1f} m, "
f"worst error at the 3 points left out {max(errors(fit(four), rest)):4.1f} m")all 7 points: residual RMS 10.2 m, leave-one-out RMS 17.5 m
4 of 7 points: residual RMS 1.9 m, worst error at the 3 points left out 23.0 m
4 of 7 points: residual RMS 1.9 m, worst error at the 3 points left out 39.7 m
Across every subset of the seven points, the fits with the smallest residuals are the ones that miss the points they left out by the most:
What a residual is, and how fewer points mislead
1. Know what a residual is
A residual is the distance between where the fitted transformation puts a control point and where you said it was. It is computed on the same points that determined the fit, so it rewards a fit for agreeing with its own inputs:
A = fit(GCPS)
for (name, *_), r in zip(GCPS, errors(A, GCPS)):
print(f"{name:24} residual {r:5.1f} m")oakham_church residual 8.8 m
egleton_church residual 5.1 m
brooke_rd_crossing residual 11.2 m
egleton_level_crossing residual 14.6 m
stamford_rd_guide_post residual 12.5 m
braunston_rd_bridleway residual 9.6 m
brooke_guide_post residual 5.7 m
2. Fewer points, lower residuals, worse map
With four points an affine fit has six coefficients and eight constraints; there is little left over to disagree. Two of the 35 four-point subsets fit each other to under 2 m:
for four in LOW_FOUR:
print([g[0] for g in four])
for g in GCPS:
if g not in four:
print(f" not used: {g[0]:24} {errors(fit(four), [g])[0]:5.1f} m")['egleton_church', 'stamford_rd_guide_post', 'braunston_rd_bridleway', 'brooke_guide_post']
not used: oakham_church 23.0 m
not used: brooke_rd_crossing 21.2 m
not used: egleton_level_crossing 20.3 m
['brooke_rd_crossing', 'egleton_level_crossing', 'braunston_rd_bridleway', 'brooke_guide_post']
not used: oakham_church 16.1 m
not used: egleton_church 28.6 m
not used: stamford_rd_guide_post 39.7 m
Both have a residual a fifth of the seven-point fit's, and one misses a point it did not use by 39.7 m, nearly twice the seven-point fit's worst leave-one-out error. The residual rewarded both fits for leaving out the points that would have argued with them.
3. More freedom, lower residuals, worse map
The same happens when the transformation, not the point count, changes. A second-order polynomial through seven points and a thin-plate spline both have lower residuals than the affine fit:
for label, opts in (("affine", ["METHOD=GCP_POLYNOMIAL", "MAX_GCP_ORDER=1"]),
("polynomial 2", ["METHOD=GCP_POLYNOMIAL", "MAX_GCP_ORDER=2"]),
("thin-plate spline", ["METHOD=GCP_TPS"])):
res = gdal_errors(GCPS, GCPS, opts)
lo = [gdal_errors(GCPS[:i] + GCPS[i + 1:], [g], opts)[0] for i, g in enumerate(GCPS)]
print(f"{label:18} residual RMS {rms(res):6.1f} m leave-one-out RMS {rms(lo):7.1f} m")affine residual RMS 10.2 m leave-one-out RMS 17.5 m
polynomial 2 residual RMS 1.9 m leave-one-out RMS 571.1 m
thin-plate spline residual RMS 0.0 m leave-one-out RMS 11.6 m
The spline's residual is zero by construction and its leave-one-out error is the best of the three. The polynomial's residual is under 2 m and its leave-one-out error is over half a kilometre. Residuals ranked them in one order; prediction ranks them in another, and only prediction is about the map.
4. An error every point shares is invisible
If all the control points carry the same error — the same datum mix-up, the same misreading of a grid — the transformation absorbs it and the residuals do not change. Take the coordinates as WGS 84 longitude and latitude but label them OSGB 1936, a common slip with British data:
right = warp_centre([gdal.GCP(lo, la, 0, px, py) for _, px, py, lo, la in GCPS_LL], "EPSG:4326")
wrong = warp_centre([gdal.GCP(lo, la, 0, px, py) for _, px, py, lo, la in GCPS_LL], "EPSG:4277")
print(f"same numbers, so the same residuals: {np.allclose(ll_residuals('EPSG:4326'), ll_residuals('EPSG:4277'))}")
print(f"but the sheet moved {geod.inv(*right, *wrong)[2]:.0f} m")same numbers, so the same residuals: True
but the sheet moved 115 m
Identical residuals, a map 115 m off. Nothing computed from the control points alone can find this; only a comparison with independent data can.
5. Report leave-one-out error, and check against independent data
Leave-one-out error answers the question the residual cannot: how far is a point the fit did not see? An independent check — a feature you did not use as a control point — answers it again, and also catches errors shared by every point. On the Oakham sheet the honest figure is a leave-one-out RMS of 17.5 m, not a residual of 10.2 m — and the three church checks on neighbouring sheets in the mosaic guide agree with it.
Residuals in pixels, and an automatic warning
The data used on this page
The Ordnance Survey six-inch sheet Rutland 009-NW (1884–1891) from Wikimedia Commons (public domain), with seven control points matched to OpenStreetMap (Geofabrik, 1 January 2026; © OpenStreetMap contributors, ODbL) in the control-point guide.
import itertools
from pathlib import Path
import numpy as np
import requests
from osgeo import gdal
from pyproj import Geod, Transformer
gdal.UseExceptions()
Path("data").mkdir(exist_ok=True)
UA = {"User-Agent": "spatialworkflow.io guide (https://www.spatialworkflow.io)"}
if not Path("data/os_009-NW.jpg").exists():
Path("data/os_009-NW.jpg").write_bytes(requests.get(
"https://upload.wikimedia.org/wikipedia/commons/a/a0/"
"Map_of_Rutland_OS_Map_name_009-NW%2C_Ordnance_Survey%2C_1884-1891.jpg", headers=UA, timeout=300).content)
GCPS = [ # (id, pixel, line, easting, northing)
("oakham_church", 3141, 487, 486065.6, 308906.1),
("egleton_church", 4824, 2031, 487598.5, 307540.3),
("brooke_rd_crossing", 2956, 1159, 485892.0, 308306.8),
("egleton_level_crossing", 3998, 2481, 486838.1, 307119.0),
("stamford_rd_guide_post", 5127, 988, 487872.3, 308486.1),
("braunston_rd_bridleway", 1263, 1675, 484382.4, 307839.2),
("brooke_guide_post", 1547, 3292, 484648.2, 306376.1),
]
to_ll = Transformer.from_crs(27700, 4326, always_xy=True)
GCPS_LL = [(k, px, py, *to_ll.transform(e, n)) for k, px, py, e, n in GCPS]
geod = Geod(ellps="WGS84")
def fit(points):
X = np.array([[1, g[1], g[2]] for g in points], float)
return np.linalg.lstsq(X, np.array([[g[3], g[4]] for g in points]), rcond=None)[0]
def errors(A, points):
return [float(np.hypot(*(np.array([1, g[1], g[2]]) @ A - (g[3], g[4])))) for g in points]
def rms(v):
return float(np.sqrt(np.mean(np.square(v))))
def loo(points):
return [errors(fit(points[:i] + points[i + 1:]), [g])[0] for i, g in enumerate(points)]
def gdal_errors(fit_points, test_points, opts):
v = gdal.Translate("", "data/os_009-NW.jpg", format="VRT", bandList=[1], outputSRS="EPSG:27700",
GCPs=[gdal.GCP(e, n, 0, c, r) for _, c, r, e, n in fit_points])
t = gdal.Transformer(v, None, opts)
return [float(np.hypot(*(np.array(t.TransformPoint(0, c, r)[1][:2]) - (e, n)))) for _, c, r, e, n in test_points]
def warp_centre(gcps, srs):
v = gdal.Translate("", "data/os_009-NW.jpg", format="VRT", bandList=[1], outputSRS=srs, GCPs=gcps)
w = gdal.Warp("", v, format="VRT", dstSRS="EPSG:4326", polynomialOrder=1)
gt = w.GetGeoTransform()
return gt[0] + gt[1] * w.RasterXSize / 2, gt[3] + gt[5] * w.RasterYSize / 2
def ll_residuals(srs):
"""Residuals of an affine fit in the declared CRS's own units: the numbers are the same either way."""
X = np.array([[1, g[1], g[2]] for g in GCPS_LL], float)
W = np.array([[g[3], g[4]] for g in GCPS_LL])
return X @ np.linalg.lstsq(X, W, rcond=None)[0] - W
LOW_FOUR = [list(s) for s in itertools.combinations(GCPS, 4) if rms(errors(fit(list(s)), list(s))) < 2.0]
print(len(GCPS), "control points")7 control points
Example 1: residuals in pixels, the way QGIS reports them
The QGIS georeferencer lists residuals in pixels of the scan. Converted with the fitted scale, they are the same numbers:
scale = np.hypot(*fit(GCPS)[1]) # metres per pixel
for (name, *_), r in zip(GCPS, errors(fit(GCPS), GCPS)):
print(f"{name:24} {r / scale:5.1f} px = {r:5.1f} m")oakham_church 9.8 px = 8.8 m
egleton_church 5.7 px = 5.1 m
brooke_rd_crossing 12.4 px = 11.2 m
egleton_level_crossing 16.2 px = 14.6 m
stamford_rd_guide_post 13.8 px = 12.5 m
braunston_rd_bridleway 10.6 px = 9.6 m
brooke_guide_post 6.3 px = 5.7 m
A georeferencer's "mean error" is a summary of these. Whatever the units, it is a residual, with all the blind spots on this page.
Example 2: an automatic warning
A simple ratio flags the fits to distrust: leave-one-out RMS divided by residual RMS. For a reasonable affine fit it is a little above 1; far above that, the residuals are flattering the fit:
for label, pts in (("all 7", GCPS), ("4 of 7", LOW_FOUR[0]), ("4 of 7", LOW_FOUR[1])):
ratio = rms(loo(pts)) / max(rms(errors(fit(pts), pts)), 1e-9)
print(f"{label:6} leave-one-out / residual = {ratio:5.1f}")all 7 leave-one-out / residual = 1.7
4 of 7 leave-one-out / residual = 4.8
4 of 7 leave-one-out / residual = 4.3
With only four points, removing one leaves three — an exactly determined affine fit — so the leave-one-out error is driven by the geometry of the remaining three, and the ratio jumps. The warning is the point: four points cannot check each other.
Fitting and testing are different
Fitting a transformation and testing it are different questions answered with the same data. A residual tells you how consistent your control points are with the transformation you chose; it gets smaller whenever the transformation gains freedom or the points lose the ability to contradict each other. Neither makes the map better.
Leave-one-out error treats each point in turn as a check, fitted without it. It goes up when the transformation has too much freedom (it starts honouring individual errors), and it goes up when there are too few points to constrain the fit away from them. That is why on the Oakham sheet it ranks the three transformations in the opposite order from their residuals, and why it is the number that should be reported.
Some errors escape both, because they are shared by every point: a wrong datum, a misread grid, a consistently mis-identified kind of feature. A transformation absorbs a shared error completely. Only data that did not go into the fit — an independent feature, a neighbouring sheet, a trusted overlay — can expose it.
Three points, splines and removing the worst
- Three points have zero residual for an affine fit and cannot be checked at all.
- A thin-plate spline always has zero residual; never compare its residuals with another method's.
- Removing the worst-residual point usually lowers the RMS and can make the map worse; check leave-one-out before and after.
- Pixel or map units do not change the problem; a residual is a residual.
- A constant offset in every point is invisible to residuals and leave-one-out alike.
- Report the method with the number: "leave-one-out RMS 17.5 m, 7 points, affine".
Internal links
- Ground control points explained — spread, number and leverage
- How to compare georeferencing transformations — leave-one-out in practice
- Georeferenced map lands in the wrong place — the shared datum error at full size
- How accurate is a historical map? — what the honest number means
- Old map is offset from modern data — residuals cannot see an offset
FAQ
Why are my georeferencing residuals low but the map wrong?
Residuals are measured at the points the fit was made from. Fewer points or a more flexible transformation lower them without improving the map, and an error shared by every point never appears in them.
What should I use instead of residual RMS?
Leave-one-out error — predict each control point from the others — and, where possible, independent check points not used in the fit.
Is an RMS error of under one pixel good?
Only if there are enough points to check each other and the transformation is not flexible enough to pass through them all. A spline reports zero whatever the map looks like.
Should I delete the control point with the largest residual?
Not on that basis alone. Check its leave-one-out error and look at the scan; the point with the largest residual is often the one holding the edge of the sheet in place.