Troubleshooting GCP and Coordinate Errors

Ground control point failures rarely announce themselves as errors. More often the pipeline runs to completion, the orthomosaic looks plausible, and only a check against surveyed truth reveals that the whole model is offset by a UTM zone, tilted by a vertical datum confusion, or warped because one mis-tagged marker dragged the least-squares solution toward it. This page is the aggregator for that failure class: it gives a repeatable diagnostic methodology that turns a vague “the georeferencing is wrong” into a specific, code-checkable root cause — axis-order swap, wrong UTM zone or flipped hemisphere, ellipsoidal-versus-orthometric height mismatch, GCP mis-tagging or bad image observations, high residuals from an over-weighted marker, a malformed gcp_list.txt, or too few and poorly distributed control points. It belongs to the ground control point optimization and coordinate sync stage, and it deliberately routes each confirmed cause to the deep-dive that fixes it rather than duplicating those fixes here.

Audience and prerequisites. This guide targets Python 3.10+ on a 64-bit OS with at least 8 GB RAM. You should be comfortable reading a residual report, know your project’s target CRS and vertical datum, and understand the difference between a horizontal error (metres of offset in easting/northing) and a vertical one (metres of offset in height). The methodology is ordered deliberately: cheap, deterministic checks first (axis order, zone) before expensive ones (re-tagging markers), because the earlier a failure is located the less compute it wastes. The detection side of this work is covered in automating GCP detection with Python, and the transform side in coordinate transformation workflows in PyProj.

# Prerequisites

The diagnostic snippets below read residual reports, re-validate CRS declarations, and compare coordinate bands. They need only the core geospatial stack; no reconstruction engine is invoked during diagnosis.

Keep one pinned PROJ across the environment; a diagnosis that assumes the field laptop and the processing server resolve datum shifts identically is only valid if their PROJ grids match.

# Conceptual architecture: the diagnostic decision tree

Effective GCP debugging is a sequence of elimination gates, ordered cheapest-first. Each gate asks one yes/no question about the coordinate data; a “no” routes you to a specific, named cause and its fix guide, while a “yes” advances to the next, more expensive question. Only when every gate passes should you conclude the residuals are irreducible and the control network itself needs densifying. The tree below is the exact order the numbered steps implement.

The numbered steps below walk the spine top to bottom. Run them in order; skipping ahead to re-tagging markers before you have ruled out an axis swap is the classic way to waste an afternoon “fixing” data that was never wrong.

# Step 1: Reproduce the failure and build a symptom signature

Before touching any fix, quantify what “wrong” means. Read the residual report or compare a handful of transformed control coordinates against their surveyed truth, and record the shape of the error: is it a constant offset (whole model shifted), a systematic tilt (grows across the site), a single spike (one marker), or a vertical-only discrepancy? These four shapes map almost one-to-one onto the four gates. The helper below turns a control-versus-computed table into a signature.

import numpy as np

def error_signature(surveyed: np.ndarray, computed: np.ndarray) -> dict:
    """Classify the shape of a control-vs-computed coordinate error.

    surveyed / computed: (N, 3) arrays of easting, northing, height (metres).
    """
    d = computed - surveyed                       # per-point residual vectors
    horiz = np.hypot(d[:, 0], d[:, 1])            # horizontal error magnitude
    vert = np.abs(d[:, 2])                        # vertical error magnitude
    mean_h = float(horiz.mean())
    return {
        "mean_horizontal_m": round(mean_h, 3),
        "mean_vertical_m": round(float(vert.mean()), 3),
        "max_single_horiz_m": round(float(horiz.max()), 3),
        # a near-constant offset points at zone/axis; a lone spike at tagging.
        "looks_like_constant_offset": bool(horiz.std() < 0.2 * mean_h + 1e-9),
        "vertical_dominates": bool(vert.mean() > 3 * horiz.mean() + 1e-9),
    }

A looks_like_constant_offset of True with a large magnitude sends you straight to Steps 2 and 3; vertical_dominates jumps you to Step 4’s height gate; a single dominant max_single_horiz_m points at a mis-tagged marker in Step 5. The residual report from OpenDroneMap is the richest source for this signature — reading it marker by marker is covered in interpreting GCP residual reports in ODM.

# Step 2: Rule out axis-order and CRS misdeclaration first

The cheapest and most common failure is an axis swap: coordinates passed as (lat, lon) into a transformer expecting (lon, lat), or a CRS object honouring an authority’s latitude-first order. The signature is a large, structured offset where easting and northing look transposed. Re-validate the declared CRS and confirm the transform direction before anything else.

