RTK and PPK Geotagging Workflows in Python

An ordinary drone writes a GPS position into each image’s EXIF that is good to a few metres. An RTK or PPK-equipped aircraft records a trajectory good to a couple of centimetres — but that trajectory is a separate file, sampled on its own clock, referenced to the antenna rather than to the camera. Turning it into a per-image camera station is a short, unglamorous pipeline that every survey-grade workflow depends on, and every step of it has a way to be quietly wrong.

This page builds that pipeline: parse the trajectory, pair each exposure with the trajectory instant it actually occurred at, correct for the offset between antenna and sensor, and write positions with honest uncertainties. Done correctly it reduces the ground control needed for a given accuracy class; done carelessly it produces a survey that is precisely wrong, which is worse than one that is imprecisely right.

The accuracy this feeds is measured against the thresholds in setting accuracy thresholds for survey projects, and the positions it produces are consumed by the same solver that uses ground control points.

Audience and prerequisites. Python 3.10+, a post-processed trajectory from the flight, and the camera’s EXIF. Familiarity with coordinate reference systems as handled in coordinate transformation workflows in PyProj is assumed.

# Prerequisites

# Conceptual architecture

Four quantities have to be reconciled: where the antenna was, when the shutter fired, where the sensor is relative to the antenna, and in what frame the answer should be expressed. Each is a separate source of error and each has a distinct fix.

Figure 1 — The chain from trajectory to camera station. The magnitudes matter: skipping interpolation costs more than the lever arm, and getting the frame wrong costs more than both together.

# Step 1: Pair each exposure with a trajectory instant

The trajectory is sampled on a fixed rate; the shutter fires whenever it fires. Rounding to the nearest trajectory sample introduces a position error equal to the ground speed times half the sample interval — at 12 m/s and 5 Hz, that is 1.2 m, which discards the entire benefit of the RTK receiver.

Linear interpolation between the bracketing samples reduces that to the trajectory’s own noise, because over a tenth of a second an aircraft’s motion is very nearly linear.

import numpy as np
import pandas as pd


def interpolate_at_events(traj: pd.DataFrame, events: pd.Series) -> pd.DataFrame:
    """Antenna position at each exposure instant.

    traj: columns gps_time, east, north, up (already in the working CRS)
    events: GPS times of the shutter, one per image, index = image name
    """
    t = traj["gps_time"].to_numpy(dtype=float)
    if not np.all(np.diff(t) > 0):
        raise ValueError("trajectory times are not strictly increasing")

    te = events.to_numpy(dtype=float)
    if te.min() < t.min() or te.max() > t.max():
        raise ValueError("an exposure falls outside the trajectory span")

    out = {axis: np.interp(te, t, traj[axis].to_numpy(dtype=float))
           for axis in ("east", "north", "up")}
    return pd.DataFrame(out, index=events.index)

The two guards are not decoration. A trajectory with a non-monotonic time column — which happens when two processing sessions are concatenated — makes np.interp return silently wrong values rather than raising. And an exposure outside the trajectory span means the logger started late or stopped early, and extrapolation there is unbounded.

# Step 2: Apply the lever arm

The antenna is on top of the aircraft; the sensor is underneath it. On a typical multirotor that offset is 20–40 cm, and it rotates with the aircraft: at a heading of 0° the antenna may be 25 cm above and 5 cm behind the sensor, and at 180° the horizontal component points the other way. Ignoring it leaves a heading-dependent error that partly cancels over a serpentine flight and does not cancel at all on a single-direction corridor.

import numpy as np


