Ingesting Multispectral Drone Imagery in Python

Agronomic significance

Surface reflectance is the physically interpretable quantity that connects sensor readings to crop physiology. Raw DNs are instrument-specific; reflectance values are comparable across flights, sensors, and dates, making longitudinal analysis of crop stress, nitrogen status, or irrigation uniformity possible. Without calibrated reflectance, a single cloud passing between flights can shift apparent NDVI by 0.15–0.25 units — large enough to trigger false prescription zones in variable-rate application workflows.

2. Step-by-step implementation

Step 1: Discover and group per-band files

Multispectral flights output one file per spectral band per capture event. Group files by capture timestamp so the five bands from each moment can be processed together.

from pathlib import Path
from typing import Dict, List
import rasterio

BAND_ORDER = ["Blue", "Green", "Red", "RedEdge", "NIR"]

def discover_flight_bands(flight_dir: str) -> Dict[str, List[Path]]:
    """
    Group single-band GeoTIFFs by capture timestamp.

    Returns a dict keyed by timestamp string, with values being a list of
    Paths sorted to BAND_ORDER (Blue, Green, Red, RedEdge, NIR).
    Raises FileNotFoundError if no TIFFs are found.
    """
    flight_path = Path(flight_dir)
    tiffs = sorted(flight_path.glob("*.tif"))
    if not tiffs:
        raise FileNotFoundError(f"No .tif files found in {flight_dir}")

    band_groups: Dict[str, List[Path]] = {}
    for tiff in tiffs:
        with rasterio.open(tiff) as src:
            timestamp = src.tags().get("DateTimeOriginal") or src.tags(ns="TIFF").get(
                "DateTimeOriginal", tiff.stem[:15]
            )
        band_groups.setdefault(timestamp, []).append(tiff)

    # Sort each group into canonical band order (sensor-specific suffix pattern)
    for ts, paths in band_groups.items():
        def band_key(p: Path) -> int:
            name = p.stem.upper()
            for i, band in enumerate(["_1", "_2", "_3", "_4", "_5"]):
                if name.endswith(band):
                    return i
            return 99
        band_groups[ts] = sorted(paths, key=band_key)

    assert all(len(v) == 5 for v in band_groups.values()), (
        "Expected exactly 5 bands per capture — check for missing or extra files."
    )
    return band_groups

Note: MicaSense RedEdge-MX files use numeric suffixes (_1 through _5); DJI P4 Multispectral uses _B, _G, _R, _RE, _NIR. Adjust the band_key comparator for your sensor. Parrot Sequoia adds an _sunshine sensor file per capture — exclude those by filtering on file size or a naming pattern before calling this function.

Step 2: Radiometric calibration — DNs to surface reflectance

The standard field calibration uses a reflectance panel with known grey-card reflectance values (typically 0.41–0.57 depending on manufacturer). Measure the mean DN of the panel region in the panel photograph, then apply the linear scaling to the scene imagery.

import numpy as np

def calibrate_to_reflectance(
    raw_array: np.ndarray,
    panel_dn: float,
    panel_reflectance: float,
    dark_offset: float = 0.0,
    exposure_ratio: float = 1.0,
) -> np.ndarray:
    """
    Convert raw DNs to surface reflectance using panel calibration.

    Args:
        raw_array:         2-D uint16 array from rasterio.read(1).
        panel_dn:          Mean DN of the calibration panel ROI in the panel photo.
        panel_reflectance: Known reflectance of the panel (0.0–1.0).
        dark_offset:       Sensor dark current at the flight temperature.
        exposure_ratio:    panel_exposure / scene_exposure (1.0 if identical).

    Returns:
        float32 array with values clipped to [0.0, 1.0].
    """
    if panel_dn <= 0:
        raise ValueError("panel_dn must be > 0; check that the panel ROI is correctly specified.")

    corrected_panel_dn = panel_dn * exposure_ratio
    scale = panel_reflectance / (corrected_panel_dn - dark_offset)
    reflectance = (raw_array.astype(np.float32) - dark_offset) * scale
    return np.clip(reflectance, 0.0, 1.0)

Panel coefficients differ per band; call this function independently for each of the five channels. MicaSense provides per-band panel reflectance values in the panel certificate that ships with each sensor — store these in a JSON sidecar alongside the flight data for auditability.

