Georeferencing Orthomosaics with Ground Control Points

TL;DR: Measure RTK ground control points, pair each world coordinate with its pixel location to build a GDAL GCP list, warp the raster with rasterio.transform.from_gcps plus reproject (or gdal.Warp for a thin-plate spline), then report a per-point and aggregate RMSE residual — and skip the whole exercise when RTK-tagged imagery already meets your accuracy budget.

Why Ground Control Points Arise in Orthomosaic Workflows

A structure-from-motion pipeline (Pix4D, OpenDroneMap, Agisoft Metashape) can stitch a geometrically consistent orthomosaic from image EXIF GPS alone, but the absolute position of that mosaic inherits the drone’s onboard GPS error. A consumer DJI Phantom 4 RGB logs image centres to ±2–5 m horizontally and far worse vertically. Drop that mosaic onto a basemap and the field boundary can sit a full swath off, so a variable-rate prescription clipped to it applies product to the neighbour’s headland. Ground control points (GCPs) — surveyed markers with known RTK coordinates — pin the mosaic to the real world and collapse that error to a couple of centimetres.

The cost of skipping GCPs is not just a constant shift. SfM blocks accumulate a low-frequency “bowl” or “dome” distortion across a flight, so the error is spatially variable: fine near the block centre, tens of centimetres at the edges. A single-shift correction cannot fix that, which is why georeferencing uses a transform fitted to several distributed points. Getting the target CRS right is a prerequisite — read understanding CRS in precision agriculture first, because a GCP list in the wrong datum produces a clean-looking fit that is systematically metres off.

The diagram below shows the flow from field markers to a residual-checked GeoTIFF.

Ground control point georeferencing pipeline Four-stage flow: RTK-measured field markers, pixel-to-world GCP pairing, polynomial or thin-plate-spline warp, and an RMSE-validated georeferenced GeoTIFF. RTK field markers ±1–3 cm world XY Pair pixel ↔ world GDAL GCP list Fit + warp polynomial / TPS Georeferenced GeoTIFF RMSE-checked 1. Measure 2. Match 3. Transform 4. Validate

Prerequisites

Only what differs from the Orthomosaic Stitching Workflows guide:

TEXT 
rasterio==1.3.10
gdal==3.8.4
numpy==1.26.4

Install with:

BASH 
pip install rasterio==1.3.10 numpy==1.26.4
# GDAL Python bindings are easiest via conda to match the system libgdal:
conda install -c conda-forge gdal=3.8.4

Input requirements:

  • An un-georeferenced or coarsely-georeferenced orthomosaic (GeoTIFF, any band count).
  • A GCP table with, per point, the marker world coordinates (easting, northing) in a known projected CRS and the corresponding image location (col, row in pixels).
  • Marker world coordinates from an RTK rover (±1–3 cm). Field GPS (±1–3 m) is not accurate enough to improve an SfM block.
  • All GCP world coordinates in the same projected CRS you intend to write — mixing WGS84 lat/lon and UTM in one list is the most common silent failure.

Step-by-Step

What each level of ground control buys Bars comparing absolute positional error for four georeferencing setups: no ground control with consumer GNSS, no ground control with an on-board RTK receiver, three clustered control points, and five or more points spread across the block with one held back as a check. No ground control, consumer GNSS 1.5–3 m absolute error No ground control, RTK on board 5–10 cm absolute error 3 points, clustered unconstrained away from them 5+ points, spread and one central 2–4 cm, checked by a held-out point Bars are relative error, not to scale in metres. The last row is the only one whose accuracy has actually been measured rather than assumed.

Step 1 — Place and measure the control points

Lay durable, high-contrast targets (checkerboard or chevron aeroplates, ~0.5 m for a 100 m flight so the centre resolves to several pixels) before the flight. Distribute them: the four corners plus the centre catch most of the SfM bowl distortion; add mid-edge points on fields larger than ~20 ha. Measure each target centre with an RTK rover and record the CRS the rover reports (often geographic WGS84, so you will project it). Reserve one or two markers as checkpoints — never fed to the fit — so the accuracy you report is independent.

Step 2 — Match world coordinates to pixel locations

Each control point needs its pixel column and row in the orthomosaic. In production this comes from your SfM report or from clicking targets in QGIS; here we load a CSV. Pixel row/col are (x, y) in GDAL’s GCP object where x is column and y is row.

Step 3 — Build the GDAL GCP list and pick a transform order

A first-order (affine) polynomial handles shift, scale, rotation, and shear — the right default when the mosaic is already internally consistent and you only need to seat it in absolute space. A thin-plate spline (TPS) warps locally to honour every point exactly, which fixes residual SfM doming but will happily contort the raster around a single mislabelled marker. Use the polynomial unless per-point residuals prove you need local warping.

