Thermal Orthomosaic Processing in Python

A thermal orthomosaic looks like any other raster and behaves like nothing else in a photogrammetry pipeline. Its pixel values are temperatures, which means they are physical quantities with units, derived from radiance through assumptions that can each be wrong by a degree. Its geometry is poor — an uncooled sensor is typically 640 × 512 pixels with weak texture — so reconstruction from thermal frames alone frequently fails. And its calibration drifts measurably while the aircraft is in the air.

The result is a product where the two halves of the problem pull in opposite directions. The geometric half wants many overlapping frames, which means a long flight; the radiometric half wants a short flight, because the sensor is warming throughout.

This page covers the chain from radiometric frames to a temperature orthomosaic with a stated accuracy: the radiance-to-temperature conversion and what it assumes, drift estimation and correction, mosaicking that preserves the measurement, and the geometric strategy that avoids reconstructing from thermal imagery at all. It is the thermal branch of multispectral, thermal and index mapping pipelines.

Audience and prerequisites. Python 3.10+, radiometric thermal frames — the colourised JPEG output of a consumer camera is not usable — and, ideally, one or more ground reference temperatures measured during the flight.

# Prerequisites

# Conceptual architecture

The sensor measures radiance in a broad thermal band. Converting that to the temperature of the surface requires four quantities the sensor does not know: the surface’s emissivity, the temperature of whatever is being reflected off it — usually the sky — the atmospheric transmission over the path, and the temperature of the air in that path.

The sensitivity to each is worth internalising because it determines where effort belongs. Emissivity dominates: an error of 0.02 in emissivity is roughly 1 °C on a typical scene. Reflected temperature matters most on low-emissivity surfaces such as metal roofs, where a large share of what the sensor sees is reflected sky. Atmospheric terms are small at drone altitudes — under a hundred metres of path, transmission is above 0.98 — and are frequently and defensibly ignored.

Figure 1 — Where a degree comes from. The ordering decides what is worth measuring on site.

# Step 1: Extract the radiometric data

A radiometric thermal JPEG embeds the raw values in a vendor tag alongside the visible preview. Reading the preview instead of the embedded data is the single most common mistake, and it produces a file of display values that look like temperatures and are not.

import numpy as np
import exiftool


def read_radiometric(path: str) -> dict:
    """Raw thermal values and the conversion parameters from a radiometric JPEG.

    The embedded raw data is a separate stream from the visible image. A
    reader that decodes the JPEG normally gets the colourised preview, whose
    values are display levels rather than radiance.
    """
    with exiftool.ExifToolHelper() as et:
        tags = et.get_metadata(path)[0]
        raw = et.execute("-b", "-RawThermalImage", path, raw_bytes=True)

    def get(name, default=None):
        hit = next((v for k, v in tags.items() if k.endswith(name)), default)
        if hit is None:
            raise KeyError(f"{path}: missing thermal parameter {name}")
        return hit

    return {
        "raw": raw,
        "emissivity": float(get("Emissivity", 0.95)),
        "reflected_temp_c": float(get("ReflectedApparentTemperature", 20.0)),
        "planck": {k: float(get(f"Planck{k.upper()}")) for k in ("r1", "b", "f", "o", "r2")},
        "datetime": get("DateTimeOriginal"),
    }

# Step 2: Convert radiance to temperature

The conversion is the sensor’s Planck calibration, inverted, with the emissivity and reflected-temperature correction applied to the radiance first.

import numpy as np


def raw_to_celsius(raw: np.ndarray, *, planck: dict, emissivity: float,
                   reflected_temp_c: float) -> np.ndarray:
    """Invert the sensor's Planck calibration to surface temperature.

    The reflected component is removed from the measured radiance before
    inversion, which is why a low-emissivity surface is so sensitive to the
    reflected temperature: a large share of what the sensor measured came
    from the sky rather than from the surface.
    """
    r1, b, f, o, r2 = (planck[k] for k in ("r1", "b", "f", "o", "r2"))

    def radiance_of(temp_c: float) -> float:
        return r1 / (r2 * (np.exp(b / (temp_c + 273.15)) - f)) - o

    measured = raw.astype(np.float64)
    reflected = radiance_of(reflected_temp_c)
    surface = (measured - (1.0 - emissivity) * reflected) / emissivity

    # Invert: solve the Planck expression for temperature.
    inner = r1 / (r2 * (surface + o)) + f
    with np.errstate(invalid="ignore", divide="ignore"):
        kelvin = b / np.log(np.maximum(inner, 1e-9))
    celsius = kelvin - 273.15
    return np.where(np.isfinite(celsius), celsius, np.nan).astype("float32")

# Step 3: Estimate and remove sensor drift