Step 3: Spatial alignment — reproject and co-register bands

All five bands must share an identical affine transform, resolution, and CRS before stacking. Use the first band as the reference geometry and reproject the remaining four into it. To work with metric coordinates (required for accurate area and distance operations), reproject to the appropriate UTM zone rather than leaving data in WGS 84 (EPSG:4326).

from rasterio.warp import reproject, Resampling, calculate_default_transform
from rasterio.crs import CRS
from rasterio.transform import Affine
from typing import Tuple

def align_bands_to_utm(
    band_files: List[Path],
    target_epsg: int = 32610,
) -> Tuple[np.ndarray, Affine, CRS]:
    """
    Reproject and co-register a list of single-band GeoTIFFs into a common UTM CRS.

    Returns:
        stacked:   float32 ndarray of shape (n_bands, height, width)
        transform: Affine transform of the output grid
        crs:       Output CRS object
    """
    target_crs = CRS.from_epsg(target_epsg)

    # Derive output grid from first band
    with rasterio.open(band_files[0]) as ref:
        dst_transform, dst_width, dst_height = calculate_default_transform(
            ref.crs, target_crs, ref.width, ref.height, *ref.bounds
        )

    stacked = np.zeros((len(band_files), dst_height, dst_width), dtype=np.float32)

    for i, band_path in enumerate(band_files):
        with rasterio.open(band_path) as src:
            reproject(
                source=rasterio.band(src, 1),
                destination=stacked[i],
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=dst_transform,
                dst_crs=target_crs,
                resampling=Resampling.bilinear,
            )

    # Sanity check: all bands must share the same shape
    assert stacked.shape[1] == dst_height and stacked.shape[2] == dst_width, (
        "Unexpected shape after reprojection — likely a transform mismatch."
    )
    return stacked, dst_transform, target_crs

UTM zone selection: For Northern Hemisphere fields use EPSG:326xx (e.g. 32610 = UTM Zone 10N for California). Southern Hemisphere fields use EPSG:327xx. Pass the correct EPSG code for your flight location; see understanding CRS in precision agriculture for a zone-lookup approach.

Step 4: Stack, validate, and export

Write the calibrated, aligned stack as a tiled, LZW-compressed GeoTIFF. Tiling (256×256 blocks) is essential for windowed read access during downstream index computation on large orthomosaics.

def export_reflectance_stack(
    stacked: np.ndarray,
    transform: Affine,
    crs: CRS,
    output_path: Path,
    band_names: List[str] | None = None,
) -> None:
    """
    Write a calibrated float32 multi-band GeoTIFF with validation guards.

    Raises AssertionError if value range or shape invariants are violated.
    """
    if band_names is None:
        band_names = BAND_ORDER[: stacked.shape[0]]

    # Pre-export validation
    assert stacked.ndim == 3, "stacked must be 3-D (bands, height, width)"
    assert stacked.dtype == np.float32, "Expected float32 array"
    assert stacked.min() >= 0.0 and stacked.max() <= 1.0, (
        f"Reflectance out of [0,1]: min={stacked.min():.4f}, max={stacked.max():.4f}. "
        "Check calibration coefficients."
    )
    nan_count = np.isnan(stacked).sum()
    assert nan_count == 0, f"{nan_count} NaN values found — check resampling step."

    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "count": stacked.shape[0],
        "height": stacked.shape[1],
        "width": stacked.shape[2],
        "crs": crs,
        "transform": transform,
        "compress": "lzw",
        "nodata": 0.0,
        "tiled": True,
        "blockxsize": 256,
        "blockysize": 256,
        "interleave": "band",
    }

    output_path.parent.mkdir(parents=True, exist_ok=True)
    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(stacked)
        for i, name in enumerate(band_names, start=1):
            dst.set_band_description(i, name)

    print(
        f"Exported {output_path.name} | "
        f"shape={stacked.shape} | CRS={crs.to_epsg()} | "
        f"reflectance range=[{stacked.min():.4f}, {stacked.max():.4f}]"
    )

Putting it together

from pathlib import Path