Step 4 — Warp the raster

Two routes. rasterio.transform.from_gcps collapses the GCPs into a single affine Affine transform and is the lean choice for order-1 fits; you then reproject into a fresh array. For TPS or order-2/3 polynomials, drive gdal.Warp with tps=True or polynomialOrder=N, because those transforms are not a single affine and cannot be expressed as one Affine.

Step 5 — Report the RMSE residual budget

The number that certifies the job: for each control point, transform its pixel location forward and compare to its surveyed world coordinate. Report per-point residuals and the aggregate RMSE, and assert the aggregate against a budget tied to your ground sample distance. The complete, directly runnable script:

PYTHON 
import numpy as np
import rasterio
from rasterio.control import GroundControlPoint
from rasterio.transform import from_gcps
from rasterio.warp import reproject, Resampling

# ── Config ────────────────────────────────────────────────────────────────
SRC_PATH = "ortho_raw.tif"
DST_PATH = "ortho_georef.tif"
TARGET_CRS = "EPSG:32615"          # UTM zone 15N — replace with your zone
GSD = 0.03                          # ground sample distance, metres/pixel
RMSE_BUDGET_M = 2 * GSD             # accept if horizontal RMSE < 2 GSD

# GCPs: (col, row) pixel  ->  (easting, northing) in TARGET_CRS, RTK-measured.
# The last entry is held out as an independent checkpoint.
GCP_TABLE = [
    # col,    row,      easting,      northing
    (  120.5,  95.2,  512340.11,  4790120.44),
    ( 8210.0, 110.7,  512610.83,  4790118.02),
    (  140.3, 6050.9,  512338.55,  4789930.71),
    ( 8190.6, 6042.1,  512609.20,  4789931.98),
    ( 4160.2, 3075.4,  512474.66,  4790025.13),
]
CHECKPOINTS = [
    (2100.8, 1520.6, 512406.19, 4790073.55),
]

# ── 1. Build the rasterio/GDAL GCP list ───────────────────────────────────
gcps = [
    GroundControlPoint(row=r, col=c, x=e, y=n)
    for (c, r, e, n) in GCP_TABLE
]

# ── 2. Fit a first-order (affine) transform from the GCPs ──────────────────
transform = from_gcps(gcps)

# ── 3. Warp the source raster into the target grid ─────────────────────────
with rasterio.open(SRC_PATH) as src:
    # Size the destination to the source dimensions; from_gcps preserves scale.
    dst_height, dst_width = src.height, src.width
    profile = src.profile.copy()
    profile.update(
        crs=TARGET_CRS,
        transform=transform,
        driver="GTiff",
        compress="deflate",
    )
    dst_arr = np.zeros((src.count, dst_height, dst_width), dtype=src.dtypes[0])

    for b in range(src.count):
        reproject(
            source=rasterio.band(src, b + 1),
            destination=dst_arr[b],
            src_crs=TARGET_CRS,
            src_transform=transform,
            dst_crs=TARGET_CRS,
            dst_transform=transform,
            resampling=Resampling.bilinear,
        )

with rasterio.open(DST_PATH, "w", **profile) as dst:
    dst.write(dst_arr)

# ── 4. Residual (RMSE) reporting ───────────────────────────────────────────
def forward(col, row, aff):
    """Pixel (col,row) -> world (x,y) via the fitted affine transform."""
    x, y = aff * (col, row)
    return x, y

def residuals(rows, aff):
    out = []
    for c, r, e, n in rows:
        px, py = forward(c, r, aff)
        out.append(np.hypot(px - e, py - n))
    return np.array(out)

fit_res = residuals(GCP_TABLE, transform)
chk_res = residuals(CHECKPOINTS, transform)

fit_rmse = float(np.sqrt(np.mean(fit_res ** 2)))
print("Per-point fit residuals (m):", np.round(fit_res, 3).tolist())
print(f"Fit RMSE:        {fit_rmse:.3f} m")
if len(chk_res):
    print(f"Checkpoint RMSE: {float(np.sqrt(np.mean(chk_res ** 2))):.3f} m")

# ── 5. Verification assert on GCP residuals ────────────────────────────────
assert fit_rmse < RMSE_BUDGET_M, (
    f"GCP RMSE {fit_rmse:.3f} m exceeds budget {RMSE_BUDGET_M:.3f} m — "
    "check for a mislabelled point, wrong CRS, or over-fit transform."
)
print(f"OK: georeferenced {DST_PATH} within {RMSE_BUDGET_M:.3f} m budget")

For a thin-plate spline or a higher-order polynomial, swap Step 3 for gdal.Warp, which handles those non-affine transforms natively:

