Parsing Sentinel-2 vs Drone Multispectral Bands in Python

TL;DR: Use (DN - boa_offset) / 10000 for Sentinel-2 L2A reflectance (reading the offset from MTD_MSIL2A.xml) and map drone bands by name rather than positional index; both paths converge on the same float32 reflectance array in [0.0, 1.0] before any vegetation index calculation.

Why this task arises in ag-GIS workflows

Agronomic workflows routinely fuse Sentinel-2 satellite scenes — which cover entire farms at 10–60 m resolution — with drone orthomosaics at 2–5 cm resolution captured during the same growth stage. The goal is to cross-validate plant health signals or to fill cloud-gaps in the satellite time series with drone data. Without careful band parsing, the fusion pipeline silently produces wrong reflectance values: NDVI computed from unscaled Sentinel-2 DN (values near 3000–8000) will always return values close to 1.0, making healthy and stressed crop appear identical.

The problem compounds when comparing across platforms. Sentinel-2 L2A applies a radiometric offset that changed in January 2022, and drone orthomosaics from Pix4D, Agisoft Metashape, or WebODM embed calibration in ways that vary by firmware version. A single hard-coded parsing function breaks on either the old or new product format.

This page covers the exact parsing step within the broader ingesting multispectral drone imagery workflow. Once bands are parsed into valid reflectance, the output feeds directly into band math and raster algebra for NDVI, NDRE, and SAVI computation.

Prerequisites

This page assumes you have completed the environment setup in the parent guide. The only additional requirement here is:

  • rasterio >= 1.3
  • numpy >= 1.22
  • An XML parser from the standard library (xml.etree.ElementTree) for reading Sentinel-2 metadata

No additional packages are needed beyond what the parent cluster installs.

Platform differences at a glance

Sentinel-2 vs Drone Multispectral Band Architecture A side-by-side comparison table showing resolution, band count, file format, radiometric scaling, and band ordering for Sentinel-2 L2A versus drone platforms such as MicaSense RedEdge-MX, DJI P4 Multispectral, and Parrot Sequoia. Sentinel-2 L2A Drone Multispectral (MicaSense / DJI P4 / Sequoia)

RESOLUTION 10–60 m per band (fixed by ESA) RESOLUTION 1–5 cm orthomosaic (flight-dependent)

BAND COUNT 13 bands (B02–B8A, B11, B12…) BAND COUNT 4–6 bands (Blue, Green, Red, RE, NIR)

FILE FORMAT .jp2 or .tif per band + XML manifest FILE FORMAT Stacked .tif or per-band GeoTIFFs

RADIOMETRIC SCALING (DN − BOA_OFFSET) / 10 000 Offset = −1000 for PB ≥ 04.00, else 0 RADIOMETRIC SCALING DN × panel factor or SCALE TIFF tag Stitcher may pre-apply calibration

BAND ORDER Fixed by ESA specification BAND ORDER Firmware/stitcher-dependent — never assume

Step-by-step implementation

Step 1 — Parse the Sentinel-2 L2A BOA offset from XML

Processing Baseline 04.00 (applied to scenes acquired on or after 25 January 2022) added a BOA_ADD_OFFSET of −1000 to every band. Older scenes carry an implicit offset of 0. Read the value from MTD_MSIL2A.xml rather than hardcoding it, so the same function works on both old and new tiles.

PYTHON 
import xml.etree.ElementTree as ET
from pathlib import Path

def read_s2_boa_offsets(safe_dir: Path) -> dict[str, int]:
    """
    Parse BOA_ADD_OFFSET values from a Sentinel-2 SAFE tile directory.
    Returns a dict mapping band_id (e.g. 'B02') to integer offset.
    Older products have no offset node; default is 0.
    """
    xml_path = next(safe_dir.glob("MTD_MSIL2A.xml"), None)
    if xml_path is None:
        raise FileNotFoundError(f"MTD_MSIL2A.xml not found in {safe_dir}")

    tree = ET.parse(xml_path)
    ns = {"n1": "https://psd-14.sentinel2.eo.esa.int/PSD/S2_PDI_Level-2A_Tile_Metadata.xsd"}
    offsets: dict[str, int] = {}

    # BOA_ADD_OFFSET elements appear under General_Info/Product_Image_Characteristics
    for elem in tree.iter("BOA_ADD_OFFSET"):
        band_id = elem.attrib.get("band_id", "")
        offsets[band_id] = int(elem.text or 0)

    return offsets  # empty dict for pre-PB04.00 tiles — treat as all-zero

