Cloud Masking for Agricultural Imagery

Cloud contamination is one of the most persistent bottlenecks in precision agriculture workflows. When optical sensors capture crop canopies, atmospheric interference introduces spectral noise that corrupts vegetation indices, skews yield models, and compromises field-level decision support. This guide walks through a production-ready cloud masking pipeline that takes a raw Sentinel-2 L2A scene or multispectral drone orthomosaic as input and produces a validated Cloud Optimized GeoTIFF (COG) binary mask — every pixel flagged 1 for cloud or shadow, 0 for valid surface reflectance — ready to feed directly into band math and raster algebra workflows.

This page is part of the Drone Imagery Processing & Vegetation Index Workflows guide. See that page for the full pipeline context from ingest through temporal aggregation.

Step 1 Step 2 Step 3 Step 4 Step 5

Prerequisites

Python packages (exact versions tested):

• rasterio>=1.3.9
• numpy>=1.24
• scipy>=1.10
• scikit-image>=0.21
• GDAL>=3.6 compiled with Cloud Optimized GeoTIFF support

Input data requirements:

• Multi-band GeoTIFF in a projected CRS (e.g. EPSG:32632 for UTM zone 32N). Geographic CRS (EPSG:4326) will produce incorrect pixel-size calculations in the morphological stages — reproject first using the approach described in understanding CRS in precision agriculture.
• Band order: blue (B02), green (B03), red (B04), NIR (B08), SWIR1 (B11) for Sentinel-2 L2A; or equivalent mapped bands from a MicaSense RedEdge-MX or DJI P4 Multispectral flight.
• Radiometric calibration complete — raw digital numbers will produce inconsistent thresholds across acquisition dates. For drone imagery, consult the ingesting multispectral drone imagery guide for calibration steps before running this pipeline.
• Nodata value set in the GeoTIFF profile (commonly 0 or 65535 for uint16 Sentinel-2 products).
• Memory: agricultural rasters frequently exceed 10 GB. The windowed approach below keeps RAM usage under 2 GB for any scene size.

1. Concept & Algorithm

Effective cloud masking in agricultural imagery requires more than a single brightness threshold. Opaque cumulus clouds, thin cirrus layers, cloud shadows, and high-reflectance surfaces — dry bare soil, plastic mulch, greenhouse roofs, standing water — all occupy overlapping spectral regions. Misidentifying any of these as valid canopy injects false values into the NDVI and NDRE time series used for crop stress detection and variable-rate prescription generation.

The algorithm here follows a simplified Fmask logic adapted for agricultural landscapes:

Spectral stage. Cloud pixels have elevated blue-band reflectance (>0.20), elevated NIR (>0.15), and elevated SWIR1 (>0.10). The combination of all three separates clouds from most agricultural surfaces. The Normalized Difference Snow Index (NDSI = (Green − SWIR1) / (Green + SWIR1)) is computed to prevent snow-covered fields in winter acquisitions from being misclassified as cloud; NDSI > 0.4 indicates snow, not cloud.

Shadow stage. Cloud shadows appear as pixels with low reflectance across all visible and NIR bands (typically below 0.08). Because shadows can be mistaken for water bodies, flooded paddies, or saturated dark soils — each of which carries genuine agronomic information — a geometric displacement step projects the shadow zone along the solar azimuth vector before finalising the shadow class.

Morphological stage. Isolated bright speckles (specular water reflections on irrigation canals, white-painted structures) and isolated dark pixels (tractor shadows, access road patches) survive the spectral stage. Connected-component filtering with a minimum object size threshold removes these false positives before the mask is written to disk.

Why does this matter agronomically? A single un-masked cloud or shadow pixel over a crop field can shift the NDVI at that location by 0.15–0.40 reflectance units — far exceeding the 0.05–0.10 change that triggers a management zone reclassification in a variable-rate application map. Consistent masking is therefore a non-negotiable prerequisite for repeatable prescription generation, as detailed in threshold mapping for crop health detection.

2. Step-by-step Implementation

Step 1 — Inspect band metadata and confirm CRS

import rasterio
import numpy as np

INPUT_PATH = "sentinel2_l2a_T32UNU_20240715.tif"