An uncooled microbolometer has no active temperature stabilisation, so its response changes as the airframe warms. Over a twenty-minute flight the shift is commonly two to three degrees, and mosaicking averages it into visible banding.

import numpy as np


def estimate_drift(times_s: np.ndarray, reference_temps_c: np.ndarray,
                   *, expected_real_trend_c_per_hour: float = 0.0) -> dict:
    """Sensor drift from repeated observations of a stable ground target.

    A concrete pad observed at intervals through a flight warms slowly and
    genuinely; subtracting an expected physical trend leaves the sensor's own
    drift, which is what the correction should remove.
    """
    A = np.column_stack([times_s, np.ones_like(times_s)])
    (slope, intercept), *_ = np.linalg.lstsq(A, reference_temps_c, rcond=None)
    physical = expected_real_trend_c_per_hour / 3600.0
    drift_per_s = slope - physical

    residual = reference_temps_c - (slope * times_s + intercept)
    return {"drift_c_per_min": float(drift_per_s * 60),
            "total_drift_c": float(drift_per_s * (times_s.max() - times_s.min())),
            "residual_sd_c": float(np.std(residual, ddof=1)),
            "observations": int(times_s.size)}


def apply_drift_correction(temps_c: np.ndarray, frame_time_s: float,
                           reference_time_s: float, drift_c_per_min: float) -> np.ndarray:
    """Refer every frame back to the sensor's state at a reference time."""
    minutes = (frame_time_s - reference_time_s) / 60.0
    return temps_c - drift_c_per_min * minutes

The reference target is what makes this possible, and it does not need to be elaborate: a square metre of concrete, a water tray, or a painted board, flown over at the start, middle and end of the survey. Correcting thermal drift across a flight covers the practice in detail.

Figure 2 — The drift, measured from a square metre of concrete flown over three times.

# Step 4: Get the geometry from the RGB, not the thermal

Thermal frames reconstruct poorly. The resolution is low, the texture is weak, and matching between frames is unreliable — which is why so many thermal mosaics have geometric errors far larger than their radiometric ones.

The robust approach is to fly an RGB camera alongside, reconstruct from the RGB imagery using the whole of the rest of this site’s machinery, and project the thermal frames through the resulting camera model. That requires a rigid transform between the two cameras, which is measured once and reused, and gives the thermal product the geometric quality of the RGB survey.

import numpy as np


def thermal_to_rgb_camera(thermal_pose: np.ndarray, extrinsic: np.ndarray) -> np.ndarray:
    """Place a thermal frame using the RGB camera's solved pose.

    `extrinsic` is the fixed transform between the two cameras on the rig,
    measured once from a calibration flight over a target visible in both.
    Reconstructing thermal frames independently would give two geometries
    that disagree by more than either's own error.
    """
    return thermal_pose @ extrinsic

The matching problem this introduces — establishing the extrinsic in the first place — is covered in matching thermal and RGB orthomosaics.

# Step 5: Mosaic without averaging the measurement away

General-purpose orthomosaic software blends overlapping frames to hide seams, which for a visual product is exactly right and for a temperature raster is destructive. Blending averages two frames that, after drift correction, should agree — and where they do not, it smears the disagreement across the overlap instead of exposing it.

Three settings make a mosaic a measurement.

No blending. Each output pixel takes its value from one frame, chosen by a seamline, rather than from a weighted average. The seams are then visible where the frames disagree, which is information rather than a defect.

No colour balancing. Any per-frame level adjustment is undoing the calibration by definition. It must be off.

Choose the frame, deliberately. Where several frames cover a pixel, prefer the one whose centre is closest — the frame that saw the pixel most nearly nadir, with the shortest path and the least off-axis falloff.

import numpy as np


def select_source_frame(candidates: list[dict], pixel_xy: tuple[float, float]) -> int:
    """Pick which frame supplies a mosaic pixel, for a measurement product.

    Nearest-to-centre rather than a blend: the value must come from one
    observation, so that a disagreement between frames shows as a seam
    instead of being averaged into something neither frame measured.
    """
    x, y = pixel_xy
    best, best_d = None, float("inf")
    for i, frame in enumerate(candidates):
        cx, cy = frame["centre"]
        d = np.hypot(x - cx, y - cy)
        if d < best_d:
            best, best_d = i, d
    return best


def seam_disagreement(frame_a: np.ndarray, frame_b: np.ndarray,
                      overlap: np.ndarray) -> dict:
    """How far apart two frames are over their shared area, in degrees."""
    a = frame_a[overlap & np.isfinite(frame_a) & np.isfinite(frame_b)]
    b = frame_b[overlap & np.isfinite(frame_a) & np.isfinite(frame_b)]
    if a.size < 200:
        return {"note": "insufficient overlap to judge"}
    diff = a - b
    return {"median_c": float(np.median(diff)),
            "p95_abs_c": float(np.percentile(np.abs(diff), 95)),
            "acceptable": bool(abs(np.median(diff)) < 0.5)}

