Multispectral, Thermal & Index Mapping Pipelines

An RGB survey answers questions about shape. A multispectral or thermal survey answers questions about state — how much chlorophyll is in that canopy, where the irrigation is failing, which panel in the array is running hot. The geometry is the same problem solved in the rest of this site; the difference is that the pixel values now have to mean something physical, and almost everything in a standard photogrammetry pipeline is designed to make them look good instead.

That single distinction drives everything on this page. An orthomosaic assembled with colour balancing, exposure matching and seamline blending is exactly what you want for a visual product and exactly what destroys a reflectance measurement. A vegetation index computed from uncalibrated digital numbers varies with the cloud cover at the moment each frame was taken, which is why index values from two flights over the same field routinely disagree by more than the effect anybody is trying to measure.

This section treats radiometry as a first-class pipeline concern: calibrate before mosaicking, co-register the bands before computing anything across them, mask before summarising, and keep every step reproducible so a number from March is comparable with a number from September.

Figure 1 — The pipeline. The annotation on the orthomosaic stage is the one that most often has to be enforced against the default behaviour of general-purpose software.

# What the sensors actually record

The hardware shapes every decision downstream, so it is worth being precise about what a multispectral or thermal rig produces.

A multispectral camera is several narrow-band sensors in one housing, typically five: blue, green, red, red edge and near-infrared. Each has its own lens, its own optical axis and its own exposure. The bands are narrow — 20 to 40 nm — which is what makes them useful for discriminating vegetation state, and it also means each band collects far less light than an RGB channel, so exposures are longer and motion blur is a real constraint on flying speed.

A downwelling light sensor sits on top of the airframe and measures incident irradiance in the same bands, once per frame. It is what makes a flight under broken cloud recoverable: the irradiance record says how much light was available when each frame was taken, so the reflectance calculation can divide it out.

A reflectance panel is a calibrated Lambertian surface with a certificate giving its reflectance per band. Photographed before and after a flight, it ties the whole dataset to a physical scale. A panel alone is sufficient under stable light; under changing cloud it is a two-point calibration of a quantity that varied continuously, and the downwelling sensor is what fills the gap.

A radiometric thermal camera records not a picture but a per-pixel radiance from which temperature is derived, given assumptions about emissivity, reflected sky temperature, atmospheric transmission and distance. Consumer thermal cameras often output a colourised image with the radiometric data discarded, which is unusable for measurement; the radiometric variants embed the raw values, usually in a proprietary tag inside an otherwise ordinary JPEG.

Reading that table as a checklist for the flight is more effective than any processing step: most of the failures below are cheaper to prevent in the air than to correct on the ground.

# Why an uncalibrated index is not a measurement

A camera records digital numbers. What reaches the sensor is the product of the surface’s reflectance, the illumination falling on it, the atmosphere between, the lens transmission at that point in the frame, and the exposure the camera chose. Only the first of those is the quantity of interest, and the others vary between frames of the same flight.

The most significant is illumination. A flight under broken cloud has frames taken in full sun and frames taken in shadow, and the digital numbers differ by a factor of two or more for identical ground. An index built as a ratio of two bands cancels some of this — which is why ratio indices are used at all — but only if the two bands were exposed identically, and they were not.

Calibration converts digital numbers to reflectance: a dimensionless, physical property of the surface that does not depend on the light. That is what makes a value from one flight comparable with a value from another, and it is the entire justification for the reflectance panel and the downwelling light sensor that multispectral rigs carry.

import numpy as np


def digital_number_to_reflectance(dn: np.ndarray, *, dark: float,
                                  exposure_s: float, gain: float,
                                  vignette: np.ndarray,
                                  irradiance: float,
                                  panel_factor: float) -> np.ndarray:
    """The full chain from a raw frame to surface reflectance.

    Each term removes one thing the camera introduced. Dark current is the
    sensor's own floor; exposure and gain are the camera's choices; the
    vignette map is the lens; irradiance is the light at that instant; and
    the panel factor ties the result to a known reflectance standard.
    """
    corrected = (dn.astype(np.float64) - dark) / (exposure_s * gain)
    corrected /= np.maximum(vignette, 1e-6)
    radiance = corrected
    return panel_factor * radiance / max(irradiance, 1e-9)