Verify the result before proceeding:

PYTHON 
offsets = read_s2_boa_offsets(Path("/data/S2A_MSIL2A_20230815T103031_N0509_R108_T32UME_20230815T135912.SAFE"))
print(offsets)  # {'0': -1000, '1': -1000, ...} for PB04.00+; {} for older

Step 2 — Load and scale Sentinel-2 bands to reflectance

PYTHON 
import rasterio
import numpy as np
from pathlib import Path

def load_s2_band(jp2_path: Path, boa_offset: int = 0) -> np.ndarray:
    """
    Load a single Sentinel-2 L2A band and return float32 reflectance in [0, 1].
    boa_offset: integer from MTD_MSIL2A.xml (0 for pre-PB04.00 tiles, -1000 for newer).
    """
    with rasterio.open(jp2_path) as src:
        assert src.crs is not None, f"No CRS on {jp2_path} — check the .jp2 file integrity"
        dn = src.read(1).astype(np.float32)

    reflectance = (dn + boa_offset) / 10_000.0
    reflectance = np.clip(reflectance, 0.0, 1.0)

    assert reflectance.max() <= 1.0, "Reflectance ceiling exceeded after scaling"
    return reflectance

Step 3 — Parse drone multispectral bands by name, not index

Never assume a positional index for drone bands. MicaSense RedEdge-MX writes band descriptions into the TIFF; DJI P4 Multispectral and Parrot Sequoia may use file-naming conventions instead. The function below tries band descriptions first, then falls back to a user-supplied mapping.

PYTHON 
def load_drone_bands(
    tif_path: Path,
    band_name_map: dict[str, int] | None = None,
    scale: float | None = None,
) -> dict[str, np.ndarray]:
    """
    Load a stacked drone multispectral GeoTIFF and return a dict mapping
    band names (e.g. 'Red', 'NIR') to float32 reflectance arrays.

    band_name_map: explicit override {name: 1-based rasterio index}, e.g.
        {'Blue': 1, 'Green': 2, 'Red': 3, 'RedEdge': 4, 'NIR': 5}
    scale: if None, the function reads the SCALE TIFF tag or defaults to 1.0
           (assuming the stitcher already produced 0–1 reflectance).
    """
    with rasterio.open(tif_path) as src:
        assert src.crs is not None, f"No CRS embedded in {tif_path}"

        # Resolve band mapping
        if band_name_map is not None:
            mapping = band_name_map
        elif src.descriptions and all(src.descriptions):
            # rasterio returns a tuple of strings; indices are 0-based here
            mapping = {name: idx + 1 for idx, name in enumerate(src.descriptions)}
        else:
            raise ValueError(
                f"Cannot determine band order for {tif_path}. "
                "Pass an explicit band_name_map."
            )

        # Resolve scale factor
        if scale is None:
            tags = src.tags()
            scale = float(tags.get("SCALE", tags.get("scale", 1.0)))

        bands: dict[str, np.ndarray] = {}
        for name, rasterio_idx in mapping.items():
            dn = src.read(rasterio_idx).astype(np.float32)
            refl = np.clip(dn * scale, 0.0, 1.0)
            bands[name] = refl

    return bands


# Quick verification — median NIR reflectance over a healthy canopy should be > 0.3
result = load_drone_bands(
    Path("/flights/field_a/orthomosaic.tif"),
    band_name_map={"Blue": 1, "Green": 2, "Red": 3, "RedEdge": 4, "NIR": 5},
)
assert 0.0 < float(np.nanmedian(result["NIR"])) < 1.0, "NIR out of reflectance range"
print({k: round(float(np.nanmedian(v)), 3) for k, v in result.items()})

