Warped Map Is Folded or Stretched

Problem statement

The georeferenced map is in the right place, but part of it is wrong in a way no offset explains: a corner smeared out to a long streak, a field drawn twice, roads that run backwards for a stretch, a margin swollen to several times its size. GDAL reported no error. What happened is that the transformation folded — somewhere on the sheet it maps two different scan pixels to the same place, or turns part of the map over — or stretched one region far beyond the rest.

This page produces both on real maps, measures how much of the sheet is affected with one test that works for any transformation, and shows the fix for each cause.

Quick answer

A transformation folds wherever its Jacobian determinant changes sign. Push a fine grid of scan pixels through it, measure the signed area of every cell, and count the cells whose sign differs from the rest:

for label, (points, scan, size, opts) in CASES.items():
    folded, stretch = fold_check(points, scan, size, opts)
    print(f"{label:44} folded {folded:6.1%}   99th-percentile cell area {stretch:5.2f}x the median")

OS sheet, affine, 7 points                   folded   0.0%   99th-percentile cell area  1.00x the median
OS sheet, thin-plate spline, 7 points        folded   0.0%   99th-percentile cell area  1.03x the median
OS sheet, spline + one conflicting click     folded   5.3%   99th-percentile cell area  1.60x the median
Speed page, polynomial 2, 33 villages        folded   0.6%   99th-percentile cell area  2.21x the median
Speed page, polynomial 3, 33 villages        folded  38.0%   99th-percentile cell area  5.14x the median

Zero folded cells and a stretch near 1 is a healthy warp. Two things produced the bad ones here: a thin-plate spline given one conflicting control point, and a high-order polynomial asked to cover a whole atlas page from points in its middle.

Measure folding, find conflicting points, lower the order

1. Measure folding instead of looking for it

A fold can hide in a corner nobody looks at. The signed area of each cell of a pixel grid, after transformation, is positive for the map the right way up and negative where it has been turned over:

P, folded = pixel_grid(transformer(OS, "data/os_009-NW.jpg", ["METHOD=GCP_TPS"]), (5343, 3573), 90)
print(P.shape[:2], "grid nodes;", int(folded.sum()), "folded cells with the seven good points")

(60, 90) grid nodes; 0 folded cells with the seven good points

2. Look for conflicting control points

A thin-plate spline passes through every control point exactly. Two points a few pixels apart on the scan with world coordinates far apart — a double click on the same church with the second click's coordinates mistyped — force it to stretch a small region to cover a large distance, and it folds to do so:

def conflicts(points, px_radius=200, ratio=3):
    """Pairs of points whose world distance is far from what the median scale predicts."""
    P = np.array([[c, r] for _, c, r, _, _ in points], float)
    W = np.array([[e, n] for *_, e, n in points], float)
    d_px = np.hypot(*(P[:, None] - P[None]).transpose(2, 0, 1))
    d_w = np.hypot(*(W[:, None] - W[None]).transpose(2, 0, 1))
    scale = np.median(d_w[d_px > 0] / d_px[d_px > 0])
    out = []
    for i, j in zip(*np.triu_indices(len(points), 1)):
        if d_px[i, j] < px_radius and d_w[i, j] / max(d_px[i, j] * scale, 1e-9) > ratio:
            out.append((points[i][0], points[j][0], round(d_px[i, j]), round(d_w[i, j])))
    return out

print(conflicts(OS + [CONFLICT]))

[('oakham_church', 'oakham_church_again', 50, 192)]

Fifty pixels on the scan is about 45 m on the ground; the two points claim 192 m. Remove one and the fold disappears.

3. Lower the order, or limit the area

A polynomial of order 3 has 20 coefficients and bends freely everywhere. Fitted to villages in the middle of Speed's county map, it does reasonably inside the villages and wildly beyond them:

