Georeferenced Map Lands in the Wrong Place

Problem statement

You georeferenced a scanned map, the warp ran without an error, and the result is somewhere it should not be: in the sea off Ireland, in the Indian Ocean, at a coordinate of 485,649 degrees, or almost right but a hundred metres off. GDAL rarely complains about any of these. It fits whatever numbers it is given, and every one of the mistakes below produces a valid-looking GeoTIFF.

This page reproduces each one with the same seven control points on the Ordnance Survey six-inch sheet of Oakham (1884–1891), measures where it puts the map, and shows the check that catches it.

Quick answer

Compute where the warped sheet's centre lands and compare it with where you know the map should be. The size of the miss identifies the mistake:

for label, (gcps, srs) in CASES.items():
    try:
        lon, lat = centre(gcps, srs)
        miss = geod.inv(TRUE[0], TRUE[1], lon, lat)[2] / 1000 if abs(lat) <= 90 and abs(lon) <= 180 else float("nan")
        print(f"{label:34} centre {lon:12.4f} {lat:10.4f}   " + (f"{miss:9.3f} km off" if miss == miss else "not a coordinate"))
    except RuntimeError as err:
        print(f"{label:34} GDAL error: {str(err)[:60]}")

correct                            centre      -0.7351    52.6606       0.000 km off
pixel and world swapped            centre     -15.6838    52.1319    1017.425 km off
latitude and longitude swapped     centre      52.6606    -0.7351    7710.586 km off
metres declared as degrees         centre  485649.8698 307741.3027   not a coordinate
WGS 84 declared as OSGB 1936       centre      -0.7367    52.6610       0.116 km off

Argument order, axis order, the declared CRS

1. Check the order of the arguments to gdal.GCP

The constructor is gdal.GCP(x, y, z, pixel, line) — world first. Control points are usually picked pixel first, so swapping the two is the commonest mistake of all:

right = gdal.GCP(486065.6, 308906.1, 0, 3141, 487)
wrong = gdal.GCP(3141, 487, 0, 486065.6, 308906.1)
for g in (right, wrong):
    print(f"world ({g.GCPX:,.1f}, {g.GCPY:,.1f})  pixel ({g.GCPPixel:,.1f}, {g.GCPLine:,.1f})")

world (486,065.6, 308,906.1)  pixel (3,141.0, 487.0)
world (3,141.0, 487.0)  pixel (486,065.6, 308,906.1)

With the roles reversed, GDAL fits the inverse of the transformation you meant and the sheet lands about a thousand kilometres away — west of Ireland for this sheet. Print one GCP back in the form above and the mistake is obvious.

2. Check the axis order of geographic coordinates

Control points from a web map, a GPS or OpenStreetMap usually arrive as longitude and latitude — but a spreadsheet column order, or a library that follows the EPSG axis order, can hand you latitude first. GDAL's GCPs are always x then y, longitude then latitude:

lon, lat = to_wgs84.transform(486065.6, 308906.1)
print(f"Oakham church: lon {lon:.5f}, lat {lat:.5f}")
print("as GCP x, y:", (round(lon, 5), round(lat, 5)), "- swapped:", (round(lat, 5), round(lon, 5)))

Oakham church: lon -0.72865, lat 52.67104
as GCP x, y: (-0.72865, 52.67104) - swapped: (52.67104, -0.72865)

Latitude 52.66 read as longitude and longitude −0.73 read as latitude puts the map just south of the equator in the Indian Ocean, off Somalia — the miss in the quick answer.

3. Check the CRS you declared matches the numbers

outputSRS on the translate step does not convert anything; it labels the GCP coordinates. Declare British National Grid metres as EPSG:4326 and GDAL accepts them as degrees:

v = gdal.Translate("", "data/os_009-NW.jpg", format="VRT", bandList=[1], outputSRS="EPSG:4326",
                   GCPs=[gdal.GCP(e, n, 0, px, py) for _, px, py, e, n in GCPS])
print(gdal.Warp("", v, format="VRT", dstSRS="EPSG:4326", polynomialOrder=1).GetGeoTransform()[:4])

