Preparing Ground Control Point Files

This page builds and validates the ground control file a photogrammetry run consumes — the CRS header, the image observations, the distribution that actually constrains a solution, the check points held back from it, and the validation that catches a swapped easting, a typo in a height or a point observed in only one image, all against EPSG:25832+7837.

Why you hit this

Control points are the only mechanism that removes systematic deformation from a reconstruction, and the file that carries them is plain text with no schema. A single transposed digit moves one point by a hundred metres, the solver either rejects it or bends the whole model towards it, and the result looks like a processing problem rather than a data-entry one. Worse, the failure often appears as a small, plausible tilt rather than an obvious break. Validating the file before a six-hour run costs a minute and prevents re-running. The role control plays in the pipeline is set out in photogrammetry processing pipelines.

Prerequisites

  • Surveyed points in the target compound CRS — EPSG:25832+7837 here — with their survey accuracy recorded.
  • The flight’s imagery, with EXIF positions, so observations can be checked against where the camera was.
  • Python 3.10+ with numpy>=1.24, pyproj>=3.6, pillow>=10, exifread>=3.0.
  • A marking tool for picking image coordinates: the ODM or WebODM GCP interface, or any viewer that reports pixel coordinates with the origin at the top left.

Step-by-Step

1. Understand the format

text
EPSG:25832+7837
691204.412 5335818.221 519.310 2104 1436 DJI_0123.JPG gcp_A
691204.412 5335818.221 519.310 1876 2210 DJI_0124.JPG gcp_A
691204.412 5335818.221 519.310 3012 1104 DJI_0125.JPG gcp_A
691286.901 5335836.978 519.290 980 1502 DJI_0141.JPG gcp_B

The first line is the coordinate reference system, as an EPSG code or a PROJ string. Every following line is one observation: the point’s three world coordinates, then its pixel coordinates in one image, then the image name, then an optional label. A point observed in four images appears on four lines with identical world coordinates and different pixel coordinates.

Two conventions cause most of the confusion. World coordinates are easting, northing, height in the declared CRS — not longitude, latitude — and pixel coordinates count from the top-left corner of the image, x to the right and y downward. Some survey exports give northing first, and some marking tools report pixels from the bottom left; both produce a file that parses perfectly and reconstructs wrongly.

2. Write the file from a survey table and a marking table

python
from pathlib import Path
import numpy as np

CRS = "EPSG:25832+7837"

survey = {                      # label → (easting, northing, orthometric height)
    "gcp_A": (691204.412, 5335818.221, 519.310),
    "gcp_B": (691286.901, 5335836.978, 519.290),
    "gcp_C": (691275.305, 5335887.874, 519.330),
    "gcp_D": (691192.799, 5335869.108, 519.300),
    "gcp_E": (691240.118, 5335852.664, 521.045),
}
observations = [                # (label, image, px, py) from the marking session
    ("gcp_A", "DJI_0123.JPG", 2104, 1436),
    ("gcp_A", "DJI_0124.JPG", 1876, 2210),
    ("gcp_A", "DJI_0125.JPG", 3012, 1104),
    ("gcp_B", "DJI_0141.JPG", 980, 1502),
    # …
]

def write_gcp_list(path, crs, survey, observations):
    lines = [crs]
    for label, image, px, py in observations:
        e, n, h = survey[label]
        lines.append(f"{e:.3f} {n:.3f} {h:.3f} {px} {py} {image} {label}")
    Path(path).write_text("\n".join(lines) + "\n")
    return len(lines) - 1

n_obs = write_gcp_list("gcp_list.txt", CRS, survey, observations)
print(f"{n_obs} observations of {len({o[0] for o in observations})} points written")

Generating the file from two tables rather than editing it by hand is what keeps the world coordinates consistent across a point’s observations. A hand-edited file where one of four lines for gcp_A has a different height is accepted by the solver and quietly inconsistent.

3. Validate before running anything

python
import exifread
from pyproj import CRS as PCRS, Transformer