import pyproj

def check_axis_and_crs(sample_lon: float, sample_lat: float, target_epsg: str) -> None:
    tgt = pyproj.CRS.from_user_input(target_epsg)
    assert tgt.is_projected, f"{target_epsg} is not projected — cannot hold metres"

    # always_xy=True forces (lon, lat) in; without it EPSG:4326 is lat-first.
    tf = pyproj.Transformer.from_crs("EPSG:4326", target_epsg, always_xy=True)
    e, n = tf.transform(sample_lon, sample_lat)

    # A swapped call would put a longitude-scaled value where northing belongs.
    assert 100_000 <= e <= 900_000, f"easting {e:.1f} outside UTM band — check axis order"
    print(f"axis/CRS OK: E={e:.1f} N={n:.1f} via {tgt.name}")

If this gate fails, the fix is not in this page — it is the always_xy / axis-mapping discipline detailed in coordinate transformation workflows in PyProj. Rebuild every transformer with the correct axis contract and re-derive the control coordinates before proceeding.

# Step 3: Confirm the UTM zone and hemisphere

A model that is internally consistent but globally displaced by hundreds of kilometres is almost always a wrong UTM zone or a flipped hemisphere — a southern-hemisphere site processed as northern, so northings are off by 10,000,000 m (the false-northing offset). The check derives the correct zone from a representative longitude/latitude and compares it to the declared EPSG.

import math

def expected_utm_epsg(lon: float, lat: float) -> str:
    """Return the EPSG code for the UTM zone that actually contains this point."""
    zone = int(math.floor((lon + 180) / 6) % 60) + 1
    # 326xx = northern hemisphere, 327xx = southern.
    base = 32600 if lat >= 0 else 32700
    return f"EPSG:{base + zone}"

def check_zone(lon: float, lat: float, declared_epsg: str) -> None:
    want = expected_utm_epsg(lon, lat)
    assert declared_epsg.upper() == want, (
        f"declared {declared_epsg} but point belongs in {want} "
        f"(wrong zone or hemisphere flip)")
    print(f"zone OK: {want}")

A hemisphere flip is the sub-case that catches people out: the horizontal position can look right in easting but be exactly 10,000 km wrong in northing. If the declared code and the derived code disagree, reproject with the correct one — the transform mechanics are again in the coordinate transformation guide.

# Step 4: Separate vertical-datum errors from horizontal ones

When vertical_dominates is true, heights are off by a roughly constant tens-of-metres bias while horizontal position is fine. This is the ellipsoidal-versus-orthometric mismatch: GNSS reports ellipsoidal height, the survey control is orthometric (geoid-referenced), and no geoid model bridged them. The bias equals the local geoid undulation. The check flags a suspiciously constant vertical offset.

import numpy as np

def check_vertical_datum(surveyed_z: np.ndarray, computed_z: np.ndarray) -> None:
    dz = computed_z - surveyed_z
    bias, spread = float(dz.mean()), float(dz.std())
    # A near-constant multi-metre bias with tight spread == geoid undulation,
    # not random noise. Typical undulations run tens of metres.
    if abs(bias) > 2.0 and spread < 1.0:
        raise AssertionError(
            f"constant vertical bias {bias:.2f} m — ellipsoidal/orthometric "
            f"mismatch; apply a geoid model via a compound CRS")
    print(f"vertical datum OK: bias {bias:.2f} m, spread {spread:.2f} m")

The remedy is a compound or 3D target CRS that pairs the horizontal grid with a geoid-based height, so PROJ applies the correct undulation. Confirm the geoid grid is installed before trusting the fix, exactly as the vertical-datum troubleshooting in the PyProj transformation guide describes.

# Step 5: Inspect GCP tagging, observations, distribution, and file format

If the first four gates pass and residuals are still high, the problem is in the control network itself: a mis-clicked image observation, a marker with too few images, a malformed gcp_list.txt, or control clustered in one corner leaving the rest of the block unconstrained. Parse the GCP file and check both its format and its geometry before re-running.

from pathlib import Path