PANEL_COEFFICIENTS = {
    # Per-band: (panel_dn, panel_reflectance, dark_offset)
    # Values from MicaSense factory certificate — replace with your panel's actual values
    "Blue":    (38200, 0.53, 4100),
    "Green":   (39500, 0.53, 4200),
    "Red":     (37800, 0.53, 4050),
    "RedEdge": (36100, 0.53, 4000),
    "NIR":     (40200, 0.53, 4300),
}

def ingest_flight(flight_dir: str, output_dir: str, target_epsg: int = 32610) -> None:
    band_groups = discover_flight_bands(flight_dir)
    output_path = Path(output_dir)

    for timestamp, band_files in band_groups.items():
        # Step 1: Align bands spatially
        stacked, transform, crs = align_bands_to_utm(band_files, target_epsg)

        # Step 2: Calibrate each band independently
        calibrated = np.empty_like(stacked)
        for i, band_name in enumerate(BAND_ORDER):
            panel_dn, panel_refl, dark_offset = PANEL_COEFFICIENTS[band_name]
            calibrated[i] = calibrate_to_reflectance(
                stacked[i], panel_dn, panel_refl, dark_offset
            )

        # Step 3: Export validated stack
        safe_ts = timestamp.replace(":", "-").replace(" ", "_")
        export_reflectance_stack(
            calibrated, transform, crs,
            output_path / f"reflectance_{safe_ts}.tif"
        )

3. Key parameters & tuning

4. Handling edge cases & failure modes

Missing or partial band captures. Multi-lens sensors occasionally drop a band if one detector triggers at a different moment (e.g. due to GPS-triggered capture jitter on DJI P4 Multispectral). The discover_flight_bands assertion (len(v) == 5) will catch this early. Log the timestamp and skip it rather than proceeding with a 4-band stack, which would shift all band indices and corrupt any code that relies on band_names.

UTM zone boundary crossings. Fields that straddle a UTM zone boundary (every 6° of longitude) will have different source CRS values for images captured on each side. Always derive target_epsg from the centroid of the flight area rather than from the first image’s CRS. Use pyproj.database.query_utm_crs_info with the field centroid to determine the correct zone programmatically.

Calibration panel captured at different exposure. Forgetting to set exposure_ratio when the pilot used a different exposure for the panel photograph (common on overcast days) causes band-specific offsets of 0.05–0.15 in reflectance. Check panel photograph metadata before every batch run and compute exposure_ratio = panel_exp_time / scene_exp_time.

NaN values introduced by bilinear resampling at edges. During reprojection, pixels outside the source raster bounds are set to nodata (0.0 in this pipeline). Bilinear resampling creates a partial-data transition zone of 1–2 pixels around the raster edge. These pixels have reflectance values below the physical minimum and should be masked before index computation. Apply a 2-pixel inward buffer when clipping to field boundaries.

Orthomosaic vs per-capture workflow. This pipeline processes per-capture single-band images. If your photogrammetry software (Agisoft Metashape, DJI Terra) has already produced a 5-band orthomosaic GeoTIFF, skip Steps 1–3 and apply radiometric calibration directly to each band of the mosaic. Verify band order against the software’s export documentation — Metashape reverses the MicaSense band order in some export configurations.

5. Verification & output validation

After running the pipeline, confirm the output stack is correct before passing it downstream.

def validate_reflectance_stack(stack_path: Path) -> None:
    """Spot-check a multi-band reflectance GeoTIFF for common ingestion errors."""
    with rasterio.open(stack_path) as src:
        assert src.count == 5, f"Expected 5 bands, got {src.count}"
        assert src.crs is not None, "CRS is None — geospatial metadata missing"
        assert src.crs.is_projected, (
            f"CRS {src.crs.to_epsg()} is geographic, not projected. "
            "Reproject to a UTM zone for metric analysis."
        )
        assert src.nodata == 0.0, f"Unexpected nodata value: {src.nodata}"

        # Sample a 256x256 window from the image centre
        cx, cy = src.width // 2, src.height // 2
        window = rasterio.windows.Window(cx - 128, cy - 128, 256, 256)
        sample = src.read(window=window).astype(np.float32)

        # Exclude nodata pixels
        valid = sample[sample > 0.0]
        assert valid.min() >= 0.0 and valid.max() <= 1.0, (
            f"Reflectance out of [0,1]: [{valid.min():.4f}, {valid.max():.4f}]"
        )
        print(f"PASS: {stack_path.name} | valid pixel range [{valid.min():.4f}, {valid.max():.4f}]")