def apply_lever_arm(east, north, up, yaw_deg, pitch_deg, roll_deg,
                    arm_body=(0.05, 0.0, -0.28)):
    """Shift antenna positions to the sensor.

    arm_body is the sensor relative to the antenna in the body frame,
    metres, as (forward, right, down). Rotation is yaw-pitch-roll applied
    in that order, which matches the convention of most flight logs.
    """
    y, p, r = (np.radians(np.asarray(a, dtype=float))
               for a in (yaw_deg, pitch_deg, roll_deg))
    fx, fy, fz = arm_body

    # Body → local level (east, north, up).
    de = (fx * (np.cos(y) * np.cos(p))
          + fy * (np.cos(y) * np.sin(p) * np.sin(r) - np.sin(y) * np.cos(r))
          + fz * (np.cos(y) * np.sin(p) * np.cos(r) + np.sin(y) * np.sin(r)))
    dn = (fx * (np.sin(y) * np.cos(p))
          + fy * (np.sin(y) * np.sin(p) * np.sin(r) + np.cos(y) * np.cos(r))
          + fz * (np.sin(y) * np.sin(p) * np.cos(r) - np.cos(y) * np.sin(r)))
    du = (-fx * np.sin(p) + fy * np.cos(p) * np.sin(r) + fz * np.cos(p) * np.cos(r))
    return east + de, north + dn, up + du

Getting the sign of the vertical component right is the most common error, and it is worth verifying empirically rather than from a diagram: compare the resulting camera heights against a known ground elevation plus the flight height, and confirm the correction moved them the way you expected.

# Step 3: Declare the frame, including its epoch

An RTK correction service delivers positions in the frame of its reference stations, which for a global service is a realisation of ITRF at the epoch of observation. A national grid is usually a plate-fixed frame at a fixed epoch. Between the two lies plate motion of a few centimetres per year — irrelevant at metre accuracy and dominant at centimetre accuracy.

Figure 2 — Why the epoch belongs in the manifest. A correctly RTK-positioned survey delivered in the wrong epoch is wrong by more than the receiver’s precision, and the error grows every year.

The practical rule is to keep two things explicit: the frame the trajectory was computed in, and the frame the deliverable is contracted in. Where they differ, the transformation between them is time-dependent and pyproj will apply it if — and only if — both epochs are supplied.

# Step 4: Write positions with honest uncertainties

The solver weights camera positions by their stated accuracy. A pipeline that writes centimetre positions with a default metre uncertainty gets no benefit from the receiver, and one that writes centimetre uncertainties for a fix that was actually a float solution over-trusts them.

Post-processing software reports a solution quality per epoch. Carry it through: fixed-integer epochs get the receiver’s specification, float epochs get several times that, and single-point epochs are excluded from the weighting entirely rather than being pretended into the block.

QUALITY_SIGMA_M = {1: 0.02, 2: 0.15, 5: 2.0}   # fixed, float, single


def station_sigma(quality: int, base_h: float = 0.02, base_v: float = 0.03):
    """Horizontal and vertical one-sigma for a camera station."""
    scale = QUALITY_SIGMA_M.get(int(quality), 5.0) / QUALITY_SIGMA_M[1]
    return base_h * scale, base_v * scale

Vertical uncertainty is deliberately larger than horizontal: GNSS geometry constrains height about half as well as plan position, and stating them as equal is a claim the receiver does not make.

# How much ground control an RTK survey still needs

The usual reason to fit an RTK receiver is to reduce the ground control burden, and the usual overcorrection is to eliminate it. Accurate camera positions constrain the translation of a block extremely well and constrain two other things poorly.

The first is scale in the vertical, which is weakly determined by nadir imagery regardless of how well the camera positions are known, because the rays that intersect at a ground point are close to parallel. A small error in the assumed focal length or in the lever arm’s vertical component appears as a scale error in height, and camera positions alone cannot separate it from a genuine elevation difference.

The second is systematic deformation — the gentle bowl that appears when self-calibration absorbs radial distortion imperfectly. Camera positions constrain where the block sits, not whether it is flat, so a bowl of several centimetres can coexist with camera stations that all agree with the trajectory to a centimetre.

Both are exactly what a small number of well-placed check points detect, which is why the honest reduction is fewer control points rather than none. A workable pattern for a well-configured RTK block is three to five surveyed points: one near each end of the longest dimension and one near the centre, used as check points rather than as control. If they agree with the reconstruction, the block is validated and no control was needed. If they do not, they are already in place to be promoted to control and the flight does not need repeating.