Every term in that function is covered in radiometric calibration in Python, and the panel step in particular in applying reflectance panel calibration in Python.

# Band alignment: a sub-pixel problem with visible consequences

Most multispectral rigs use a separate sensor per band, physically offset by a few centimetres. At 60 m altitude those offsets project to several pixels on the ground, and the offset depends on the distance to the surface — so a single calibration cannot correct it everywhere in a scene with relief.

The consequence appears in any index that combines bands. Where red and near-infrared are misaligned by two pixels, every edge in the scene produces a bright or dark fringe in the index: a hedge acquires a halo, a bare patch acquires a rim. Those artefacts are then dutifully included in per-plot statistics.

Figure 2 — A two-pixel offset, and the fringe it writes into every index raster.

The alignment workflow — estimating a per-frame homography from matched features, checking the residual, and refitting where it fails — is covered in band alignment and stacking for multispectral sets.

# Indices: choosing one, and computing it safely

An index is a ratio, and ratios have two failure modes that ordinary raster arithmetic handles badly: division by values near zero, and propagation of NoData. Both produce values that are not merely wrong but unbounded, which then dominate any statistic computed over them.

import numpy as np


def safe_ratio_index(a: np.ndarray, b: np.ndarray, *, nodata_mask: np.ndarray,
                     min_sum: float = 0.02) -> np.ndarray:
    """A normalised difference index with the two failure modes handled.

    The denominator floor is not arbitrary: a pixel whose two bands sum to
    almost nothing is a shadow or a NoData edge, and the index there is
    meaningless rather than extreme. Masking it is more honest than clipping.
    """
    a = a.astype("float32")
    b = b.astype("float32")
    total = a + b
    valid = (~nodata_mask) & np.isfinite(total) & (np.abs(total) > min_sum)
    out = np.full(a.shape, np.nan, dtype="float32")
    np.divide(a - b, total, out=out, where=valid)
    return out

Which index to choose is a question about what is being measured and what the sensor can see: choosing between NDVI, NDRE and GNDVI works through the common cases, and writing a safe NDVI with NoData handling covers the arithmetic in full.

# Thermal: a different measurement wearing the same clothes

A thermal orthomosaic looks like any other raster and behaves completely differently. The pixel values are temperatures derived from radiance, which depend on the emissivity of the surface, the reflected temperature of the sky, the atmospheric path and — crucially — the sensor’s own temperature, which drifts through a flight as the airframe warms.

That drift is the defining problem. An uncooled microbolometer can move by two or three degrees over a twenty-minute flight, so the first and last frames of a survey disagree about a surface that never changed. Mosaicking averages that into visible banding, and any absolute temperature read from the result is wrong by an amount nobody can estimate after the fact.

import numpy as np


def drift_from_repeat_passes(times_s: np.ndarray, temps_c: np.ndarray) -> dict:
    """Estimate sensor drift from repeated observations of a stable target.

    A patch of hardstanding observed at the start, middle and end of a flight
    should read the same temperature, allowing for its own slow warming. What
    remains after removing a physical trend is the sensor's drift.
    """
    A = np.column_stack([times_s, np.ones_like(times_s)])
    (slope, intercept), *_ = np.linalg.lstsq(A, temps_c, rcond=None)
    residual = temps_c - (slope * times_s + intercept)
    return {"drift_c_per_min": float(slope * 60.0),
            "residual_sd_c": float(np.std(residual, ddof=1)),
            "total_drift_c": float(slope * (times_s.max() - times_s.min()))}

The full treatment, including converting radiometric JPEGs to temperature and correcting the drift across a mosaic, is in thermal orthomosaic processing in Python.

# Masking before summarising

A per-plot index statistic is only as meaningful as the pixels it averages. A plot polygon drawn around a crop row contains crop, soil between the rows, shadow cast by the canopy, and often a strip of track at the edge. Averaging all of it produces a number that moves with row spacing and sun angle as much as with crop health.

Masking is therefore not an optional refinement — it is what makes the statistic a measurement of the crop rather than of the plot’s geometry. Three masks cover most cases, and all three are cheap.

A vegetation mask removes soil. The usual construction is a threshold on the index itself, which sounds circular and is not: an NDVI below about 0.2 is bare ground under any reasonable calibration, and excluding it leaves the distribution of the canopy.