PYTHON 
from osgeo import gdal

gcp_list = [gdal.GCP(e, n, 0, c, r) for (c, r, e, n) in GCP_TABLE]  # x,y,z,pixel,line
src_ds = gdal.Open("ortho_raw.tif", gdal.GA_Update)
src_ds.SetGCPs(gcp_list, gdal.osr.GetUserInputAsWKT("EPSG:32615"))

gdal.Warp(
    "ortho_georef_tps.tif",
    src_ds,
    tps=True,                      # thin-plate spline; or polynomialOrder=2
    dstSRS="EPSG:32615",
    resampleAlg="bilinear",
    xRes=0.03, yRes=0.03,
    targetAlignedPixels=True,
)
src_ds = None

Inline verification: confirm the output carries the intended CRS and covers a sane spatial extent:

PYTHON 
with rasterio.open("ortho_georef.tif") as chk:
    print(chk.profile)
    assert chk.crs.to_epsg() == 32615, "Output CRS is not the target UTM zone"
    left, bottom, right, top = chk.bounds
    assert right > left and top > bottom, "Degenerate extent — transform failed"
    print(f"Extent: {right - left:.1f} m × {top - bottom:.1f} m")

Gotchas and Edge Cases

  • Mixed CRS in the GCP list. If some world coordinates are WGS84 lat/lon (small decimal degrees) and others are UTM (six/seven-digit metres), from_gcps still returns a transform and the fit looks plausible — but every residual is nonsense and the raster lands off the planet. Project all points into TARGET_CRS before building the list, and sanity-check the coordinate magnitudes.
Why a low residual is not the same as an accurate mosaic Two panels on ground control point placement. Points clustered in one corner constrain the solve locally and produce excellent residuals while the far end of the block is extrapolated; points spread across the block with one held back give a residual that is an honest accuracy estimate. Points clustered in one corner All five markers near the access track The solve is well constrained there. The far end of the block is extrapolated. Residuals look excellent and mean nothing. Error grows with distance from the cluster. Points spread, one held back Four near the corners, one in the middle The solve is constrained across the block. The held-out point is never fitted. Its residual is an honest accuracy estimate. This is the number worth recording.

  • Thin-plate spline over-fitting. TPS honours every point exactly, so a single mislabelled marker (rows/cols swapped, or clicked on the wrong target) drags the local surface into a visible warp near that point while the aggregate RMSE stays deceptively low. Always read the per-point residual list, not just the mean, and prefer an affine fit unless doming genuinely demands local warping.

  • Row/col versus x/y confusion. GDAL’s GCP(x, y, z, pixel, line) takes world coordinates first then pixel column (pixel) and row (line); rasterio’s GroundControlPoint(row=, col=, x=, y=) names them. Passing pixel column where a row is expected mirrors or transposes the whole mosaic. The residual assert catches it — the RMSE explodes — but name the arguments explicitly to avoid the trap.

  • Trusting RTK image tags blindly. RTK-tagged imagery is only as good as the base station’s known position. If the base sat on an autonomous (unsurveyed) fix, every image shares a common absolute offset of up to a metre or two even though relative accuracy is centimetric. One surveyed checkpoint exposes that offset; without it the error is invisible.

Frequently Asked Questions

How many ground control points do I need for a field survey?

A first-order polynomial needs at least three points, but for a production field survey use five to ten well-distributed points including the corners and the centre. A thin-plate spline needs more, roughly one every few hectares, because it warps locally between points. Always keep one or two points out of the fit as independent checkpoints so your reported accuracy is not self-referential.

Do I still need ground control points if my drone has RTK GPS?

For most relative measurements such as plant counts, canopy area, and vegetation index maps, RTK-tagged images give one to three centimetre relative accuracy and ground control points are unnecessary. You still want at least one or two checkpoints to verify absolute accuracy and detect a vertical datum offset. Add full ground control when you must tie the map to a legal boundary, stack multi-date flights to the millimetre, or when the RTK base was on an unknown or autonomous position.

What RMSE should I expect from a good GCP georeferencing run?

With RTK-measured points at one to three centimetre accuracy and a clean fit, horizontal RMSE should land within one to two ground sample distances, typically a few centimetres for a low-altitude flight. If residuals blow past a tenth of a metre, suspect a mislabelled point, a wrong CRS, or an over-fit high-order polynomial reacting to a single bad marker. Inspect per-point residuals rather than trusting the aggregate number.

Parent Guide

This guide is part of Orthomosaic Stitching Workflows — see there for the full stitching pipeline, from image ingestion and block adjustment through to the georeferenced mosaic this page validates.

Parent section Orthomosaic Stitching for Agricultural Drone Imagery