Applying Rubber-Sheeting Patches with GDAL

When a scanned map or an older survey layer sits a few metres off its true position and the offset is not constant, a rubber-sheeting patch built from ground control points corrects it — part of automated patching for minor geometry shifts within collaborative editing workflows.

Concept & Context

Not every misalignment is a coordinate reference system problem. A layer can carry the right EPSG:4326 tag, project cleanly, and still drift because the underlying source was distorted before anyone digitized it — paper stretch on a scanned cadastral sheet, uneven georeferencing of an orthophoto, or accumulated digitizing error that grows toward one edge. These offsets are local and non-uniform: two metres north-east in one corner, half a metre south-west in another. A global reprojection cannot remove them because there is no single formula that describes the distortion.

Rubber-sheeting solves this by treating the layer as an elastic sheet pinned at a set of matched ground control points (GCPs). Each GCP pairs a coordinate in the distorted source with the coordinate where that feature actually belongs, taken from an authoritative reference. GDAL then stretches everything between the pins with a thin-plate spline — the mathematical analogue of bending a thin metal sheet through fixed heights with minimum bending energy. Because rubber-sheeting is a targeted correction rather than a wholesale rewrite, it fits naturally alongside attribute reconciliation for tabular spatial data as a class of small, reviewable patches that teams apply and commit like any other change. This guide keeps the raster and its companion vectors aligned by driving both from one shared GCP table.

Core Rubber-Sheeting Pipeline

  1. Collect matched control points. Identify at least six unambiguous features that appear in both the distorted layer and the reference — road intersections, building corners, survey monuments. Record the from coordinate (source) and the to coordinate (reference) for each. Distribute them across the whole extent and concentrate a few where the distortion is worst.
  2. Attach GCPs to the raster. Embed the control points into the raster header with gdal_translate -gcp. This step is lossless — no pixels are resampled yet; it only records the tie points and the target CRS.
  3. Warp the raster with a thin-plate spline. Run gdalwarp -tps to elastically resample the raster so each control point lands on its reference position. Use -r bilinear for continuous imagery or -r near for categorical rasters.
  4. Transform the companion vectors. Feed the identical GCP set through a thin-plate spline in Python (osgeo.gdal.Transformer) and push every vertex of the vector layer through it, so parcels and roads move with the raster.
  5. Check residuals and commit. Confirm the control-point residuals are within tolerance, sanity-check a few known features, then commit the patched outputs and the GCP manifest so the correction is reproducible and auditable.

Working Implementation

The script below drives the whole patch from a single CSV of control points (pixel,line,x,y for the raster tie points plus src_x,src_y,dst_x,dst_y map pairs for the vectors). It shells out to GDAL for the raster warp and reuses the same map-space GCPs to build a thin-plate spline Transformer for the vector layer, guaranteeing both products share one displacement field.

#!/usr/bin/env python3
"""
Rubber-sheeting patch: warp a raster and its companion vectors from one GCP set.

Usage:
    python rubber_sheet.py gcps.csv source.tif source.gpkg patched/

CSV columns (header row required):
    pixel,line,map_x,map_y   -> raster tie points (image space -> reference map)
    src_x,src_y,dst_x,dst_y  -> vector control pairs (distorted map -> reference map)
Rows use only the raster columns OR only the vector columns; blanks are ignored.
"""

import csv
import subprocess
import sys
from pathlib import Path

from osgeo import gdal, ogr

gdal.UseExceptions()
TARGET_SRS = "EPSG:4326"


def read_gcps(csv_path: Path):
    """Split the control-point CSV into raster tie points and vector pairs."""
    raster_gcps, vector_pairs = [], []
    with csv_path.open(newline="") as fh:
        for row in csv.DictReader(fh):
            if row.get("pixel"):
                raster_gcps.append(gdal.GCP(
                    float(row["map_x"]), float(row["map_y"]), 0.0,
                    float(row["pixel"]), float(row["line"]),
                ))
            if row.get("src_x"):
                vector_pairs.append((
                    float(row["src_x"]), float(row["src_y"]),
                    float(row["dst_x"]), float(row["dst_y"]),
                ))
    if len(raster_gcps) < 3 or len(vector_pairs) < 3:
        raise ValueError("Need at least 3 raster GCPs and 3 vector pairs for a TPS.")
    return raster_gcps, vector_pairs


def warp_raster(src_tif: Path, gcps, out_tif: Path) -> None:
    """Attach GCPs then thin-plate-spline warp the raster onto the reference frame."""
    # Step 2: embed GCPs without resampling (lossless header edit).
    tagged = gdal.Translate(str(out_tif.with_suffix(".gcp.vrt")), str(src_tif),
                            outputSRS=TARGET_SRS, GCPs=gcps)
    tagged = None  # flush the intermediate VRT
    # Step 3: elastic warp. -tps selects the thin-plate spline; bilinear for imagery.
    subprocess.run([
        "gdalwarp", "-overwrite", "-tps", "-r", "bilinear",
        "-t_srs", TARGET_SRS, "-co", "COMPRESS=DEFLATE",
        str(out_tif.with_suffix(".gcp.vrt")), str(out_tif),
    ], check=True)
    out_tif.with_suffix(".gcp.vrt").unlink(missing_ok=True)