A shadow mask removes pixels where the illumination differed from the rest of the plot. Shadowed canopy has a genuinely different reflectance spectrum, and its inclusion widens the distribution without adding information about the plant.

A boundary buffer removes the plot’s own edge, where mixed pixels and any residual band misalignment concentrate. Buffering inward by two or three pixels typically removes a few percent of the area and a disproportionate share of the variance.

import numpy as np


def canopy_mask(ndvi: np.ndarray, nir: np.ndarray, *,
                veg_threshold: float = 0.2,
                shadow_percentile: float = 10.0) -> np.ndarray:
    """Pixels that are canopy, lit, and inside the plot's reliable interior.

    The shadow test uses the near-infrared band rather than brightness:
    vegetation is bright in NIR under any illumination, so a low NIR value
    inside a vegetated pixel means the pixel is shaded rather than sparse.
    """
    vegetated = np.isfinite(ndvi) & (ndvi > veg_threshold)
    if not vegetated.any():
        return vegetated
    cutoff = np.percentile(nir[vegetated], shadow_percentile)
    return vegetated & (nir > cutoff)


def plot_statistics(index: np.ndarray, mask: np.ndarray) -> dict:
    """Robust summary of an index over a masked plot."""
    values = index[mask & np.isfinite(index)]
    if values.size < 50:
        return {"pixels": int(values.size), "note": "too few valid pixels to summarise"}
    return {"pixels": int(values.size),
            "median": float(np.median(values)),
            "p10": float(np.percentile(values, 10)),
            "p90": float(np.percentile(values, 90)),
            "iqr": float(np.percentile(values, 75) - np.percentile(values, 25))}

Reporting the median and an interquartile range rather than a mean and a standard deviation is deliberate. Index distributions within a plot are routinely skewed — a patch of poor establishment pulls a long tail — and the median tracks the bulk of the crop while the spread describes the variability that agronomists actually act on. The full workflow is in computing zonal index statistics per plot polygon.

# Keeping two flights comparable

Everything above exists to serve one requirement: that a number from this flight can be compared with a number from the last one. Four practices deliver it, and they are largely organisational rather than technical.

Fly under similar conditions. Calibration removes much of the illumination effect and not all of it — sun angle changes the proportion of shadow within a canopy, which no radiometric correction addresses. Flying within an hour of solar noon, consistently, removes a source of variation that no amount of processing can.

Use the same index, computed the same way. An NDVI with a denominator floor of 0.02 and one with no floor produce different plot means on the same data. Pin the computation, version it, and record which version produced each result.

Keep a permanent target in every flight. A concrete pad, a painted board, an unchanging roof — anything whose reflectance does not vary seasonally. Its index value across flights is a direct, continuous check on the calibration, and it is the single most informative diagnostic available for a multi-flight programme.

Record the provenance with the numbers. Panel serial and certificate values, irradiance source, alignment residuals, index version, mask thresholds. A plot mean without those is a number; with them it is a measurement.

from datetime import datetime, timezone


def index_provenance(*, panel_id: str, panel_reflectance: dict[str, float],
                     irradiance_source: str, alignment_residual_px: float,
                     index_name: str, index_version: str,
                     mask_thresholds: dict) -> dict:
    """The record that makes a plot statistic comparable across flights."""
    return {
        "computed_at": datetime.now(timezone.utc).isoformat(),
        "panel": {"id": panel_id, "reflectance": panel_reflectance},
        "irradiance_source": irradiance_source,
        "alignment_residual_px": alignment_residual_px,
        "index": {"name": index_name, "version": index_version},
        "masks": mask_thresholds,
    }

# Where this section sits in the wider pipeline

Nothing here replaces the geometry. A multispectral survey still needs the flight planning, EXIF handling, ground control and reconstruction covered elsewhere on this site — a reflectance value is only useful if it is in the right place, and the positional work is identical to an RGB survey’s.

Three interactions are worth naming, because they are where the two halves of the pipeline disagree.

Overlap requirements are higher, not lower. Narrow-band frames have less texture than RGB, so feature matching is weaker, and the reconstruction needs more overlap to converge. The geometry side of that is in calculating optimal flight overlap for Python processing; the practical consequence is that a flight plan copied from an RGB job will under-perform.