(483228.1353661027, 0.9011104928400812, 0.0, 309370.0598893213)

An origin at 483,000 "degrees" is not a place, and no error says so. Any coordinate outside ±180, ±90 in a geographic CRS is this mistake.

4. Check the datum when the miss is about 100 m

The subtle one. Longitude and latitude from GPS or OpenStreetMap are WGS 84. Old British maps and many British datasets use OSGB 1936. The numbers look alike — both are degrees around 52.7° N, 0.7° W — but declaring WGS 84 values as EPSG:4277 (OSGB 1936 geographic) moves the map:

a = centre([gdal.GCP(lo, la, 0, px, py) for _, px, py, lo, la in GCPS_LL], "EPSG:4326")
b = centre([gdal.GCP(lo, la, 0, px, py) for _, px, py, lo, la in GCPS_LL], "EPSG:4277")
print(f"same numbers, two datums: {geod.inv(*a, *b)[2]:.0f} m apart")

same numbers, two datums: 115 m apart

About a hundred metres, mostly east–west — large enough to put roads in the wrong fields, small enough to be mistaken for the old map's own error. A consistent shift of that size across the whole sheet is a datum, not the surveyor.

A check to run after every fit

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.

from pathlib import Path

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_wgs84 = Transformer.from_crs(27700, 4326, always_xy=True)
GCPS_LL = [(k, px, py, *to_wgs84.transform(e, n)) for k, px, py, e, n in GCPS]   # lon, lat in WGS 84
geod = Geod(ellps="WGS84")

def centre(gcps, srs):
    """Longitude and latitude of the centre of the sheet after an affine warp to EPSG:4326."""
    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

CASES = {
    "correct":                        ([gdal.GCP(e, n, 0, px, py) for _, px, py, e, n in GCPS], "EPSG:27700"),
    "pixel and world swapped":        ([gdal.GCP(px, py, 0, e, n) for _, px, py, e, n in GCPS], "EPSG:27700"),
    "latitude and longitude swapped": ([gdal.GCP(la, lo, 0, px, py) for _, px, py, lo, la in GCPS_LL], "EPSG:4326"),
    "metres declared as degrees":     ([gdal.GCP(e, n, 0, px, py) for _, px, py, e, n in GCPS], "EPSG:4326"),
    "WGS 84 declared as OSGB 1936":   ([gdal.GCP(lo, la, 0, px, py) for _, px, py, lo, la in GCPS_LL], "EPSG:4277"),
}
TRUE = centre(*CASES["correct"])
print("sheet centre:", [round(v, 4) for v in TRUE])

sheet centre: [-0.7351, 52.6606]

Example 1: a check to run after every georeferencing

Four assertions catch three of the four mistakes before the file is used — the sheet should be in roughly the right place, on land, the right size, and not mirrored:

gdal.PushErrorHandler("CPLQuietErrorHandler")       # the bad cases emit thousands of PROJ warnings

def check(gcps, srs, expect=(-0.73, 52.66), km=5):
    lon, lat = centre(gcps, srs)
    if not (-180 <= lon <= 180 and -90 <= lat <= 90):
        return ["centre is not a coordinate"]
    problems = []
    if geod.inv(*expect, lon, lat)[2] > km * 1000:
        problems.append(f"centre {geod.inv(*expect, lon, lat)[2] / 1000:,.0f} km from where it should be")
    v = gdal.Translate("", "data/os_009-NW.jpg", format="VRT", bandList=[1], outputSRS=srs, GCPs=gcps)
    try:
        w = gdal.Warp("", v, format="VRT", dstSRS="EPSG:27700", polynomialOrder=1)
    except RuntimeError:
        return problems + ["cannot be projected to British National Grid"]
    gt = w.GetGeoTransform()
    width_km = w.RasterXSize * gt[1] / 1000
    if not 3 < width_km < 7:
        problems.append(f"sheet is {width_km:,.1f} km wide, expected about 4.8")
    if gt[1] < 0 or gt[5] > 0:
        problems.append("mirrored")
    return problems or ["ok"]