def build_vector_transformer(pairs):
    """Build a TPS transformer from map-space control pairs (distorted -> reference)."""
    gcps = [gdal.GCP(dst_x, dst_y, 0.0, src_x, src_y)  # 'pixel/line' hold source map coords
            for src_x, src_y, dst_x, dst_y in pairs]
    ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
    ds.SetGCPs(gcps, gdal.osr.SpatialReference().ExportToWkt() or "")
    # method=2 forces a thin-plate spline instead of a polynomial fit.
    return gdal.Transformer(ds, None, ["METHOD=GCP_TPS"])


def warp_vectors(src_gpkg: Path, tps, out_gpkg: Path) -> None:
    """Push every vertex of every geometry through the shared TPS transform."""
    src = ogr.Open(str(src_gpkg))
    drv = ogr.GetDriverByName("GPKG")
    out_gpkg.unlink(missing_ok=True)
    dst = drv.CreateDataSource(str(out_gpkg))
    for i in range(src.GetLayerCount()):
        layer = src.GetLayer(i)
        out_layer = dst.CreateLayer(layer.GetName(), geom_type=layer.GetGeomType())
        out_layer.CreateFields(layer.schema)
        for feat in layer:
            geom = feat.GetGeometryRef().Clone()
            for gi in range(geom.GetGeometryCount() or 1):
                ring = geom.GetGeometryRef(gi) if geom.GetGeometryCount() else geom
                for vi in range(ring.GetPointCount()):
                    x, y, _ = ring.GetPoint(vi)
                    ok, (nx, ny, _) = tps.TransformPoint(0, x, y)  # 0 = forward
                    if not ok:
                        raise RuntimeError(f"TPS failed at vertex {vi} of feature {feat.GetFID()}")
                    ring.SetPoint_2D(vi, nx, ny)
            out_feat = ogr.Feature(out_layer.GetLayerDefn())
            out_feat.SetFrom(feat)
            out_feat.SetGeometry(geom)
            out_layer.CreateFeature(out_feat)
    dst = None


def main(csv_path, src_tif, src_gpkg, out_dir):
    out = Path(out_dir); out.mkdir(parents=True, exist_ok=True)
    raster_gcps, vector_pairs = read_gcps(Path(csv_path))
    warp_raster(Path(src_tif), raster_gcps, out / "patched.tif")
    tps = build_vector_transformer(vector_pairs)
    warp_vectors(Path(src_gpkg), tps, out / "patched.gpkg")
    print(f"Patched raster and vectors written to {out}/")


if __name__ == "__main__":
    if len(sys.argv) != 5:
        print("Usage: python rubber_sheet.py gcps.csv source.tif source.gpkg patched/")
        sys.exit(1)
    main(*sys.argv[1:])

Validation & Output Verification

The value of a rubber-sheeting patch is only as good as its residuals — the leftover distance between each warped control point and its reference target. Ask gdalinfo to report the tie points and inspect the transform metadata:

gdalinfo patched/patched.tif | grep -A2 "GCP\["

For a numeric residual report, run the standalone gdaltransform in inverse mode against each reference coordinate and confirm the round-trip error is sub-pixel:

echo "36.8219 -1.2921" | gdaltransform -tps -i patched/patched.tif
# compare the returned pixel/line against the GCP you supplied for that feature

Confirm the vectors moved with the raster by checking that a known feature now overlays correctly. A quick ogrinfo extent comparison flags gross errors — the patched extent should shift by the same order of magnitude as your measured offset, not by kilometres:

ogrinfo -so -al patched/patched.gpkg | grep Extent

Finally, commit the patch as a reviewable unit. Version the GCP CSV alongside the outputs so the correction is reproducible; anyone can re-run the script and get byte-comparable geometry. Treat the GCP file as the source of truth, exactly as automated conflict detection in merge requests treats validation rules — the patch is data, not a one-off manual nudge.

Failure Modes

  • Warped raster ripples or folds between control points — too many GCPs, or several that disagree with their neighbours, over-constrain the thin-plate spline and force it to bend sharply. Fix: drop the noisiest points, keep six to twelve well-distributed pins, and re-run; a smoother point set yields a smoother surface.
  • Vectors drift away from the raster after the patch — the raster and vector transforms were built from different GCP sets or opposite src/dst orientations. Fix: drive both from one CSV and confirm the vector pairs are ordered distorted→reference to match the raster tie points.
  • gdalwarp -tps produces a hugely oversized output — a single mislabelled control point pulls one corner far off, ballooning the output bounding box. Fix: inspect residuals with gdaltransform, find the outlier with an error orders of magnitude larger than the rest, and correct or remove it.
  • Sub-pixel offset remains after warping — the distortion is actually a uniform datum or CRS shift, not local warping. Fix: stop rubber-sheeting and reproject instead; a thin-plate spline wastes control points modelling a shift a single transform removes cleanly.

Back to Automated Patching for Minor Geometry Shifts