A local model has an order

Software that georeferences an image asks for control points and a transformation order: first, second, or third. First order is an affine transformation with six parameters, second is quadratic with twelve, third is cubic with twenty. The usual advice is to pick the order the number of control points can support, which is a statement about the fitting and not about the map.

What the order actually buys is a rate. And the first thing a first-order model cannot represent is the second derivative of the map — which is the quantity this ladder has spent five essays measuring, under the names flexion and skewness.

Why the exponent is one more than the degree #

The map is smooth, so near a point it is a Taylor series: a constant, a linear part, a quadratic part, and so on. A polynomial model of degree d reproduces the first d+1 of those exactly, and its residual over a patch of radius r is governed by the first one it misses, which is of degree d+1 and therefore of size r^(d+1).

An affine model captures the linear part — which is precisely Tissot’s construction, the object Tissot stops at the first derivative is about — so its residual is the size of the quadratic part, which grows as r². A quadratic model captures the quadratic part and leaves the cubic, so r³.

Measured: 2.001, 3.001 and 4.002.

That is a satisfying agreement and it is not the interesting part, because the exponents are forced. What is not forced is the coefficient, and the coefficient is where the ladder’s own quantity enters.

The coefficient, and the claim that failed #

The essay was written expecting a clean result: the affine model’s residual should be the second derivative’s magnitude times a universal constant, since the constant comes from least squares over a disc and not from any map.

It is not universal.

The reason is that they are two different norms of the same object.

The second derivative of a map from a surface to a plane is a six-component thing: two output coordinates, each with three independent second partials. Flexion and skewness pick out its components along and across a direction of travel and are then averaged over directions, which weights those six components one way. The least-squares residual over a disc weights them another way. No single number extracted from the first can predict the second exactly, and the spread of 1.68 is the size of the disagreement between two reasonable weightings.

What survives is the statement worth having. The practitioner’s residual is the ladder’s quantity, times a factor between one fifth and one third — so a projection with twice the flexion needs twice as many control points’ worth of model over the same patch, and the flexion is computable without fitting anything.

What each order is worth in kilometres #

The exponents and the coefficients together answer the question a practitioner actually has: how large a patch can each order carry?

Five kilometres for an affine model at a metre’s tolerance is a small patch, and it is the correct order of magnitude for the practice it describes: an affine georeference of a single air photo is fine, and an affine georeference of a satellite scene is not.

The scaling is the useful part. At a tolerance ε the affine model reaches √ε, the quadratic ε^(1/3) and the cubic ε^(1/4), so tightening the tolerance by a hundred shrinks the affine patch by ten, the quadratic by 4.6 and the cubic by 3.2. A high-order model is not merely better; it degrades more slowly.

The relation to two things this site already computed #

The exponent 2 has appeared twice before in this site under different names, and they are the same statement.

How small is flat enough measures the departure of a spherical cap from its tangent plane and finds it growing as the square of the radius. A tangent plane is an affine model of the sphere, so that measurement is the r² law for a degree-one model, computed geometrically instead of by fitting.

The same law turns up in the practice of surveying, where the smallest scale distortion any map of a circular patch can have grows with the patch’s radius at a slope of exactly two, and ten parts per million is reached at forty kilometres. That is the boundary between plane surveying and geodesy, and it is an affine model’s reach at a surveyor’s tolerance rather than at a mapping one.

A straight segment is a claim about a plane measures the departure of a drawn chord from a geodesic and finds κL²/8. A chord is an affine model of an arc, so that is the same law again.

It turns up a third time in stored geometry. Six east-west segments from 50 to 1,600 kilometres, held as two points and drawn straight, depart from the geodesic at a fitted slope of 2.001, and a segment of fifty kilometres is already forty-nine metres out. Three parts of this subject have measured the $r^2$ law independently; what this essay adds is what the constant in front of it is made of.

Three measurements, three fields, one exponent — and none of them previously connected to the second derivative that produces it.

The same measurement on a different map, and why it looks identical #

A result that came out of one projection would be a result about that projection.

That the exponents are identical and the intercepts are not is the whole structure of the result in one picture. A polynomial model’s rate is a property of polynomials; its constant is a property of the map, and the constant is what this ladder can predict from a quantity computed without any fitting at all.

What was computed, and how #

