Troubleshooting Multispectral and Thermal Failures

Multispectral and thermal pipelines fail quietly more often than any other part of a photogrammetry stack. A geometric failure is loud — the reconstruction breaks, the mosaic is visibly wrong — while a radiometric failure produces a smooth raster of plausible numbers that are wrong by a constant factor.

That asymmetry shapes how to look for problems. The diagnostic order here is not “what broke” but “what would look exactly like this if it were wrong”, and almost every check is a comparison against something independent: a permanent target, a second band, a previous flight.

This page aggregates the failure modes from across multispectral, thermal and index mapping pipelines into one sequence, with the cheap checks first.

Audience and prerequisites. Python 3.10+, a suspect product, and ideally a permanent ground target in the survey. Without an independent reference, most of what follows becomes an argument rather than a measurement.

# Prerequisites

# Conceptual architecture

The cascade runs from metadata to values to comparisons, cheapest first.

Metadata is free and catches the largest class: missing band descriptions, an absent formula tag, a calibration flag that says the product is uncalibrated when the report claims otherwise.

Distribution is cheap and catches the second largest: values outside a physical range, a spike at exactly ±1, a masked fraction that has changed since last month.

Independent comparison is the expensive one and catches what the others cannot: a permanent target whose value moved, a band pair that disagrees where they should agree, two flights whose composition differs.

Figure 1 — The order, and the reason the third level cannot be skipped on a calibrated product.

# Step 1: The metadata check

import rasterio


REQUIRED_INDEX_TAGS = ("FORMULA", "MIN_SUM", "UNITS")


def metadata_diagnostics(path: str, *, expect_bands: list[str] | None = None,
                         is_index: bool = False) -> dict:
    """Everything wrong that can be seen without reading pixels."""
    problems = []
    with rasterio.open(path) as src:
        descriptions = [d or "" for d in src.descriptions]
        band_tags = {i: src.tags(i) for i in range(1, src.count + 1)}
        file_tags = src.tags()
        dtype = src.dtypes[0]
        nodata = src.nodata

    if not any(descriptions):
        problems.append("no band descriptions — this product cannot be read by name")
    if expect_bands and descriptions != expect_bands:
        problems.append(f"band order is {descriptions}, expected {expect_bands}")
    if dtype != "float32":
        problems.append(f"dtype is {dtype} — reflectance and indices should be float32")
    if nodata is not None and nodata == 0:
        problems.append("nodata is 0, which is a valid reflectance value")

    if is_index:
        missing = [t for t in REQUIRED_INDEX_TAGS if t not in band_tags.get(1, {})]
        if missing:
            problems.append(f"index tags absent: {missing} — not comparable with anything")

    return {"bands": descriptions, "dtype": dtype, "nodata": nodata,
            "file_tags": file_tags, "problems": problems}

The missing-tags branch is the one that determines whether a product has a future. An index raster without its formula and denominator floor cannot be compared with next month’s, and that is not recoverable later — the parameters are gone.

# Step 2: The distribution check

import numpy as np
import rasterio


def distribution_diagnostics(path: str, *, kind: str = "index") -> dict:
    """Value-range and shape checks that catch the common silent failures."""
    with rasterio.open(path) as src:
        arr = src.read(1, masked=True).filled(np.nan)

    finite = arr[np.isfinite(arr)]
    if finite.size == 0:
        return {"problems": ["no finite values at all"]}

    problems = []
    masked = 1 - finite.size / arr.size

    if kind == "index":
        if np.any(np.abs(finite) > 1.0001):
            problems.append("values outside [-1, 1] — a band was negative")
        at_extreme = float(np.count_nonzero(np.abs(np.abs(finite) - 1.0) < 1e-4)
                           / finite.size)
        if at_extreme > 0.001:
            problems.append(f"{at_extreme:.2%} of pixels sit at exactly ±1 — "
                            "NoData is being treated as data")
    elif kind == "reflectance":
        over = float(np.count_nonzero(finite > 1.1) / finite.size)
        if over > 0.01:
            problems.append(f"{over:.1%} of pixels exceed reflectance 1.1 — "
                            "check the panel measurement")
        if float(np.count_nonzero(finite < 0) / finite.size) > 0.01:
            problems.append("a percent or more of pixels are negative — dark level")
    elif kind == "temperature":
        if finite.max() - finite.min() > 120:
            problems.append("temperature range exceeds 120 °C — conversion fault")

    if masked > 0.3:
        problems.append(f"{masked:.0%} of the raster is masked")

    return {"masked_fraction": masked, "median": float(np.median(finite)),
            "p01": float(np.percentile(finite, 1)),
            "p99": float(np.percentile(finite, 99)), "problems": problems}