def validate(path, image_dir, expect_epsg=25832, max_camera_dist_m=400.0):
    lines = Path(path).read_text().strip().splitlines()
    header, rows = lines[0], [l.split() for l in lines[1:] if l.strip()]
    problems, counts = [], {}

    crs = PCRS.from_user_input(header.split("+")[0] if header.startswith("EPSG") else header)
    if crs.to_epsg() != expect_epsg:
        problems.append(f"header CRS {header} is not the expected EPSG:{expect_epsg}")
    if header.startswith("EPSG") and "+" not in header:
        problems.append("header has no vertical CRS: heights will be treated as ellipsoidal")

    to_wgs = Transformer.from_crs(header, "EPSG:4326", always_xy=True)
    coords_by_label = {}

    for i, r in enumerate(rows, start=2):
        if len(r) < 6:
            problems.append(f"line {i}: expected at least 6 fields, got {len(r)}")
            continue
        e, n, h, px, py, image = float(r[0]), float(r[1]), float(r[2]), float(r[3]), float(r[4]), r[5]
        label = r[6] if len(r) > 6 else f"line{i}"
        counts[label] = counts.get(label, 0) + 1
        coords_by_label.setdefault(label, set()).add((e, n, h))

        if not (200_000 < e < 900_000) or not (4_000_000 < n < 6_500_000):
            problems.append(f"line {i}: {e}, {n} is not a plausible UTM easting/northing — swapped?")
        if not (-100 < h < 3000):
            problems.append(f"line {i}: height {h} is implausible")

        img_path = Path(image_dir) / image
        if not img_path.exists():
            problems.append(f"line {i}: image {image} not found")
            continue
        with open(img_path, "rb") as f:
            tags = exifread.process_file(f, details=False)
        w = int(str(tags.get("EXIF ExifImageWidth", 0)) or 0)
        hgt = int(str(tags.get("EXIF ExifImageLength", 0)) or 0)
        if w and hgt and not (0 <= px <= w and 0 <= py <= hgt):
            problems.append(f"line {i}: pixel ({px}, {py}) is outside the {w}×{hgt} image")

    for label, coords in coords_by_label.items():
        if len(coords) > 1:
            problems.append(f"{label}: inconsistent world coordinates across observations")
    for label, n_obs in counts.items():
        if n_obs < 3:
            problems.append(f"{label}: only {n_obs} observation(s); three or more are needed")

    return problems, counts

problems, counts = validate("gcp_list.txt", "/data/odm/block_09/images")
print(f"{len(counts)} points, observations per point: {counts}")
for p in problems:
    print("PROBLEM:", p)
assert not problems, "fix the control file before running"

Each check maps to a real failure. The plausibility bounds catch swapped easting and northing, which is otherwise undetectable in a file of large numbers. The pixel-bounds check catches a marking tool whose coordinate origin differs from the expected one, because a y flipped about the image centre usually still lands inside the image but a systematically offset one does not. The consistency check catches hand edits, and the observation count catches points that will contribute nothing.

What each validation check catches A table of checks. The CRS header check catches a missing vertical datum. Coordinate plausibility catches swapped easting and northing. Height plausibility catches typos. Pixel bounds catch a flipped image origin. Coordinate consistency catches hand edits. The observation count catches points that constrain nothing. check failure it catches compound CRS in the headerheights silently treated as ellipsoidal easting/northing boundscoordinates written in the wrong order height boundsa transposed digit in a survey height pixel inside the imagea marking tool with a flipped y origin one coordinate per labela hand-edited line that drifted
None of these are exotic; all of them produce a file that parses and a reconstruction that is subtly wrong.
Anatomy of one observation line A single line of a control file broken into its seven fields: easting, northing and height in the declared coordinate reference system, then the pixel x and pixel y in the image with the origin at the top left, then the image file name and an optional label. The header line above declares the compound coordinate reference system. EPSG:25832+7837 691204.412 5335818.221 519.310 2104 1436 DJI_0123.JPG gcp_A easting, northing, height in the declared CRS pixels, origin top-left image, then label The header is the only place the vertical datum appears — omit it and heights become ellipsoidal.
Seven fields, two coordinate conventions and one header: nearly every control-file defect is a confusion between two of them.

4. Check the distribution, not just the count

python
from shapely.geometry import MultiPoint, Point
import numpy as np

pts = np.array([survey[l] for l in survey])
hull = MultiPoint([Point(p[0], p[1]) for p in pts]).convex_hull
site = MultiPoint([Point(x, y) for x, y in flight_footprint_coords]).convex_hull   # from image positions

coverage = hull.area / site.area
centroid_offset = Point(pts[:, 0].mean(), pts[:, 1].mean()).distance(site.centroid)
z_spread = pts[:, 2].max() - pts[:, 2].min()

print(f"control hull covers {coverage * 100:.0f}% of the site, "
      f"centroid offset {centroid_offset:.0f} m, height spread {z_spread:.1f} m")
assert coverage > 0.6, "control points do not span the site: expect a tilt"

Five points in a cluster constrain position and nothing else. What removes deformation is spread: points near the corners of the flight footprint, one near the middle, and — where the site has relief — points at different heights, because a control set that is flat in z leaves the vertical scale poorly determined. A hull covering more than about 60% of the site area is the practical rule, and the centroid offset catches the case where all the control sits on one side.