with rasterio.open(INPUT_PATH) as src:
    assert src.crs is not None, "CRS missing — reproject before masking"
    assert src.crs.is_projected, (
        f"Expected projected CRS, got {src.crs.to_epsg()}. "
        "Reproject to UTM before running this pipeline."
    )
    print(f"CRS: EPSG:{src.crs.to_epsg()}")
    print(f"Bands: {src.count}, Dtype: {src.dtypes[0]}")
    print(f"Dimensions: {src.width} x {src.height} px")
    print(f"Resolution: {src.res[0]:.2f} m/px")
    print(f"Nodata: {src.nodata}")
    # Confirm enough bands exist for spectral masking
    assert src.count >= 5, "Need at least 5 bands: B02, B03, B04, B08, B11"

Step 2 — Define reflectance scaling and spectral mask functions

Sentinel-2 L2A products store surface reflectance as uint16 values scaled by 10 000. Landsat Collection 2 Level-2 uses a different factor (0.0000275) with an additive offset of −0.2. The functions below accept a scale_factor and additive_offset argument so they work for either sensor.

def scale_to_reflectance(
    band: np.ndarray,
    scale_factor: float = 0.0001,
    additive_offset: float = 0.0,
    nodata: float = 0.0
) -> np.ndarray:
    """Convert raw DN to surface reflectance [0, 1]. Nodata pixels become NaN."""
    out = band.astype("float32") * scale_factor + additive_offset
    out = np.clip(out, 0.0, 1.0)
    out[band == nodata] = np.nan
    return out


def compute_spectral_mask(
    blue: np.ndarray,
    green: np.ndarray,
    nir: np.ndarray,
    swir1: np.ndarray,
    blue_thresh: float = 0.20,
    nir_thresh: float = 0.15,
    swir1_thresh: float = 0.10,
    shadow_thresh: float = 0.08,
    ndsi_snow_thresh: float = 0.40,
) -> np.ndarray:
    """
    Return a boolean mask: True = cloud or shadow, False = clear.
    NaN pixels (nodata) are always masked True.
    """
    nodata_mask = np.isnan(blue) | np.isnan(nir) | np.isnan(swir1)

    # Snow guard: high NDSI pixels are snow, not cloud
    with np.errstate(invalid="ignore", divide="ignore"):
        ndsi = np.where(
            (green + swir1) > 0,
            (green - swir1) / (green + swir1),
            0.0,
        )
    is_snow = ndsi > ndsi_snow_thresh

    # Cloud: bright blue + elevated NIR + elevated SWIR1, not snow
    cloud = (
        (blue > blue_thresh)
        & (nir > nir_thresh)
        & (swir1 > swir1_thresh)
        & ~is_snow
    )

    # Shadow: uniformly dark across visible + NIR
    shadow = (blue < shadow_thresh) & (nir < shadow_thresh)

    return cloud | shadow | nodata_mask

Step 3 — Morphological refinement

from scipy.ndimage import binary_dilation, generate_binary_structure
from skimage.morphology import remove_small_objects


def refine_mask(
    raw_mask: np.ndarray,
    min_cloud_size_px: int = 500,
    dilation_radius_px: int = 3,
) -> np.ndarray:
    """
    Remove isolated noise and buffer cloud/shadow edges.

    Parameters
    ----------
    raw_mask : bool array, True = cloud/shadow
    min_cloud_size_px : minimum connected-component size to keep
    dilation_radius_px : buffer added around every cloud/shadow pixel
    """
    # Remove tiny isolated detections (speckle, roof highlights)
    cleaned = remove_small_objects(raw_mask.copy(), min_size=min_cloud_size_px)

    # Buffer edges to capture cloud fringe and shadow penumbra
    struct = generate_binary_structure(2, 1)
    for _ in range(dilation_radius_px):
        cleaned = binary_dilation(cleaned, structure=struct)

    return cleaned

Step 4 — Windowed processing and COG output

from rasterio.windows import Window
from rasterio.enums import Resampling
import os


BAND_INDEX_MAP = {
    # Sentinel-2 L2A band order in a standard 10/20 m stack
    "blue": 1,   # B02
    "green": 2,  # B03
    "nir": 4,    # B08
    "swir1": 5,  # B11 (resampled to 10 m)
}