def audit_gcp_list(path: Path, min_obs: int = 3) -> list[str]:
    """Validate ODM-style gcp_list.txt structure and per-marker observation count."""
    problems: list[str] = []
    lines = [ln.strip() for ln in Path(path).read_text().splitlines() if ln.strip()]
    if not lines:
        return ["empty gcp_list.txt"]
    if not (lines[0].upper().startswith("EPSG:") or "+proj" in lines[0].lower()):
        problems.append("first line must be a CRS (EPSG:xxxxx or a PROJ string)")

    obs: dict[str, int] = {}                       # image-observations per marker
    for ln in lines[1:]:
        cols = ln.split()
        # geo_x geo_y geo_z im_x im_y image_name [gcp_name]
        if len(cols) < 6:
            problems.append(f"row has {len(cols)} cols, need >= 6: {ln[:40]}")
            continue
        name = cols[6] if len(cols) > 6 else f"{cols[0]}_{cols[1]}"
        obs[name] = obs.get(name, 0) + 1

    for name, count in obs.items():
        if count < min_obs:
            problems.append(f"marker {name} has {count} observations (< {min_obs})")
    return problems

Under-observed or clustered control is a distribution problem: how those residuals then get spread across the model is the subject of distributing GCP errors across orthomosaics, and the acceptance thresholds that decide whether the remaining residuals are tolerable at all come from setting accuracy thresholds for survey projects.

# Symptom-signature reference

Read the observed symptom in the left column, match it to the signature, and jump to the step that isolates it.

# Verification and output inspection

After applying a fix, re-run the signature helper on the corrected control set and assert the errors dropped below the project tolerance. The verification is the same regardless of which gate you fixed — a clean signature is the only proof the fix worked.

import numpy as np

def verify_fixed(surveyed: np.ndarray, computed: np.ndarray,
                 h_tol: float = 0.05, v_tol: float = 0.10) -> None:
    d = computed - surveyed
    horiz = np.hypot(d[:, 0], d[:, 1])
    vert = np.abs(d[:, 2])
    assert horiz.max() <= h_tol, f"horizontal residual {horiz.max():.3f} m > {h_tol} m"
    assert vert.max() <= v_tol, f"vertical residual {vert.max():.3f} m > {v_tol} m"
    # No single marker should carry an outsized share of the total error.
    assert horiz.max() < 3 * horiz.mean() + 1e-9, "one marker still dominates — re-check tagging"
    print(f"verified: max horiz {horiz.max():.3f} m, max vert {vert.max():.3f} m")

If the per-marker dominance assertion still fires after re-tagging, the offending observation is genuinely bad and should be dropped rather than re-weighted; the residual report tells you which one in interpreting GCP residual reports in ODM.

# Troubleshooting

My orthomosaic is offset by a constant amount everywhere but the shape is correct. A uniform offset across the whole model is a datum or zone problem, not a tagging one. Run Step 2 and Step 3: confirm always_xy is set and that the declared UTM zone matches the one derived from a representative coordinate. A hemisphere flip shows up as a northing that is wrong by the 10,000,000 m false-northing offset.

One control point shows a huge residual and the rest are fine. That marker’s image observations are almost certainly mis-clicked, or its surveyed coordinate was transcribed wrong. Do not re-weight it — inspect its residual in the ODM report, correct the observation if the image pick was off, and if the surveyed value itself is suspect, drop the marker and re-run. Step 5’s audit flags markers with too few observations that make a bad pick impossible to detect.

Horizontal accuracy is fine but every height is off by the same ~30 metres. This is the ellipsoidal-versus-orthometric mismatch from Step 4. GNSS gives ellipsoidal height, your control is orthometric, and the constant bias is the local geoid undulation. Fix it with a compound/3D target CRS that applies a geoid model, and confirm the geoid grid is installed rather than silently skipped.

ODM rejects my gcp_list.txt before processing starts. The file format is strict: the first line must be a CRS (an EPSG:xxxxx code or a PROJ string), and each subsequent row needs at least six whitespace-separated columns — geo x/y/z, image x/y, and the image filename. Run the audit_gcp_list helper from Step 5; a common cause is a missing CRS header line or rows padded with tabs that split into the wrong column count.

Residuals are low on average but corner tiles are visibly warped. Low mean RMSE hides clustered residual. Your control is too sparse or bunched in one region, leaving the extremities unconstrained. Densify control at the corners and along the perimeter, then spread the remaining error with the weighting approach in distributing GCP errors across orthomosaics; a good average is not the same as a well-conditioned network.

Two machines produce different coordinates for the same GCP set. The installs resolve different PROJ grids or versions, so datum shifts differ. Pin pyproj and GDAL to one shared PROJ, point both at the same PROJ_DATA directory, and re-run the diagnosis; until parity is confirmed, any cross-machine residual comparison is meaningless.

# Related

• Automating GCP Detection with Python
• Coordinate Transformation Workflows in PyProj
• Distributing GCP Errors Across Orthomosaics
• Setting Accuracy Thresholds for Survey Projects
• Interpreting GCP Residual Reports in ODM

← Ground Control Point Optimization & Coordinate Sync