Measuring the seam disagreement across every overlap and reporting its distribution is the strongest single quality statement a thermal mosaic can carry. A survey whose frames agree to within half a degree over their overlaps has a well-corrected drift; one where they disagree by two degrees has not, and the number says so without needing any ground truth.

# Step 6: State an accuracy, and say what it is an accuracy of

Thermal deliverables are routinely handed over with no accuracy statement, or with the sensor’s datasheet figure, which describes neither the conversion assumptions nor the drift. A defensible statement has three parts.

Absolute accuracy — how close a pixel is to the true surface temperature — is limited by emissivity and reflected-temperature assumptions, and is rarely better than ±2 °C on a mixed scene. It is the figure clients assume they are getting and the hardest to deliver.

Relative accuracy within the flight — how comparable two pixels in the mosaic are — is limited by the residual drift after correction, and is commonly ±0.5 °C. It is what most thermal applications actually need: finding the hottest panel in an array, or the coolest part of a field, is a relative question.

Repeatability between flights is limited by everything above plus the differences between the two days, and is the weakest of the three. Two thermal surveys a month apart are not comparable in absolute terms without ground references on both.

def accuracy_statement(seam_p95_c: float, reference_bias_c: float | None,
                       reference_spread_c: float | None) -> dict:
    """A three-part accuracy statement, derived from what was actually measured."""
    relative = round(max(seam_p95_c, 0.3), 1)
    absolute = (round(abs(reference_bias_c) + 2 * (reference_spread_c or 0.5), 1)
                if reference_bias_c is not None else None)
    return {
        "relative_within_flight_c": relative,
        "absolute_c": absolute,
        "absolute_note": ("from ground references" if absolute is not None
                          else "not measurable — no ground references on this flight"),
        "between_flight_note": "requires ground references on both flights",
    }

Writing the absolute figure as unavailable when no ground reference was flown is more useful than quoting a datasheet number. It tells the client what to commission next time, and it prevents a relative product being used for an absolute decision.

# Parameter deep-dive

# Verification and output inspection

import numpy as np


def thermal_accuracy_report(measured_c: dict[str, float],
                            reference_c: dict[str, float]) -> dict:
    """Compare mosaic temperatures against ground measurements."""
    common = sorted(set(measured_c) & set(reference_c))
    if len(common) < 2:
        return {"note": "at least two ground references are needed"}

    diffs = np.array([measured_c[k] - reference_c[k] for k in common])
    bias = float(np.mean(diffs))
    spread = float(np.std(diffs, ddof=1)) if diffs.size > 1 else float("nan")

    problems = []
    if abs(bias) > 1.5:
        problems.append(f"bias {bias:+.1f} °C — check emissivity and reflected temperature")
    if spread > 1.5:
        problems.append(f"spread {spread:.1f} °C — drift correction may be incomplete")
    return {"targets": common, "bias_c": bias, "spread_c": spread,
            "problems": problems}

Two ground references are the practical minimum, and they should differ in temperature — a shaded surface and a sunlit one — because a single reference constrains the offset and says nothing about the scale.

Figure 3 — Four stages, each adding an error the next cannot remove.

# Troubleshooting

The mosaic has stripes along flight lines. Sensor drift, uncorrected. Fit it from a reference target and apply per frame before mosaicking.

Every temperature is several degrees high or low. Emissivity. A default of 0.95 applied to water (0.99) or to bare metal (0.2) is wrong by degrees.

Metal roofs read close to sky temperature. Correct behaviour: a low-emissivity surface reflects the sky, and the sensor sees mostly that. Thermal measurement of bare metal is unreliable at any accuracy.

The thermal mosaic does not line up with the RGB. The extrinsic between the cameras is wrong or was not applied. Do not attempt to warp the finished mosaics into agreement; fix the projection.

Values are display levels rather than temperatures. The colourised preview was read instead of the embedded radiometric stream. The two live in the same file, and a normal JPEG decode returns the wrong one.

The geometry is poor and the reconstruction failed. Expected from thermal frames alone. Reconstruct from RGB and project the thermal through it, using the measured rigid transform between the two cameras.

Two flights over the same site disagree by several degrees. Almost certainly the conditions rather than the processing: wind speed, cloud, and time since sunrise all change surface temperatures by more than any calibration error. Absolute comparison between flights needs ground references on both days, and relative comparison within each flight is usually the better product.

# Related

• Converting Radiometric JPEGs to Temperature Rasters
• Correcting Thermal Drift Across a Flight
• Fixing Thermal Mosaic Seam Lines
• Matching Thermal and RGB Orthomosaics
• Radiometric Calibration in Python

← Multispectral, Thermal & Index Mapping Pipelines