Masking Clouds with the Sentinel-2 Scene Classification Layer in Python

An unmasked cloud is not noise — it is a bright, systematic, seasonally clustered bias, and in a forest NDVI series it looks exactly like a disturbance. Sentinel-2 Level-2A products ship a Scene Classification Layer (SCL) that labels every 20 m pixel as vegetation, cloud, shadow, snow and so on, and it is the cheapest defensible mask available. Using it well is a matter of choosing classes deliberately, dilating the result, and knowing what it misses. This guide covers that, inside Vegetation Index Calculation in Python, part of Ecological GIS Data Foundations in Python.

When to use this approach

The practical recipe for forest work is SCL plus a dilation, followed by a temporal median that absorbs whatever leaked through. Scene-level cloud cover belongs in the STAC query as a loose pre-filter only — as covered in Searching STAC Catalogs for Sentinel-2 with pystac-client.

The classes, and which to keep

SCL is a single band of small integers. The choice of which values count as usable is the whole method.

Keeping class 5 matters more than it looks. A clearcut, a fresh burn, or a heavily thinned stand classifies as “not vegetated”, and dropping it removes exactly the pixels a disturbance analysis exists to find — producing a time series in which harvest looks like missing data.

Minimal reproducible example

The mask itself is one isin. Everything that makes it work is the dilation around it — SCL boundaries are conservative, and a cloud’s real influence extends a pixel or two beyond where the classifier draws its edge.

import numpy as np
import rasterio
from scipy.ndimage import binary_dilation

CLEAR_CLASSES = (4, 5, 6, 7)          # vegetation, bare, water, unclassified
CLOUDY_CLASSES = (3, 8, 9, 10, 11)    # shadow, cloud, cirrus, snow


def clear_mask(scl: np.ndarray, dilate_px: int = 3,
               keep_water: bool = True) -> np.ndarray:
    """Boolean mask of usable pixels from a Sentinel-2 SCL band."""
    keep = set(CLEAR_CLASSES)
    if not keep_water:
        keep.discard(6)
    clear = np.isin(scl, list(keep))

    if dilate_px:
        # Grow the cloud/shadow footprint outward; SCL edges are optimistic.
        bad = np.isin(scl, CLOUDY_CLASSES)
        structure = np.ones((3, 3), dtype=bool)
        bad = binary_dilation(bad, structure=structure, iterations=dilate_px)
        clear &= ~bad
    return clear


def masked_ndvi(red_path: str, nir_path: str, scl_path: str,
                dilate_px: int = 3) -> tuple[np.ndarray, dict]:
    """NDVI with cloud, shadow and cirrus set to NaN."""
    with rasterio.open(red_path) as r, rasterio.open(nir_path) as n, \
            rasterio.open(scl_path) as s:
        red = r.read(1).astype("float32")
        nir = n.read(1).astype("float32")
        scl = s.read(1, out_shape=red.shape,
                     resampling=rasterio.enums.Resampling.nearest)
        profile = r.profile

    clear = clear_mask(scl, dilate_px=dilate_px)
    denom = nir + red
    ndvi = np.where(denom != 0, (nir - red) / denom, np.nan).astype("float32")
    ndvi[~clear] = np.nan
    profile.update(dtype="float32", nodata=np.nan, count=1)
    return ndvi, profile

Two details are load-bearing. SCL is delivered at 20 m while red and NIR are 10 m, so it must be resampled — and with nearest resampling, because averaging class codes 4 and 9 produces 6, which is water. And the guard on denom != 0 prevents a division that would otherwise emit warnings and infinities across every no-data pixel.

Parameter reference

Expected output and verification

The number that decides whether a composite is trustworthy is not the mask fraction on any one date — it is the count of clear observations per pixel across the series. A median over three observations is a very different statistic from a median over thirty.

import numpy as np


def clear_observation_count(scl_stack: np.ndarray, dilate_px: int = 3) -> np.ndarray:
    """Per-pixel count of clear dates from a stacked SCL array (time, y, x)."""
    counts = np.zeros(scl_stack.shape[1:], dtype="int16")
    for t in range(scl_stack.shape[0]):
        counts += clear_mask(scl_stack[t], dilate_px=dilate_px).astype("int16")
    return counts


counts = clear_observation_count(scl_stack)
print({
    "dates": int(scl_stack.shape[0]),
    "min_clear": int(counts.min()),
    "median_clear": int(np.median(counts)),
    "pixels_under_5": int((counts < 5).sum()),
})
assert np.median(counts) >= 5, "fewer than 5 clear dates at the median pixel — widen the window"

Ship that count raster alongside every composite. It is the honest uncertainty layer: a pixel with two clear observations and a pixel with thirty carry the same NDVI value and very different confidence, and without the count nothing downstream can tell them apart.

The second check is that masking actually changed the distribution in the direction it should. Clouds are bright in red and dark in NIR relative to vegetation, so removing them should raise the NDVI mean and shrink its spread:

def mask_effect(ndvi_raw: np.ndarray, ndvi_masked: np.ndarray) -> dict:
    """Confirm the mask moved the distribution the way clouds imply."""
    raw = ndvi_raw[np.isfinite(ndvi_raw)]
    msk = ndvi_masked[np.isfinite(ndvi_masked)]
    return {
        "kept_frac": round(float(msk.size / raw.size), 3),
        "mean_raw": round(float(raw.mean()), 3),
        "mean_masked": round(float(msk.mean()), 3),
        "sd_raw": round(float(raw.std()), 3),
        "sd_masked": round(float(msk.std()), 3),
    }


eff = mask_effect(ndvi_raw, ndvi_masked)
assert eff["mean_masked"] >= eff["mean_raw"], f"masking lowered mean NDVI: {eff}"
assert eff["sd_masked"] <= eff["sd_raw"], f"masking widened the spread: {eff}"

If masking lowers the mean, the mask is inverted — a mistake that is otherwise invisible because the output still looks like a plausible NDVI raster.

Common pitfalls

• Resampling SCL bilinearly. Interpolating class codes invents classes. Always nearest.
• Dropping class 5. Harvested and burned stands classify as not-vegetated, so excluding it deletes the disturbances a monitoring series exists to detect.
• No dilation. SCL cloud edges are tight, and the pixels immediately outside them are contaminated by adjacency scattering. Three pixels of dilation is cheap insurance.
• Trusting the mask over snow. Snow and bright cloud are routinely confused in both directions. In winter series, exclude both and expect a genuine seasonal data gap.
• Compositing without reporting the clear count. A median over two observations is not comparable with a median over twenty, and nothing in the output raster records which one you have.

Frequently Asked Questions

Related

• Vegetation Index Calculation in Python — the wider index workflow this mask protects
• Calculating NDVI from Sentinel-2 with rasterio — the index computation this feeds
• Chunked Raster Processing with xarray and dask — applying this mask across a lazy stack
• Measuring Post-Fire NDVI Recovery Trajectories in Python — an analysis that an unmasked series would ruin

Up: Vegetation Index Calculation · Ecological GIS Data Foundations in Python