Masking Clouds Before Computing NDVI

TL;DR Answer

Never compute NDVI on raw surface reflectance. Clouds, cloud shadows, and cirrus produce reflectance values that map to plausible-looking NDVI numbers, so they slip past range checks and silently poison every downstream feature. Build a boolean mask from the Sentinel-2 Scene Classification Layer (SCL) or the Landsat QA_PIXEL bitmask, dilate it a few pixels to catch cloud edges, set the flagged pixels to NaN, and compute NDVI only on the clean array. This is a prerequisite for reliable NDVI and EVI calculation with rasterio.

import numpy as np

def mask_and_ndvi(red, nir, mask):
    """Compute NDVI with cloud/shadow pixels set to NaN."""
    red = np.where(mask, np.nan, red.astype("float32"))
    nir = np.where(mask, np.nan, nir.astype("float32"))
    denom = nir + red
    ndvi = np.divide(nir - red, denom, out=np.full_like(red, np.nan),
                     where=denom != 0)
    return ndvi

Part of: Raster Band Math and Index Calculation

Why Unmasked Clouds Wreck NDVI Features

The Normalized Difference Vegetation Index is (NIR - RED) / (NIR + RED). It is only meaningful when the red and near-infrared reflectance actually came from the land surface. A cloud reflects strongly and fairly flatly across both bands, so its NDVI collapses toward zero — indistinguishable from bare soil or water. A cloud shadow darkens both bands unevenly and can push NDVI toward either extreme. Thin cirrus and haze are the most dangerous because they only slightly perturb reflectance, yielding an NDVI that lands squarely inside the valid [-1, 1] range while being physically meaningless.

Because these corrupted values look numerically ordinary, they defeat the naive defence of clipping NDVI to a plausible band. You cannot tell a real NDVI of 0.15 over dry rangeland from a cirrus-contaminated 0.15 over a forest. The only reliable signal is the per-pixel scene classification that the data provider already computed from the full spectral signature and viewing geometry.

The damage compounds through the pipeline:

• Training labels drift. If NDVI thresholds define your land-cover labels, cloud-driven values mislabel vegetated pixels as bare or vice versa.
• Temporal features lie. A single cloudy date in a monthly composite can drag a mean or amplitude feature far from the true phenology curve. This is why cloud masking must happen before any temporal aggregation of time-series geodata.
• Distribution shifts silently. Contaminated pixels fatten the tails of your feature distributions, so scalers fit on dirty data and models see input ranges at inference that never appeared in training.

Masking first is cheaper and more honest than any amount of post-hoc cleanup.

The Two Mask Sources You Will Actually Use

Sentinel-2 Scene Classification Layer (SCL)

Sentinel-2 Level-2A products ship an SCL band at 20 metre resolution where each pixel holds a single integer class. The classes you almost always exclude are 3 (cloud shadow), 8 (cloud medium probability), 9 (cloud high probability), 10 (thin cirrus), and usually 11 (snow/ice) and 1 (saturated/defective). Building the mask is a simple membership test — no bit arithmetic is needed because SCL is already categorical.

CLOUD_SCL_CLASSES = {1, 3, 8, 9, 10, 11}

def scl_cloud_mask(scl: np.ndarray) -> np.ndarray:
    """Boolean mask: True where the pixel should be discarded."""
    return np.isin(scl, list(CLOUD_SCL_CLASSES))

Because SCL is 20 m and the red/NIR bands are 10 m, resample the mask to the target grid with nearest-neighbour so class codes are never averaged into meaningless intermediates. If your bands and SCL live on different grids entirely, align them first as covered in reprojecting a raster with rasterio warp reproject.

Landsat QA_PIXEL bitmask

Landsat Collection 2 encodes quality as a 16-bit packed integer in QA_PIXEL. Each bit is an independent flag, so you unpack it with bitwise operations rather than equality tests. The bits that matter for cloud masking:

Extract a single flag by shifting it down to bit zero and masking with 1:

def qa_pixel_mask(qa: np.ndarray) -> np.ndarray:
    """Decode Landsat QA_PIXEL into a discard mask."""
    dilated = (qa >> 1) & 1
    cirrus  = (qa >> 2) & 1
    cloud   = (qa >> 3) & 1
    shadow  = (qa >> 4) & 1
    return (dilated | cirrus | cloud | shadow).astype(bool)

The >> shift moves the target bit into the ones position and & 1 isolates it. Combining the four flags with bitwise OR yields a single boolean array that is True wherever any unwanted condition holds.