def build_cloud_mask_cog(
    input_path: str,
    output_path: str,
    chunk_size: int = 1024,
    min_cloud_size_px: int = 500,
    dilation_radius_px: int = 3,
    scale_factor: float = 0.0001,
    nodata_dn: float = 0.0,
    band_map: dict = None,
):
    """
    Windowed cloud + shadow masking pipeline. Writes a tiled COG (uint8).
    1 = cloud or shadow, 0 = clear surface.
    """
    if band_map is None:
        band_map = BAND_INDEX_MAP

    with rasterio.open(input_path) as src:
        assert src.crs is not None and src.crs.is_projected, (
            "Input must be in a projected CRS."
        )
        profile = src.profile.copy()
        profile.update(
            dtype="uint8",
            count=1,
            compress="lzw",
            predictor=2,
            tiled=True,
            blockxsize=256,
            blockysize=256,
            nodata=255,  # 255 = unprocessed / nodata
        )

        with rasterio.open(output_path, "w", **profile) as dst:
            for row_off in range(0, src.height, chunk_size):
                for col_off in range(0, src.width, chunk_size):
                    win = Window(
                        col_off,
                        row_off,
                        min(chunk_size, src.width - col_off),
                        min(chunk_size, src.height - row_off),
                    )

                    # Read only the bands required for masking
                    blue_dn  = src.read(band_map["blue"],  window=win)
                    green_dn = src.read(band_map["green"], window=win)
                    nir_dn   = src.read(band_map["nir"],   window=win)
                    swir1_dn = src.read(band_map["swir1"], window=win)

                    blue  = scale_to_reflectance(blue_dn,  scale_factor, nodata=nodata_dn)
                    green = scale_to_reflectance(green_dn, scale_factor, nodata=nodata_dn)
                    nir   = scale_to_reflectance(nir_dn,   scale_factor, nodata=nodata_dn)
                    swir1 = scale_to_reflectance(swir1_dn, scale_factor, nodata=nodata_dn)

                    raw = compute_spectral_mask(blue, green, nir, swir1)
                    refined = refine_mask(raw, min_cloud_size_px, dilation_radius_px)

                    out = refined.astype("uint8")
                    dst.write(out, 1, window=win)

    print(f"Mask written: {output_path}")
    return output_path


if __name__ == "__main__":
    build_cloud_mask_cog(
        input_path="sentinel2_l2a_T32UNU_20240715.tif",
        output_path="cloud_mask_T32UNU_20240715.tif",
    )

3. Key Parameters & Tuning

4. Handling Edge Cases & Failure Modes

Bare soil and plastic mulch flagged as cloud. Light-coloured soils (Mollisols, Aridisols) and UV-stabilised white plastic mulch used in vegetable production have blue reflectance between 0.18 and 0.28 — overlapping with the cloud threshold. Add a SWIR2 (B12) ratio check: (swir1 - swir2) / (swir1 + swir2) < 0.1 is characteristic of cloud; soil ratios typically exceed 0.15. If a B12 band is unavailable, clip known agricultural polygon extents loaded from field boundary GeoPackages extracted with GeoPandas and relax thresholds only within those areas.

Shadow misclassified as crop stress. Cloud shadows in NIR-dark areas (stressed or senescing crops) share reflectance signatures. Because a variable-rate prescription that excludes shadow-masked pixels would incorrectly skip treating those areas, always validate shadow masks against the acquisition solar geometry. Compute the solar zenith and azimuth from the scene metadata and apply a directional dilation kernel along the shadow projection vector rather than an isotropic buffer. Sentinel-2 metadata provides MEAN_SOLAR_AZIMUTH_ANGLE and MEAN_SOLAR_ZENITH_ANGLE in the .SAFE directory.

UTM zone boundary crossings. A farm straddling two UTM zones (e.g., spanning 6°E for zone 32N/33N boundary in Germany) may be delivered in two separate tiles. Running the mask on mismatched CRS tiles and attempting to merge produces a seam artifact. Always confirm src.crs.to_epsg() matches between tiles before stitching. Use the CRS validation workflow in validating coordinate systems for variable-rate maps.

Zero-size windows at raster edges. When src.width or src.height is not evenly divisible by chunk_size, the final row or column window will have zero extent if not clipped. The min(chunk_size, src.width - col_off) guard in the implementation above prevents this, but verify it explicitly:

# Assertion to catch edge window bugs during development
for row_off in range(0, src.height, chunk_size):
    for col_off in range(0, src.width, chunk_size):
        w = min(chunk_size, src.width - col_off)
        h = min(chunk_size, src.height - row_off)
        assert w > 0 and h > 0, f"Zero-extent window at col={col_off}, row={row_off}"

DJI P4 Multispectral and MicaSense RedEdge-MX RGB-only fallback. When only three bands (RGB) are available from a Parrot Sequoia or lower-spec payload, the SWIR-based cloud check is unavailable. Use HSV brightness thresholding (V > 0.85 in the V channel from skimage.color.rgb2hsv) combined with a texture metric (local variance on a 5×5 window) to separate smooth cloud from structured canopy. Expect a 10–15% higher false-positive rate compared to a five-band pipeline.