hull = Polygon(np.array([[c, r] for _, c, r, _, _ in SPEED], float)[ConvexHull([[c, r] for _, c, r, _, _ in SPEED]).vertices])
for order in (1, 2, 3):
    t = transformer(SPEED, "data/speed_rutland.jpg", ["METHOD=GCP_POLYNOMIAL", f"MAX_GCP_ORDER={order}"])
    P, folded = pixel_grid(t, (12232, 8343), 60)
    gx, gy = np.linspace(0, 12232, 60), np.linspace(0, 8343, folded.shape[0] + 1)
    xs, ys = np.meshgrid((gx[:-1] + gx[1:]) / 2, (gy[:-1] + gy[1:]) / 2)          # cell centres
    inside = np.array([hull.contains(Point(x, y)) for x, y in zip(xs.ravel(), ys.ravel())]).reshape(folded.shape)
    print(f"order {order}: folded on the whole page {folded.mean():5.1%}, inside the villages' hull {folded[inside].mean():5.1%}")

order 1: folded on the whole page  0.0%, inside the villages' hull  0.0%
order 2: folded on the whole page  0.6%, inside the villages' hull  0.0%
order 3: folded on the whole page 37.7%, inside the villages' hull  0.0%

Two fixes, both needed: a lower order where the data cannot support a higher one, and a cutline so the warp covers only the part of the sheet the points describe (see the early county map guide).

4. Check stretch as well as folds

A transformation can stay unfolded and still distort a region badly. The ratio of each cell's area to the median cell shows it:

for opts, label in ((["MAX_GCP_ORDER=2"], "polynomial 2, all 7 points"), (["MAX_GCP_ORDER=2"], "polynomial 2, 6 points")):
    pts = OS if "7" in label else OS[1:]
    print(label, "-> 99th percentile stretch", round(fold_check(pts, "data/os_009-NW.jpg", (5343, 3573), opts)[1], 2))

polynomial 2, all 7 points -> 99th percentile stretch 1.07
polynomial 2, 6 points -> 99th percentile stretch 1.8

Dropping the one point near the top edge leaves a polynomial with exactly as many coefficients as constraints; the corners of the sheet stretch accordingly. A stretch well above 1.1 on a sheet that should be uniform is a warning even without a fold.

5. Re-warp and re-check

After removing the conflicting point, lowering the order or adding a cutline, run fold_check again. A clean warp has no folded cells and a stretch close to 1 over the area you keep.

What a fold does, and catching it in a pipeline

The data used on this page

The Ordnance Survey six-inch sheet Rutland 009-NW (1884–1891) and John Speed's Rutlandshire (1676 edition of the 1611 map, scanned by the David Rumsey Map Collection, Stanford Libraries), both from Wikimedia Commons (public domain); OpenStreetMap from Geofabrik (1 January 2026; © OpenStreetMap contributors, ODbL). Control points as in the control-point guide and the early county map guide.

from pathlib import Path

import geopandas as gpd
import numpy as np
import requests
from osgeo import gdal
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon

gdal.UseExceptions()
Path("data").mkdir(exist_ok=True)
UA = {"User-Agent": "spatialworkflow.io guide (https://www.spatialworkflow.io)"}
downloads = {
    "data/os_009-NW.jpg": "https://upload.wikimedia.org/wikipedia/commons/a/a0/"
                          "Map_of_Rutland_OS_Map_name_009-NW%2C_Ordnance_Survey%2C_1884-1891.jpg",
    "data/speed_rutland.jpg": "https://upload.wikimedia.org/wikipedia/commons/9/9c/Rutlandshire._%28insets%29_Oukham."
                              "_Stanford._%28IA_dr_rutlandshire-insets-oukham-stanford-12058072%29.jpg",
    "data/rutland.shp.zip": "https://download.geofabrik.de/europe/united-kingdom/england/rutland-260101-free.shp.zip",
}
for target, url in downloads.items():
    if not Path(target).exists():
        Path(target).write_bytes(requests.get(url, headers=UA, timeout=600).content)