Control distribution and what it constrains Three site plans. Five clustered points fix position only, leaving rotation and scale poorly determined. Five points along one edge fix position and one direction, leaving a tilt across the site. Five points at the corners and centre constrain position, rotation, scale and low-order deformation. clustered one edge corners and centre position only tilt across the site position, scale, rotation
The count in a specification is the least important property of a control set; what it spans decides what it can constrain.

5. Split control from check points

python
CHECK_LABELS = {"chk_1", "chk_2", "chk_3"}

control = {k: v for k, v in all_surveyed.items() if k not in CHECK_LABELS}
check = {k: v for k, v in all_surveyed.items() if k in CHECK_LABELS}

write_gcp_list("gcp_list.txt", CRS, control, [o for o in observations if o[0] in control])
Path("check_points.csv").write_text(
    "label,easting,northing,height\n" +
    "\n".join(f"{k},{v[0]:.3f},{v[1]:.3f},{v[2]:.3f}" for k, v in check.items()) + "\n"
)
print(f"{len(control)} control points in the run, {len(check)} check points held back")

Deciding the split in code, from one surveyed table, means the two files can never disagree and nobody can accidentally include a check point in the adjustment. Choose check points that are spread across the site too — three check points in one corner measure the accuracy of that corner.

Expected Output & Verification

text
18 observations of 5 points written
5 points, observations per point: {'gcp_A': 4, 'gcp_B': 4, 'gcp_C': 3, 'gcp_D': 4, 'gcp_E': 3}
control hull covers 71% of the site, centroid offset 14 m, height spread 1.8 m
5 control points in the run, 3 check points held back

After the reconstruction, the file’s quality shows up as residuals. Compare the two sets:

python
import numpy as np

ctrl_resid = np.loadtxt("odm_report/gcp_residuals.csv", delimiter=",", skiprows=1, usecols=(1, 2, 3))
chk_resid = np.loadtxt("check_residuals.csv", delimiter=",", skiprows=1, usecols=(1, 2, 3))
for name, r in (("control", ctrl_resid), ("check", chk_resid)):
    h = np.hypot(r[:, 0], r[:, 1])
    print(f"{name}: horizontal RMS {np.sqrt((h ** 2).mean()) * 100:.1f} cm, "
          f"vertical RMS {np.sqrt((r[:, 2] ** 2).mean()) * 100:.1f} cm")

Control residuals smaller than check residuals is normal and expected — the solver fitted the control. What matters is the ratio: check residuals two or three times the control residuals indicate a solution that is fitting its constraints rather than the site, usually from too few or badly distributed points. Check residuals close to the control residuals, and both close to the survey accuracy, is a healthy reconstruction.

A height spread of only 1.8 m across the control set is worth noting in the report: on a flat site it is unavoidable, and it means the vertical scale rests on the camera calibration rather than on the control.

Performance Notes

  • Marking is the slow part. Three observations per point across five points is fifteen careful picks; budget an hour and do it once, carefully, rather than twice.
  • Reuse control between epochs. Permanent markers surveyed once serve every future flight, which makes repeat surveys both cheaper and directly comparable.
  • Validate in CI. The validation function above runs in under a second and belongs in the same job that launches the processing run, so a bad file fails in seconds rather than after the reconstruction.
  • Keep the marking table, not just the generated file. Re-generating gcp_list.txt from tables lets you change the CRS declaration or drop a point without re-marking.

Common Errors

The solver reports it ignored a control point. Its reprojection residual exceeded the rejection threshold — usually a mismarked observation, or a point whose label is shared with a different physical marker. Check that point’s observations against each other.

The reconstruction is tilted despite good residuals at the control points. Control is clustered or confined to one edge. Residuals can be tiny at the points and large everywhere else.

Heights are systematically offset by tens of metres. The header declared a horizontal CRS only, so orthometric heights were interpreted as ellipsoidal. Always write the compound code, as in step 1.

Pixel coordinates seem to be mirrored vertically. The marking tool used a bottom-left origin. Convert with py = image_height - py and re-validate; the pixel-bounds check will not always catch this, so verify one observation visually.

Frequently Asked Questions

How many control points are enough?

Five well-distributed points plus three check points is a sound minimum for a site of a few hectares. Larger areas need more, roughly one control point per 100–150 m of site extent along each axis, and split-merge runs need points inside each submodel’s overlap.

Can I use existing survey markers instead of placing targets?

Yes, if they are identifiable in the imagery at the flight’s GSD — a painted cross, a manhole centre, a corner of a kerb. Natural features are harder to mark consistently across images, which shows up as larger residuals for those points.

Do check points need to be observed in the images?

They need to be identifiable, so their coordinates can be compared with the reconstruction at the same spot. That is usually by sampling the surface model at the point’s easting and northing, as the verification step in running OpenDroneMap in Docker does.

Back to Photogrammetry Processing Pipelines.