for label, (gcps, srs) in CASES.items():
    print(f"{label:34} {'; '.join(check(gcps, srs))}")

correct                            ok
pixel and world swapped            centre 1,018 km from where it should be
latitude and longitude swapped     centre 7,710 km from where it should be
metres declared as degrees         centre is not a coordinate
WGS 84 declared as OSGB 1936       ok

The datum mistake passes the 5 km test, as it would pass any test that does not know the answer to within 100 m. Catching it needs an overlay against data you trust, or knowing which datum each coordinate came from.

Example 2: reading the residual table for clues

The swapped-coordinates fit leaves a tell-tale in its scale. The fitted units per pixel should be the scan's ground resolution — about 0.9 m here:

import numpy as np

for label in ("correct", "pixel and world swapped"):
    gcps, _ = CASES[label]
    P = np.array([[1, g.GCPPixel, g.GCPLine] for g in gcps])
    W = np.array([[g.GCPX, g.GCPY] for g in gcps])
    A = np.linalg.lstsq(P, W, rcond=None)[0]
    print(f"{label:24} {np.hypot(*A[1]):10.4f} units per pixel")

correct                      0.9028 units per pixel
pixel and world swapped      1.1076 units per pixel

1.108 is exactly the reciprocal of 0.903: the swapped fit is mapping metres onto pixels, not pixels onto metres. A scale that is the inverse of the scan's resolution is the signature of this mistake.

Arithmetic fits any numbers you give it

A georeferencing fit is pure arithmetic: it finds the transformation that best maps one list of numbers onto another. It has no idea which list is pixels, which is metres, which is degrees, or which ellipsoid the degrees were measured on. Every mistake on this page gives the arithmetic consistent input, so it produces a consistent — and wrong — output.

The sizes of the misses are characteristic because they come from the scale of the numbers involved. Swapping pixels (thousands) with metres (hundreds of thousands) inverts the transformation and throws the sheet far out of the country. Swapping latitude and longitude reflects the point across the line where they are equal, which for Britain means near the equator. Declaring metres as degrees produces values no geographic CRS can contain. And confusing two datums for the same coordinates moves things by the distance between the datums — around 100 m in Britain between OSGB 1936 and WGS 84.

That makes the size of the miss the first diagnostic. Measure it before changing anything.

GCP order, labels and datum shifts

• gdal.GCP(x, y, z, pixel, line) is world first; the command line's -gcp pixel line x y is pixel first.
• EPSG:4326 axis order is latitude-longitude in the standard but longitude-latitude in GDAL's GCPs; pyproj's always_xy=True keeps it x, y.
• outputSRS labels, it does not convert.
• OSGB 1936 to WGS 84 is about 100 m in Rutland; the exact shift varies across Britain and needs a grid (OSTN15) for metre accuracy.
• No CRS at all leaves the coordinates as they are, and QGIS will put them in whatever CRS the project uses.
• Rows down, northings up. A mirrored sheet usually means row numbers were counted from the bottom.

Internal links

• How to georeference a scanned map with GDAL — the correct recipe
• Old map is offset from modern data — when the miss is tens of metres
• Georeferencing explained — pixels, coordinates and the fit
• How to pick ground control points — where the seven points came from
• Low GCP residuals, but the map is still wrong — when the fit looks good and is not

FAQ

Why is my georeferenced map in the ocean?

Most often because pixel and world coordinates are swapped in the control points, or latitude and longitude are. Measure how far off it is: about 1,000 km for the first, thousands of kilometres near the equator for the second on a British map.

Why are my map's coordinates huge numbers of degrees?

Projected coordinates in metres were declared as a geographic CRS such as EPSG:4326. Declare the CRS the numbers are actually in.

Why is my georeferenced map about 100 m off?

The control points' datum is not the one declared — typically WGS 84 longitude and latitude labelled as OSGB 1936, or the reverse.

What order does gdal.GCP take?

gdal.GCP(x, y, z, pixel, line): world coordinates first, then the pixel position.