OS = [  # (id, pixel, line, easting, northing) on six-inch sheet 009-NW
    ("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),
]
CONFLICT = ("oakham_church_again", 3181, 517, 485915.6, 309026.1)    # a second click, mistyped coordinates
SPEED_PX = {  # modern village name: (column, row) of Speed's church symbol
    "Oakham": (4960, 4035), "Langham": (4475, 3725), "Barleythorpe": (4755, 3920), "Egleton": (5265, 4290),
    "Braunston in Rutland": (4470, 4745), "Brooke": (4705, 4790), "Gunthorpe": (5095, 4745), "Manton": (5430, 4935),
    "Greetham": (6200, 2930), "Cottesmore": (5715, 3205), "Exton": (6165, 3500), "Empingham": (6730, 3890),
    "Whitwell": (6245, 4080), "Upper Hambleton": (5795, 4285), "Burley on the Hill": (5490, 3705), "Barrow": (5495, 2965),
    "Pickworth": (7375, 2930), "Tickencote": (7595, 3875), "Lyndon": (6005, 4880), "Edith Weston": (6415, 4555),
    "Wing": (5735, 5245), "Pilton": (6315, 5335), "North Luffenham": (6600, 5110), "South Luffenham": (6705, 5400),
    "Ridlington": (4755, 5425), "Preston": (5335, 5380), "Glaston": (5650, 5690), "Morcott": (6465, 5735),
    "Ayston": (4960, 5765), "Uppingham": (5125, 6080), "Bisbrooke": (5605, 6035), "Seaton": (5960, 6430),
    "Barrowden": (7070, 5985),
}
places = gpd.read_file("/vsizip/data/rutland.shp.zip/gis_osm_places_free_1.shp").to_crs(27700)
places = places[places.fclass.isin(["village", "town", "hamlet"])].drop_duplicates("name").set_index("name")
SPEED = [(k, c, r, places.loc[k].geometry.x, places.loc[k].geometry.y) for k, (c, r) in SPEED_PX.items()]

def transformer(points, scan, opts):
    v = gdal.Translate("", scan, format="VRT", bandList=[1], outputSRS="EPSG:27700",
                       GCPs=[gdal.GCP(e, n, 0, c, r) for _, c, r, e, n in points])
    t = gdal.Transformer(v, None, opts)
    t.keep = v                                 # the transformer needs its dataset alive
    return t

def pixel_grid(t, size, n):
    """World positions of an n-wide grid of scan pixels, and which cells are turned over."""
    w, h = size
    xs, ys = np.meshgrid(np.linspace(0, w, n), np.linspace(0, h, max(2, int(n * h / w))))
    P = np.array(t.TransformPoints(0, [(x, y, 0) for x, y in zip(xs.ravel(), ys.ravel())])[0])[:, :2]
    P = P.reshape(xs.shape + (2,))
    dx, dy = P[:-1, 1:] - P[:-1, :-1], P[1:, :-1] - P[:-1, :-1]
    det = dx[..., 0] * dy[..., 1] - dx[..., 1] * dy[..., 0]
    return P, np.sign(det) != np.sign(np.median(det))

def fold_check(points, scan, size, opts, n=120):
    t = transformer(points, scan, opts)
    P, folded = pixel_grid(t, size, n)
    dx, dy = P[:-1, 1:] - P[:-1, :-1], P[1:, :-1] - P[:-1, :-1]
    area = np.abs(dx[..., 0] * dy[..., 1] - dx[..., 1] * dy[..., 0])
    return float(folded.mean()), float(np.percentile(area / np.median(area), 99))

CASES = {
    "OS sheet, affine, 7 points":                 (OS, "data/os_009-NW.jpg", (5343, 3573), ["MAX_GCP_ORDER=1"]),
    "OS sheet, thin-plate spline, 7 points":      (OS, "data/os_009-NW.jpg", (5343, 3573), ["METHOD=GCP_TPS"]),
    "OS sheet, spline + one conflicting click":   (OS + [CONFLICT], "data/os_009-NW.jpg", (5343, 3573), ["METHOD=GCP_TPS"]),
    "Speed page, polynomial 2, 33 villages":      (SPEED, "data/speed_rutland.jpg", (12232, 8343), ["MAX_GCP_ORDER=2"]),
    "Speed page, polynomial 3, 33 villages":      (SPEED, "data/speed_rutland.jpg", (12232, 8343), ["MAX_GCP_ORDER=3"]),
}
print(len(OS), "OS points,", len(SPEED), "Speed points")

7 OS points, 33 Speed points

Example 1: what the fold does to the image

Warp the OS sheet with and without the conflicting point and compare the result around Oakham church:

for label, pts in (("seven good points", OS), ("plus the conflicting click", OS + [CONFLICT])):
    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 pts])
    w = gdal.Warp("", v, format="MEM", dstSRS="EPSG:27700", tps=True, xRes=2, yRes=2,
                  outputBounds=(485500, 308500, 486600, 309300), dstAlpha=True)
    a = w.ReadAsArray()
    print(f"{label:28} ink share in the window {(a[0][a[1] > 0] < 128).mean():.1%}")