Dilating the Mask Around Cloud Edges

Automated cloud detectors are tuned to flag opaque cloud cores confidently, but they routinely under-call the feathered rim where a cloud thins into haze. That rim is exactly where contaminated-but-plausible NDVI values come from. A small morphological dilation grows the mask outward by a fixed number of pixels so the fringe is swept up with the core.

from scipy.ndimage import binary_dilation

def dilate_mask(mask: np.ndarray, iterations: int = 2) -> np.ndarray:
    """Grow a boolean cloud mask outward to cover cloud-edge haze."""
    structure = np.array([[0, 1, 0],
                          [1, 1, 1],
                          [0, 1, 0]], dtype=bool)
    return binary_dilation(mask, structure=structure, iterations=iterations)

Two iterations of a 4-connected structuring element grow the mask by roughly two pixels in each cardinal direction — about 20 m at Sentinel-2’s 10 m grid. Treat iterations as a tunable trade-off: a wider buffer discards more clear data but leaves less haze, which is the right call for dense time series where you can afford to drop observations. When clear looks are scarce, keep the dilation tight.

Production-Ready Masked NDVI Function

The function below ties the pieces together: it validates shapes, accepts a pre-built boolean mask from either sensor, optionally dilates it, propagates NaN through the division safely, and returns a float32 array with a matching updated nodata convention. It never raises on a divide-by-zero because the where argument guards the NIR + RED == 0 case.

from __future__ import annotations

import logging

import numpy as np
from scipy.ndimage import binary_dilation

logger = logging.getLogger(__name__)


def mask_and_ndvi(
    red: np.ndarray,
    nir: np.ndarray,
    cloud_mask: np.ndarray,
    dilate_iterations: int = 2,
) -> np.ndarray:
    """Compute NDVI with cloud/shadow pixels propagated to NaN.

    Parameters
    ----------
    red, nir : np.ndarray
        Surface reflectance bands, identical 2-D shape. Any integer or
        float dtype; scaled or unscaled reflectance both work because
        NDVI is a ratio.
    cloud_mask : np.ndarray of bool
        True where a pixel should be discarded (cloud, shadow, cirrus).
        Same shape as red/nir.
    dilate_iterations : int, default=2
        Morphological dilation applied to the mask to catch cloud edges.
        Pass 0 to skip dilation.

    Returns
    -------
    np.ndarray of float32
        NDVI in [-1, 1] on clear pixels, NaN everywhere masked.
    """
    if not (red.shape == nir.shape == cloud_mask.shape):
        raise ValueError(
            f"Shape mismatch: red={red.shape}, nir={nir.shape}, "
            f"mask={cloud_mask.shape}"
        )
    if cloud_mask.dtype != bool:
        raise TypeError("cloud_mask must be a boolean array.")

    mask = cloud_mask
    if dilate_iterations > 0:
        structure = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=bool)
        mask = binary_dilation(mask, structure=structure,
                               iterations=dilate_iterations)

    red_f = np.where(mask, np.nan, red.astype("float32"))
    nir_f = np.where(mask, np.nan, nir.astype("float32"))

    denom = nir_f + red_f
    ndvi = np.divide(
        nir_f - red_f,
        denom,
        out=np.full(red_f.shape, np.nan, dtype="float32"),
        where=denom != 0,
    )

    discarded = int(mask.sum())
    logger.info(
        "Masked %d/%d pixels (%.1f%%) before NDVI.",
        discarded, mask.size, 100 * discarded / mask.size,
    )
    return ndvi

Note that np.where(mask, np.nan, ...) forces the output to floating point, which is what lets NaN survive — you cannot store NaN in an integer array. Any pixel that was cloudy, or whose red and NIR summed to zero, comes back as NaN, giving you one unambiguous nodata sentinel to carry downstream.

Walkthrough on a Sentinel-2 Tile

This end-to-end example reads a Level-2A tile with rasterio, resamples the 20 m SCL band up to the 10 m red/NIR grid, and produces a clean NDVI array. The diagram after it shows the full data flow.

import numpy as np
import rasterio
from rasterio.enums import Resampling

with rasterio.open("T33UUP_B04_10m.jp2") as src:   # red
    red = src.read(1)
    profile = src.profile
    dst_shape = (src.height, src.width)

with rasterio.open("T33UUP_B08_10m.jp2") as src:   # NIR
    nir = src.read(1)