That arrangement also protects against the failure this whole workflow is most exposed to: a systematic error in the trajectory itself. A base station entered at the wrong height, or a correction stream referenced to a different frame, biases every camera position identically — and a block constrained only by camera positions will happily reproduce that bias with excellent internal residuals. Independent ground truth is the only thing that sees it.

# Parameter deep-dive

# Verification and output inspection

Two checks catch nearly everything. The first compares the corrected stations against the raw EXIF positions: the difference should be small, structured, and consistent with the lever arm — not random and not metres.

import numpy as np


def sanity_check_stations(corrected, exif, ground_speed_ms: float = 12.0) -> None:
    """The correction should be centimetres, not metres, and not random."""
    d = np.linalg.norm(corrected[["east", "north"]].to_numpy()
                       - exif[["east", "north"]].to_numpy(), axis=1)
    median, p99 = float(np.median(d)), float(np.percentile(d, 99))
    assert median < 5.0, f"median shift {median:.2f} m — check the frame, not the arm"
    assert p99 < 15.0, f"99th percentile shift {p99:.2f} m — outliers in the pairing"

    # A residual clock offset shows up as a shift along the flight direction.
    step = np.diff(corrected[["east", "north"]].to_numpy(), axis=0)
    heading = step / (np.linalg.norm(step, axis=1, keepdims=True) + 1e-9)
    along = np.einsum("ij,ij->i", heading, (corrected[["east", "north"]]
                                            .to_numpy()[:-1]
                                            - exif[["east", "north"]].to_numpy()[:-1]))
    bias = float(np.median(along))
    assert abs(bias) < 0.5, (
        f"along-track bias {bias:.2f} m ≈ {bias / ground_speed_ms * 1000:.0f} ms "
        "of clock offset — reconcile the event markers")

The along-track projection is the useful part: a clock offset produces a shift in the direction of travel and nothing else, so projecting the correction onto the flight heading isolates it from every other error source. A 40 ms residual at 12 m/s is half a metre, and it is invisible in a plan view of the positions.

The second check is the only one that measures accuracy rather than consistency: withhold surveyed check points and compare, exactly as in control points versus checkpoints.

# Troubleshooting

Every camera station is shifted by a constant metre-scale amount. A frame or epoch mismatch, not a lever arm — the lever arm cannot produce a metre. Compare the declared frame of the trajectory against the delivery frame and check whether both epochs were supplied to the transformation.

The shift is constant in magnitude but rotates with the flight lines. That is the lever arm, uncorrected or corrected with the wrong sign. Its horizontal component points along the aircraft heading, so it reverses between adjacent strips of a serpentine pattern.

Interpolation raises “an exposure falls outside the trajectory span”. The receiver logged a shorter interval than the camera shot for. The frames outside the span have no valid position; drop them rather than extrapolating, and check whether the logger was started after the first exposure.

The reconstruction is worse with PPK positions than with EXIF. Almost always over-tight uncertainties: stations stated at 2 cm when the solution was float force the solver to fit noise. Check the fixed-solution fraction and scale the sigmas by quality.

Positions are excellent in plan and poor in height. Either the vertical uncertainty was stated equal to the horizontal, or the antenna height above the sensor was applied with the wrong sign, which produces exactly twice the lever arm as a constant vertical bias.

Event markers outnumber images, or the reverse. The camera triggered without recording, or recorded without triggering the marker. Match on time rather than on index, and reject the batch if the count difference exceeds a frame or two — an index-based pairing after a dropped frame misaligns every subsequent image.

# Related

• Coordinate transformation workflows in PyProj
• Setting accuracy thresholds for survey projects
• Applying PPK corrections to image timestamps
• Fixing camera event marker offsets in PPK logs
• How to validate EXIF GPS data before processing

← Ground Control Point Optimization & Coordinate Sync

Figure 3 — The diagnostic that separates the two constant-magnitude errors. A serpentine flight distinguishes them for free; a corridor does not, and that is worth knowing before the flight rather than after.