Visual spot-check. Load the output in QGIS or use rasterio.plot with a false-colour composite (NIR=R, Red=G, Green=B) — healthy vegetation should appear vivid red to pink. A blue or grey cast on crops indicates the NIR and Red bands may be swapped.

Histogram inspection. The Red band histogram for a mid-season crop field should peak between 0.05 and 0.15 (vegetation absorbs red light). NIR should peak between 0.35 and 0.55. Values outside these ranges are a calibration warning, not necessarily an error — but document the deviation before proceeding to index computation.

Calibration target assertion. If your flight included a calibration panel in the scene, open the exported GeoTIFF and read the mean reflectance over the panel pixels. It should match panel_reflectance ± 0.02. A larger deviation indicates an exposure-ratio error or an incorrect dark offset.

6. Integration with the broader pipeline

The output of this ingestion stage is a tiled float32 GeoTIFF with physically calibrated reflectance values and a projected CRS — the standard input format expected by every downstream stage.

Immediately downstream: Pass the reflectance stack to band math and vegetation index computation, where NDVI, NDRE, and SAVI indices are derived from the same five-band structure. Raw digital numbers must be converted to reflectance before applying band math — this ingestion step is what makes that possible.

Spatial masking: Before computing indices across a full orthomosaic, clip the raster to field boundaries using field boundary extraction with GeoPandas. This eliminates road edges, tree rows, and bare soil border strips that would skew canopy-level statistics.

Satellite interoperability: This pipeline is specific to drone sensors. If you need to work with Sentinel-2 imagery alongside drone data, the band nomenclature and atmospheric correction approach differ significantly — see parsing Sentinel-2 vs drone multispectral bands in Python for the comparison.

Cloud masking for satellite inputs: When drone imagery is combined with satellite time-series for gap-filling or validation, apply cloud masking at the satellite ingest stage — see automating cloud removal in Sentinel-2 time series for the SCL-based masking approach.

CRS consistency check across the pipeline. If downstream workflows consume both drone stacks and vector field boundaries from a farm management system, verify that both are in the same projected CRS before spatial operations. Use the approach described in validating coordinate systems for variable-rate maps to assert CRS equality before any raster–vector overlay.

This page is part of Ag-GIS Data Fundamentals & Spatial Reference Systems — the full pipeline context, including orthomosaic stitching and CRS fundamentals, is covered there.

Frequently asked questions

Why are my NDVI values above 1.0 after stacking drone bands?

Reflectance values above 1.0 almost always indicate a calibration error: the panel DN used in the denominator was read from the wrong capture, or the dark-current offset was not subtracted before dividing. Verify that your panel_dn value comes from a panel photograph taken during the same flight and that dark_offset matches the sensor’s reported dark current at that temperature.

How do I handle flights where the calibration panel was captured at the wrong exposure?

Apply an exposure-correction factor before computing the panel DN. MicaSense RedEdge-MX and DJI P4 Multispectral both embed exposure time in XMP metadata. Normalize the panel DN to a reference exposure: panel_dn_normalized = panel_dn * (ref_exposure / actual_exposure). Then use panel_dn_normalized in the calibration equation.

My bands are misaligned even after reprojecting to the same CRS — what is wrong?

CRS alignment alone does not fix physical parallax between lens positions. After reprojecting with rasterio.warp.reproject, a feature-matching step corrects sub-pixel offsets between spectral channels. MicaSense provides an alignment SDK in their imageutils Python library. Alternatively, use OpenCV SIFT keypoint matching with a homography to register each band to the reference before stacking.

Related

• Parsing Sentinel-2 vs Drone Multispectral Bands in Python — band nomenclature differences and atmospheric correction when mixing satellite and drone data
• Band Math & Raster Algebra in Python — NDVI, NDRE, and SAVI calculation on the reflectance stack produced here
• Field Boundary Extraction with GeoPandas — clip raster outputs to agricultural field polygons before index aggregation
• Understanding CRS in Precision Agriculture — UTM zone selection and projection validation for farm-scale spatial analysis
• Orthomosaic Stitching Workflows — upstream photogrammetry step that assembles per-capture images into a single orthomosaic