# SCL is 20 m; read and resample to the 10 m grid with nearest-neighbour
with rasterio.open("T33UUP_SCL_20m.jp2") as src:
    scl = src.read(1, out_shape=dst_shape, resampling=Resampling.nearest)

cloud_mask = np.isin(scl, [1, 3, 8, 9, 10, 11])
ndvi = mask_and_ndvi(red, nir, cloud_mask, dilate_iterations=2)

# Persist NDVI as float32 with NaN nodata so masking survives on disk
profile.update(dtype="float32", count=1, nodata=np.nan)
with rasterio.open("T33UUP_NDVI_masked.tif", "w", **profile) as dst:
    dst.write(ndvi.astype("float32"), 1)

The saved GeoTIFF now carries nodata=nan, so every consumer — a zonal-statistics pass, a temporal composite, or a model’s inference reader — treats the cloudy pixels as missing rather than as real low-vegetation observations.

Verification

Two assertions catch the failure modes that matter: NDVI must never fall outside [-1, 1] on the clear pixels, and every pixel the mask flagged must be NaN in the output. Run these on every tile as a cheap regression guard.

def verify_masked_ndvi(ndvi: np.ndarray, cloud_mask: np.ndarray,
                       dilate_iterations: int = 2) -> None:
    """Assert masked NDVI is physically valid and cloud pixels are NaN."""
    from scipy.ndimage import binary_dilation

    structure = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=bool)
    mask = binary_dilation(cloud_mask, structure=structure,
                           iterations=dilate_iterations) \
        if dilate_iterations else cloud_mask

    # 1. Every dilated cloud pixel must be NaN
    assert np.isnan(ndvi[mask]).all(), "Cloud pixels leaked into NDVI output."

    # 2. Clear pixels must be in-range (NaN comparisons are ignored)
    clear = ndvi[~mask]
    finite = clear[np.isfinite(clear)]
    assert np.all((finite >= -1.0) & (finite <= 1.0)), \
        "NDVI out of [-1, 1] on clear pixels — check band scaling."

    print(f"OK: {np.isnan(ndvi).sum()} NaN pixels, "
          f"clear NDVI range [{finite.min():.3f}, {finite.max():.3f}]")


verify_masked_ndvi(ndvi, cloud_mask, dilate_iterations=2)

An out-of-range assertion almost always means the red and NIR bands were on different scale factors — a common trap when one band is raw digital number and the other is already scaled reflectance. Because NDVI is a ratio, consistent scaling matters more than the absolute values; if you also feed these bands to a model, handle their distributions with robust scaling for skewed spectral bands.

FAQ

Why not just filter out extreme NDVI values after computing the index?

Clipping NDVI to a plausible range hides the symptom but keeps the disease. A thin cirrus or haze pixel produces an NDVI that sits comfortably inside [-1, 1] yet is physically wrong — the reflectance was never vegetation. Threshold filtering cannot separate a genuinely low-NDVI bare-soil pixel from a cloud-contaminated one. The scene classification and QA bands were computed from the full spectral signature and geometry, so masking on them removes the actual cause of corruption before any arithmetic runs.

How do I decode the Landsat QA_PIXEL band into a cloud mask?

QA_PIXEL is a 16-bit packed bitmask where each bit flags a condition — bit 3 is cloud, bit 4 is cloud shadow, bit 1 is dilated cloud, bit 2 is cirrus. Extract a flag with a bitwise AND against a shifted mask: (qa >> 3) & 1 gives the cloud bit for every pixel. Combine the bits you care about with logical OR to build a single boolean array, then invert it to keep the clear pixels.

Should I dilate the cloud mask, and by how much?

Yes. Cloud detection routinely misses the semi-transparent fringe one to three pixels beyond the opaque core, and that fringe carries the most damaging haze. Dilating the mask by a small binary structuring element — a 3x3 to 5x5 kernel, roughly 20 to 50 metres at 10 m resolution — buys a safety margin at the cost of discarding a few clear pixels. Widen it for dense time series where clean coverage is plentiful, and narrow it when observations are scarce.

Related

• How to Calculate NDVI and EVI with Rasterio
• Reprojecting a Raster with rasterio warp reproject
• Aggregating Daily Satellite Data to Monthly Features
• Temporal Aggregation for Time-Series Geodata
• Applying Robust Scaling to Skewed Spectral Bands

Part of: Raster Band Math and Index Calculation — Spatial Feature Engineering for Machine Learning