The source coordinates are local east and north offsets in radians, which is what a georeferencing tool has after it has been told the patch’s centre. The target is the projection’s own plane. Both are exact; no control-point noise is simulated, so the residual measured is the model’s own inadequacy and nothing else.

The fit is ordinary least squares over a ring-structured sample of the patch — twelve rings, with more points on the outer ones so that the sample is roughly uniform by area. The same lstsq the conformal solver uses, with columns scaled to unit norm.

The exponent is fitted across five patch radii spanning a factor of sixteen, which is enough to separate 2 from 3 and 3 from 4 by a wide margin. The measured values sit within 0.002 of their predictions.

The reach is found by bisection on the patch radius, re-fitting at each trial. That is slower than inverting the fitted power law and it is honest about a possibility the power law hides: that the exponent is not constant over four orders of magnitude of tolerance.

Why the residual is reported in map units and in metres #

Two different quantities appear in the figures above and confusing them would make the reach numbers meaningless.

The fit residual is in the projection’s own plane units, which for a unit-sphere projection are radians of arc at the map’s own scale. That is the quantity whose exponent is measured, and it is scale-free: multiplying a projection by a constant multiplies the residual by the same constant and changes no exponent.

The reach is in kilometres on the ground, which requires dividing the residual by the local scale factor and multiplying by the Earth’s radius. That conversion is what makes the number comparable across projections — a map that is drawn twice as large has twice the residual and exactly the same accuracy — and it is where the projection’s first derivative enters a calculation otherwise about its second.

The general rule is the one this site applies to every distortion measure: a quantity that changes when the map is enlarged is not a property of the map. The exponents pass that test by construction and the reach passes it by the conversion.

Two refusals #

Each order must be worth having. A cubic model whose residual matched a quadratic’s at the same patch size would mean the extra terms were doing nothing and the ladder of orders was a fiction. Measured at a 4° patch, each extra degree cuts the residual by more than a factor of five.

The exponent must be the model’s and not the projection’s. The same three exponents come out for every projection tested, which is what makes them a property of polynomial approximation rather than of any map. If Mollweide gave 2.0 and Mercator 2.4, the whole framing would be wrong.

What a practitioner should take from it #

Three things, in the order they matter.

The order is the decision. Between an affine and a quadratic model there is a factor of twelve in reach at a fixed tolerance, and between a quadratic and a cubic a factor of three and a half. No choice of projection moves the answer nearly that far.

The patch size is the other decision. Since the residual goes as r^(d+1), halving a patch is worth as much as adding a degree — more, for a cubic model. Georeferencing a scene in tiles with an affine model per tile is a real alternative to a cubic model over the whole scene, and it has the advantage that the residual can be checked per tile.

The flexion says which projections are hard. A projection with twice the second-order magnitude needs a patch √2 times smaller at the same order and tolerance. That is computable in closed form for every projection in the library, without fitting anything, which is what makes the ladder’s abstract quantity operational.

What the tiled alternative actually costs #

The suggestion that a scene can be georeferenced in tiles with an affine model per tile is worth costing, because the reach numbers make the arithmetic decisive and it does not come out the way the suggestion implies.

On Mercator at a one-metre tolerance the three reaches are 5, 60 and 207 kilometres. Covering the cubic model’s 207-kilometre disc with affine patches of 5 kilometres takes (207/5)² ≈ 1,714 tiles; with quadratic patches of 60 kilometres it takes 12. Each affine tile carries six parameters and needs at least three control points, so the tiled affine scheme is 10,284 parameters and upwards of five thousand control points against the single cubic model’s twenty parameters and ten points.

A factor of five hundred in parameters, for the same tolerance over the same ground. That is the cost of buying reach by subdivision rather than by degree, and it is what the r^(d+1) law says it should be: subdivision buys accuracy as a power of the tile count and degree buys it as an exponential in the degree.

There is a second cost that the parameter count does not show and that decides the matter for an image. Adjacent affine tiles do not agree on their shared edge. Each is fitted independently to its own patch, so at the seam the two models differ by something of the order of the tolerance — a metre, in this example — and the georeferenced scene is discontinuous there. A feature crossing the seam is drawn twice, offset, which is the same failure a tile is drawn without its neighbours prices for a rendering and where two zones meet prices for a national grid, arriving in a georeference.