# Step 3: The independent comparison

This is the level that catches a correctly formed product with wrong values, and it needs something the pipeline did not produce.

import numpy as np


def permanent_target_check(history: dict[str, float], *, tolerance: float = 0.03,
                           current_flight: str | None = None) -> dict:
    """A surface that cannot have changed must not change.

    This single check catches panel errors, calibration drift, band mix-ups
    and firmware changes, none of which any amount of internal consistency
    checking would reveal.
    """
    flights = sorted(history)
    values = np.array([history[f] for f in flights], dtype=float)
    if values.size < 2:
        return {"note": "need at least two flights to compare"}

    baseline = float(np.median(values[:-1])) if current_flight else float(np.median(values))
    latest = float(values[-1])
    drift = latest - baseline
    spread = float(values.max() - values.min())

    problems = []
    if abs(drift) > tolerance:
        problems.append(f"the permanent target moved {drift:+.3f} on the latest flight")
    if spread > 2 * tolerance:
        problems.append(f"target spans {spread:.3f} across the series — the archive "
                        "is not internally comparable")
    return {"flights": len(flights), "latest": latest, "baseline": baseline,
            "drift": drift, "spread": spread, "problems": problems}

# Step 4: Tracing a fault back to its stage

When a check fires, the next question is which stage produced the fault, and the chain is short enough to bisect systematically. Each stage has a signature at its own output, so keeping one intermediate product per stage turns a whole-pipeline mystery into a single comparison.

After calibration, reflectance should lie in [0, 1] with only a fraction of a percent above unity, and a permanent target should read close to its known value. A failure here is the panel, the vignette model or the irradiance record.

After band alignment, the bands’ strong edges should coincide within a pixel. A failure here is feature matching or a fixed calibration applied at the wrong altitude.

After mosaicking, frames should agree over their overlaps. A failure here is exposure drift, or a mosaicking step that applied blending or colour balancing.

After index computation, values should lie in [−1, 1] with no spike at the extremes. A failure here is NoData handling or a denominator floor.

def trace_fault(stage_reports: dict[str, dict]) -> dict:
    """Identify the earliest stage whose output already showed the problem.

    Ordered by position in the pipeline, because a fault at an early stage
    propagates to every later one and reporting all of them buries the cause.
    """
    order = ["calibration", "alignment", "mosaic", "index"]
    for stage in order:
        report = stage_reports.get(stage)
        if report and report.get("problems"):
            return {"first_failing_stage": stage,
                    "problems": report["problems"],
                    "note": "later stages inherit this; fix here first"}
    return {"first_failing_stage": None, "note": "all stage outputs are clean"}

Keeping the intermediates costs storage and saves hours. A pipeline that writes only its final product forces every investigation to start from a full reprocessing run.

# Step 5: The failures that are not faults

Several things that look like failures are the pipeline working correctly, and recognising them saves a great deal of wasted investigation.

A large masked fraction over water. Water has near-zero near-infrared reflectance, so a normalised difference index there has a denominator near zero and is correctly masked. The mask is the right answer.

Metal roofs reading close to sky temperature. Low emissivity means the sensor sees mostly reflected sky. No processing makes a thermal measurement of bare metal reliable.

No ground under closed canopy in a co-flown lidar comparison. Photogrammetry has no data there, as described in the point-cloud section; the absence is real.

An index that stops discriminating at canopy closure. NDVI saturating is physics rather than a fault, and the remedy is a different index rather than a different pipeline.

Plot medians that move when the masking policy changed. Two products computed with different masks are different products. The fault, if any, is that the change was not recorded.

EXPECTED_BEHAVIOURS = {
    "water_masked": "near-zero NIR gives a near-zero denominator; masking is correct",
    "metal_reads_sky": "low emissivity means the sensor sees reflected sky",
    "ndvi_saturates": "red is fully absorbed at canopy closure; use a red-edge index",
    "canopy_no_ground": "photogrammetry has no observation under closed canopy",
}


def is_expected(symptom: str) -> str | None:
    """Return the explanation when a reported symptom is correct behaviour."""
    return EXPECTED_BEHAVIOURS.get(symptom)

Keeping that list in the codebase rather than in institutional memory is worth doing. Every one of these is reported as a bug at least once per new team member, and a one-line explanation in the diagnostic output ends the conversation before it starts.

# Step 6: What to record so the next failure is quick

The recurring theme across this section is that radiometric faults are diagnosed by comparison, and comparison needs history. Four things, recorded per flight, make almost every investigation short.