5. Verification & Output Validation

Run the following checks immediately after build_cloud_mask_cog() completes:

import numpy as np
import rasterio

def validate_cloud_mask(mask_path: str, original_path: str) -> dict:
    """
    Confirm the cloud mask is geometrically and statistically coherent.
    Returns a dict of QA metrics.
    """
    with rasterio.open(mask_path) as msk, rasterio.open(original_path) as src:
        # 1. Spatial alignment
        assert msk.crs == src.crs, "CRS mismatch between mask and source"
        assert msk.transform == src.transform, "Transform mismatch — pixel grids not aligned"
        assert (msk.width, msk.height) == (src.width, src.height), "Dimension mismatch"

        # 2. Value range
        data = msk.read(1)
        unique_vals = set(np.unique(data))
        assert unique_vals.issubset({0, 1, 255}), (
            f"Unexpected values in mask: {unique_vals - {0, 1, 255}}"
        )

        # 3. Cloud fraction — flag if suspiciously high or low
        valid = data != 255
        cloud_frac = data[valid].mean()  # fraction of valid pixels that are cloud/shadow
        if cloud_frac > 0.80:
            print(f"WARNING: {cloud_frac:.1%} cloud fraction — thresholds may be too aggressive")
        if cloud_frac < 0.001:
            print(f"WARNING: {cloud_frac:.3%} cloud fraction — masking may have failed silently")

        return {
            "crs_ok": True,
            "transform_ok": True,
            "cloud_fraction": float(cloud_frac),
            "nodata_fraction": float((data == 255).mean()),
        }


report = validate_cloud_mask(
    "cloud_mask_T32UNU_20240715.tif",
    "sentinel2_l2a_T32UNU_20240715.tif",
)
print(report)

Visual spot-check. Overlay the mask as a semi-transparent red layer on the true-colour RGB composite in QGIS (Layer > Add Layer > Add Raster Layer; use Band Rendering → Singleband pseudocolour with a red ramp on value 1). Check that mask boundaries align with visible cloud edges and that crop field interiors with high NIR are not masked. Pay particular attention to:

• Irrigated pivot circles (high NIR reflectance can occasionally trigger cloud detection near the pivot centre — check nir_thresh)
• Riparian corridors with tall tree canopy (dense shadow can fall across adjacent fields — check dilation_radius_px)
• Field calibration tarpaulins (known-value reflectance panels used during MicaSense RedEdge-MX and DJI P4 Multispectral flights should always appear as masked=0)

6. Integration with the Broader Pipeline

The COG mask produced by this pipeline plugs directly into downstream analysis:

Upstream dependency. This pipeline assumes that ingesting multispectral drone imagery (including radiometric calibration and CRS alignment) has already been completed. Passing uncalibrated digital numbers will cause threshold failures regardless of parameter tuning.

Applying the mask before index calculation. When computing NDVI, NDRE, or SAVI using the band math and raster algebra patterns, load the mask in the same windowed loop and set np.nan on masked pixels before the ratio calculation. This prevents zero-denominator errors propagating from cloud-contaminated reflectance values:

# Inside a windowed loop, after reading red and nir arrays:
mask_window = cloud_mask_src.read(1, window=win)
red = red.astype("float32")
nir = nir.astype("float32")
red[mask_window == 1] = np.nan
nir[mask_window == 1] = np.nan

with np.errstate(invalid="ignore", divide="ignore"):
    ndvi = (nir - red) / (nir + red)

Feeding temporal aggregation. Once per-scene masks exist for a growing season, temporal aggregation of vegetation indices can composite clean pixels across dates, filling cloud-contaminated acquisitions with values from adjacent clear dates. The mask COG produced here is the direct input to that compositing workflow.

Batch automation. For monitoring programmes that ingest a new Sentinel-2 tile every 5 days, wrapping build_cloud_mask_cog() in an automated scheduler that handles metadata parsing, solar angle retrieval, and threshold selection per season is covered in automating cloud removal in Sentinel-2 time series.

Related

• Automating Cloud Removal in Sentinel-2 Time Series — batch orchestration, solar geometry computation, and gap-filling across seasonal stacks
• Band Math & Raster Algebra in Python — applying the clean reflectance stack to compute NDVI, NDRE, and SAVI
• Temporal Aggregation of Vegetation Indices — cloud-free compositing across multi-date acquisition stacks
• Threshold Mapping for Crop Health Detection — converting masked index rasters into management zone boundaries
• Ingesting Multispectral Drone Imagery — radiometric calibration and CRS alignment as the upstream prerequisite