seven good points            ink share in the window 25.0%
plus the conflicting click   ink share in the window 29.3%

The folded warp draws parts of the town twice and stretches others, which changes even a crude statistic like the share of ink in the window. Anything measured from a folded warp is wrong in that region.

Example 2: catching it in a pipeline

The check is cheap enough to run on every georeferencing:

def assert_no_folds(points, scan, size, opts, max_stretch=1.5):
    folded, stretch = fold_check(points, scan, size, opts, n=80)
    if folded > 0 or stretch > max_stretch:
        raise ValueError(f"warp folds {folded:.1%} of the sheet, stretch {stretch:.2f}x")
    return "ok"

for label, args in CASES.items():
    try:
        print(f"{label:44} {assert_no_folds(*args)}")
    except ValueError as err:
        print(f"{label:44} REJECTED: {err}")

OS sheet, affine, 7 points                   ok
OS sheet, thin-plate spline, 7 points        ok
OS sheet, spline + one conflicting click     REJECTED: warp folds 5.5% of the sheet, stretch 1.61x
Speed page, polynomial 2, 33 villages        REJECTED: warp folds 0.6% of the sheet, stretch 2.20x
Speed page, polynomial 3, 33 villages        REJECTED: warp folds 37.8% of the sheet, stretch 5.10x

Speed's second-order polynomial is rejected on the whole page but folds nothing inside the villages' hull: run the check over the area you will keep, after the cutline, not over margins you are about to throw away.

Folds come from the transformation

A georeferencing transformation is a function from the scan's plane to the map's plane. For the warp to make sense it has to be one-to-one over the area being warped: every scan pixel goes to a distinct place, and the map keeps its orientation. Affine transformations are always one-to-one. Polynomials and splines are not guaranteed to be — they bend, and if they bend far enough, they fold.

Two situations push them that far. A spline is forced through every point, so two nearby points with incompatible coordinates leave it no smooth way through: it has to cover a large distance in a small space, and overshoots into a fold. A polynomial of high order is only constrained where there are points; beyond them, its cubic terms grow without limit and turn the sheet over. Both are the transformation doing exactly what it was asked — the request was impossible or under-specified.

The signed-area test catches both, because a fold is precisely a place where the orientation flips. It needs no knowledge of what the map should look like, which makes it suitable for automated checks; the stretch percentile adds a warning for distortion short of folding.

Affine never folds, and high orders extrapolate

• Affine never folds; if an affine warp looks folded, the scan or the CRS is wrong.
• One conflicting point is enough to fold a spline. Check close pairs of points against the overall scale.
• High-order polynomials extrapolate badly. Order 3 needs at least 10 points, spread to the edges.
• Cut the warp to the area the points cover, especially on atlas pages with decorations and margins.
• Stretch without folds is a softer warning; compare with what the map should look like.
• GDAL's own error threshold for inverting transforms does not detect folds; test explicitly.

Internal links

• Georeferencing transformations explained — why polynomials and splines can bend
• How to compare georeferencing transformations — choosing an order that predicts
• How to draw a distortion grid on an old map — the same grid, used for history
• How to georeference an early county map — cutlines on an atlas page
• Ground control points explained — catching the bad point

FAQ

Why is my georeferenced map folded or smeared?

The transformation is not one-to-one over the sheet: usually a thin-plate spline with two close, conflicting control points, or a high-order polynomial extrapolating beyond its points.

How do I detect a fold in a warped image?

Transform a grid of scan pixels, compute the signed area of each cell, and look for cells whose sign differs from the majority.

How do I fix a folded warp?

Remove or correct the conflicting control point, lower the polynomial order, or restrict the warp to the area covered by the points with a cutline.

Can an affine transformation fold a map?

No. It maps straight lines to straight lines and keeps orientation unless its determinant is negative everywhere, which would mirror the whole sheet rather than fold part of it.