The permanent target’s index and reflectance, per band, per flight. This is the single most valuable number in the archive and it costs nothing to extract.

The calibration parameters actually used — panel identity, certificate values, irradiance source, vignette model version. Not the intent, the actuals.

The mask composition per plot: what share was soil, shadow, boundary and canopy. A change here explains a great many apparent crop changes.

The stage reports from the cascade above, stored rather than printed. A rising masked fraction over three months is invisible in any single run and obvious in a series.

from datetime import date


def flight_record(target_values: dict, calibration: dict,
                  mask_fractions: dict, stage_reports: dict) -> dict:
    """The per-flight record that makes future diagnosis a lookup."""
    return {"flight_date": date.today().isoformat(),
            "permanent_target": target_values,
            "calibration": calibration,
            "mask_fractions": mask_fractions,
            "stage_reports": {k: v.get("problems", []) for k, v in stage_reports.items()}}

A season of those records turns “why is this month different” from an investigation into a query, which is the difference between an afternoon and a minute.

# A worked diagnosis

A client reports that this month’s NDVI is uniformly lower across a trial and asks whether the crop has been set back by the cold week.

The metadata check passes: bands are described, the formula and denominator floor match last month’s, the product is calibrated. So it is not a processing-version difference, which rules out the easiest explanation in one free step.

The distribution check shows a median of 0.58 against 0.71 last month, no out-of-range values and no spike at the extremes. The values are well formed; they are simply lower. That rules out NoData handling.

The permanent target check fires. The concrete pad reads 0.11 this month against 0.04 in the previous five flights — a shift of 0.07, well outside the 0.03 tolerance. The target cannot have changed, so the pipeline did.

Looking at the recorded calibration parameters, the panel identity is the same but the certificate values differ: somebody re-entered them from a newer certificate for a different panel of the same model. Every reflectance this month is scaled by the ratio between the two certificates.

Elapsed time: about four minutes, none of it spent looking at the crop. The cold week may still have had an effect, and it will be measurable once the calibration is corrected and the product re-issued — which is now a rerun rather than an investigation.

# Parameter deep-dive

# Verification and output inspection

Bundling the cascade into a gate is what converts these from debugging tools into something that runs on every product.

def gate_product(path: str, *, kind: str, expect_bands=None,
                 target_history: dict | None = None) -> dict:
    """Run the whole cascade and fail the job on anything it finds."""
    report = {"metadata": metadata_diagnostics(path, expect_bands=expect_bands,
                                               is_index=(kind == "index"))}
    if report["metadata"]["problems"]:
        raise ValueError("; ".join(report["metadata"]["problems"]))

    report["distribution"] = distribution_diagnostics(path, kind=kind)
    if report["distribution"]["problems"]:
        raise ValueError("; ".join(report["distribution"]["problems"]))

    if target_history:
        report["target"] = permanent_target_check(target_history)
        if report["target"].get("problems"):
            raise ValueError("; ".join(report["target"]["problems"]))
    return report

Raising rather than warning is the right default for a measurement product. A warning in a log is read once; an exception stops the delivery, which is what should happen when a survey’s absolute scale is in doubt.

Figure 3 — Triage forward, not backward from the symptom.

Figure 4 — Four ways to produce a file that is entirely valid and entirely wrong.

# Troubleshooting

Index statistics come back as NaN. A NoData sentinel survived into the arithmetic, or the denominator floor removed everything. See diagnosing NaN and Inf in index rasters.

A stack’s bands cannot be identified. Descriptions were never written or were dropped by an intermediate tool. See fixing missing or mismatched band metadata.

Calibration fails because the panel was not found. Detection thresholds, a panel too small in frame, or shading. See resolving panel detection failures in calibration.

Reflectance exceeds one over large areas. The panel radiance was measured too low — shading, or the holder included. Re-measure over the panel interior.

Everything looks correct and two flights still disagree. Check the permanent target, then the masking composition, then the sun angle. In that order, because that is their order of cost.

The thermal mosaic reads several degrees off. Emissivity, then reflected temperature, then drift. Nothing else in the chain is that large, so checking them in that order is almost always the shortest path.

A product passes every check and the client still disputes it. Establish what they are comparing it against. Most disputes at this stage are comparisons between products that were never comparable — a different sensor, a different index version, or a published threshold from another study.

# Related

• Diagnosing NaN and Inf in Index Rasters
• Fixing Missing or Mismatched Band Metadata
• Resolving Panel Detection Failures in Calibration
• Radiometric Calibration in Python
• Computing Vegetation Indices with Rasterio

← Multispectral, Thermal & Index Mapping Pipelines