The reconstruction usually runs on one band. Aligning five bands independently wastes effort and produces five slightly different geometries. The standard approach is to reconstruct from the band with the most texture — usually green or red — and project the others through the resulting camera model, which is also what guarantees the bands share a geometry rather than merely being close.

Export is a raster problem like any other. Once reflectance and index rasters exist, everything in DEM/DSM generation and raster export automation applies unchanged: Cloud-Optimized GeoTIFFs, correct NoData, overviews built with an averaging resampler, and a declared CRS. The one difference is that index rasters are float data with a meaningful zero, so a lossy codec is never acceptable.

# Parameter reference

# Failure modes and diagnostics

• Index values differ between two flights of an unchanged field. Symptom: the mean NDVI moves by 0.1 with no agronomic cause. Cause: uncalibrated digital numbers, or a panel reading taken under different light from the flight. Detect by comparing a permanent target — a concrete pad, a roof — across flights; its index should be stable.

def target_stability(index_by_flight: dict[str, float], tol: float = 0.03) -> dict:
    """A permanent target's index must not move between flights."""
    values = list(index_by_flight.values())
    spread = max(values) - min(values)
    return {"spread": spread, "stable": spread <= tol,
            "note": ("calibration is consistent across flights" if spread <= tol
                     else "calibration differs between flights — index values are "
                          "not comparable")}

• Bright or dark fringes along every edge. Symptom: hedges and plot boundaries have halos in the index raster. Cause: band misalignment. Remediate per fixing band misregistration artifacts in index rasters.

• A visible grid or banding in the reflectance mosaic. Symptom: rectangular patches of differing brightness. Cause: either per-frame exposure variation that calibration did not remove, or a mosaicking step that applied blending. Check the orthomosaic settings before the calibration.

• NaN or infinite values in the index. Symptom: statistics come back as NaN. Cause: division where the band sum is near zero, or NoData propagated as a sentinel rather than NaN. Handle as in diagnosing NaN and Inf in index rasters.

• Thermal mosaic shows stripes along flight lines. Symptom: alternating warm and cool bands matching the flight pattern. Cause: sensor drift, uncorrected. The remedy is a drift model fitted from repeat observations, per correcting thermal drift across a flight.

• Band metadata missing after processing. Symptom: a stacked GeoTIFF whose bands are unlabelled, so nobody knows which is red edge. Cause: a writer that did not carry band descriptions. Always set them, and assert them on read.

# Integration checklist

Wire the stages together for a production run by confirming each contract below. These render as interactive toggles.

• Capture panel images before and after every flight, and record which frames they are
• Apply dark current and exposure normalisation before anything else
• Apply a vignette correction derived from the sensor, not assumed
• Use downwelling irradiance where the flight had variable cloud
• Estimate band alignment per frame and gate on the residual
• Build the reflectance orthomosaic with colour balancing and blending disabled
• Compute indices with an explicit denominator floor and NaN NoData
• Mask soil, shadow and water before computing per-plot statistics
• Record emissivity and reflected temperature with every thermal product
• Fit and apply a thermal drift model from repeat observations of a stable target
• Write band descriptions into every stacked raster and assert them on read
• Track a permanent target’s index across flights as a continuous calibration check

One further habit is worth adopting on any programme that will run for more than a season: keep the raw frames. Calibration methods improve, panel certificates get re-measured, and an index definition occasionally turns out to have been computed with the wrong band pair. All of those are recoverable from raw frames and none of them is recoverable from a delivered index raster. Storage is cheap next to a re-flight, and a season of raw multispectral frames is a few terabytes — less than the imagery an RGB programme routinely keeps without anyone questioning it.

A pipeline built to these contracts produces numbers that mean the same thing in March and September, which is the only property that makes a multispectral programme worth flying more than once.

Figure 3 — Every stage before the index exists to make the index mean something.

# Related

• Radiometric Calibration in Python
• Band Alignment and Stacking for Multispectral Sets
• Computing Vegetation Indices with Rasterio
• Thermal Orthomosaic Processing in Python
• Troubleshooting Multispectral and Thermal Failures

← Drone Photogrammetry Pipelines