Step 4 — Build a unified reflectance dict from either platform

PYTHON 
def parse_multispectral(
    filepath: Path,
    platform: str,          # "sentinel2" or "drone"
    s2_safe_dir: Path | None = None,
    s2_band_id: str | None = None,
    drone_band_map: dict[str, int] | None = None,
    drone_scale: float | None = None,
) -> dict:
    """
    Unified entry point. Returns:
        {
          'reflectance': np.ndarray (float32),   # (bands, H, W) or single band
          'crs': rasterio.crs.CRS,
          'transform': Affine,
          'platform': str,
        }
    """
    if platform == "sentinel2":
        offsets = read_s2_boa_offsets(s2_safe_dir) if s2_safe_dir else {}
        offset = offsets.get(s2_band_id, 0)
        with rasterio.open(filepath) as src:
            crs, transform = src.crs, src.transform
        refl = load_s2_band(filepath, boa_offset=offset)
        return {"reflectance": refl[np.newaxis], "crs": crs, "transform": transform, "platform": "sentinel2"}

    elif platform == "drone":
        with rasterio.open(filepath) as src:
            crs, transform = src.crs, src.transform
        bands = load_drone_bands(filepath, band_name_map=drone_band_map, scale=drone_scale)
        stack = np.stack(list(bands.values()), axis=0)
        return {"reflectance": stack, "crs": crs, "transform": transform, "platform": "drone"}

    else:
        raise ValueError(f"Unknown platform '{platform}'. Use 'sentinel2' or 'drone'.")

Gotchas and edge cases

  • The BOA offset trap. Sentinel-2 scenes from before 25 January 2022 have no BOA_ADD_OFFSET element in MTD_MSIL2A.xml. Reading an absent node returns an empty dict; treat that as zero. Hardcoding −1000 for all scenes will push pre-2022 reflectance values down by 0.1, producing negative NDVI over bare soil.

  • rasterio band indices are 1-based; NumPy array indices are 0-based. src.read(1) returns band 1 as a 2-D array at array[0] if you call src.read() (all bands). Mixing these is the most common source of misassigned Red Edge values.

  • Drone stitchers sometimes pre-apply calibration — then the SCALE tag is 1.0. Always spot-check output by asserting 0.02 < median_reflectance < 0.65 over a vegetated area. Values consistently near 1.0 or below 0.0 indicate a double- or missing-calibration application.

  • Partial reprojection before parsing. If you reproject a JP2 Sentinel-2 band with rasterio.warp.reproject before scaling, the interpolated DN values may fall outside the original integer range, causing subtle saturation after the divide. Scale first, reproject second, or handle nodata explicitly when reprojecting float32 data. See understanding CRS in precision agriculture for the recommended order of operations.

FAQ

Why are my Sentinel-2 NDVI values stuck near 1.0?

You are computing NDVI on raw DN values (3000–8000 range) rather than on scaled reflectance. Apply (DN + offset) / 10000 before any band arithmetic. After scaling, canopy NIR reflectance should sit around 0.40–0.55 and Red around 0.04–0.12.

How do I know which band index is Red Edge on a drone orthomosaic?

Read src.descriptions from rasterio.open. If the stitching software wrote band names (common in MicaSense RedEdge-MX exports via Pix4D), those names appear as a tuple of strings. If the tuple contains None values, fall back to the sensor’s file-naming convention or supply band_name_map explicitly.

Can I stack Sentinel-2 and drone bands in the same array for comparison?

Only after spatial alignment. Sentinel-2 10 m pixels and drone 2–5 cm pixels differ by three orders of magnitude in area. Reproject and resample both datasets to a shared EPSG and pixel size using rasterio.warp.reproject before stacking them. Using EPSG:32632 (or the appropriate UTM zone) is recommended; see how to convert WGS84 to UTM for farm mapping for the reprojection pattern.

This guide is part of Ingesting Multispectral Drone Imagery — see there for the full pipeline context, including radiometric calibration with ground panels, spatial alignment across multiple drone sensor lenses, and memory-efficient band stacking for large orthomosaics.

Parent section Ingesting Multispectral Drone Imagery in Python