The single high-order model has no seam by construction. So the honest comparison is not twelve tiles or one model but twelve tiles with eleven seams, or one model, and the seams are visible where the residual is not.

What the tiled scheme does buy is real and worth keeping in view: a residual that can be checked per tile against that tile’s own control points, and a failure that is local rather than global. A cubic model fitted over a whole scene has twenty parameters constrained by every point at once, so a bad control point in one corner moves the model everywhere — which is the noise half of the trade this essay otherwise sets aside, and it is the half that argues for tiles.

The reconciliation is the one the reach numbers point at directly: subdivide by a factor of two or four, not by a factor of forty. Halving a patch buys a factor of four for an affine model, which is one quarter of what a degree buys, and it costs four tiles and three seams rather than 1,714 and thousands. Two levels of subdivision with a quadratic model reaches most of what a cubic does with a residual that is checkable in sixteen places, and that is the configuration the arithmetic actually recommends.

Where the model stops #

No noise, no control points. A real fit is to measured control points with errors in them, and the residual then has two parts: the model’s inadequacy, measured here, and the propagated observation error, which grows with the number of parameters rather than shrinking. That trade — a higher-order model fits the map better and the noise worse — is the actual reason software warns against high orders, and nothing here measures it. What this essay gives is the half of the trade that is a property of the map.

One point per projection. The coefficients above are measured at 20°E 40°N. The second-order magnitude varies across a map by an order of magnitude, so the reach varies with it; the exponents do not.

Polynomials in a plane, not on a sphere. The models fitted here take plane coordinates to plane coordinates, which is what a georeferencing tool does. A model that respected the sphere — a rational function, or a polynomial in a conformal chart — would do better, and the exactly conformal series of the neighbouring ladder is the extreme case of that idea.

A model with a purpose is different from a model with a residual. A rubber-sheet transformation is often chosen to make control points match exactly, which is interpolation rather than approximation and has no convergence order at all.

The same question asked about a sphere rather than a polynomial gives the same shape of answer. The Earth’s Gaussian radius of curvature varies across a band, the global mean radius sits 2,203 parts per million from the ground there, and the best local radius fits with no bias at all — fitting a local model to a global surface is one operation wherever it appears, and it always leaves a residual set by the next derivative up.

The one place this argument meets a real practice head-on #

National grids are the case where a polynomial model of a projection is not a convenience but a published standard.

The transverse Mercator projection has no elementary closed form on an ellipsoid, so every national grid is defined by a truncated series — Krüger’s, to fourth order in most standards. That is a polynomial model of exactly the kind measured here, chosen once for a whole country and held to a tolerance in millimetres.

The difference between that case and a georeference is the variable being expanded in. Krüger’s series is a polynomial in the distance from the central meridian, so its reach is a strip rather than a patch, and the zone system exists to keep every point inside the strip where the series holds — which is the argument UTM and the zone system makes about six-degree zones. A georeference is a polynomial in position within a patch, so its reach is a disc.

Same mathematics, same exponents, different geometry — and in both cases the practical decision is where to cut the world up so that a low-order model is enough.

Who found it, and when #

The polynomial rectification of imagery is as old as aerial photogrammetry and the modern form of it is from the 1970s satellite era, when Landsat scenes had to be registered to maps. The convergence orders are standard approximation theory, older than any of it.

What is not standard is the connection. The georeferencing literature discusses order selection as a statistical question — degrees of freedom, control-point counts, overfitting — and the cartographic literature discusses second-order distortion measures as a property of projections, and the two have no citations in common. Both are describing the coefficient of r^(d+1) in the same expansion.

Goldberg and Gott’s 2007 measures put flexion and skewness into cartography with an argument about how maps look; this rung’s claim is narrower and more practical, which is that the same numbers say how large a patch a linear model can carry.

Where the ladder goes next #

The flexion ladder now runs from the definition of the second derivative to a number a georeferencing tool could use. Five rungs, one quantity, and the last of them measured against a practice that has never mentioned it.

What is left open is not this ladder’s own question: the two halves of a second-order criterion disagree with each other more than the whole disagrees with the first-order one, at Spearman 0.45. That is the largest unexplained number this field has, and a ladder that separated bending from stretching over regions rather than over the whole sphere would be the place to explain it.