Masking Clouds with the Sentinel-2 SCL Band

The Sentinel-2 Level-2A product ships a Scene Classification (SCL) band that already labels each pixel as vegetation, water, cloud, shadow, cirrus, and so on. To mask clouds and shadows, read the SCL band, select the contaminated class codes, and set matching reflectance pixels to NaN:

import numpy as np
import rasterio

CONTAMINATED = {3, 8, 9, 10, 11}  # shadow, cloud med/high, cirrus, snow

with rasterio.open("SCL_20m.tif") as scl_src:
    scl = scl_src.read(1)

mask = np.isin(scl, list(CONTAMINATED))  # True where contaminated

This is the fast, label-driven path inside Cloud and Shadow Masking Strategies. It is distinct from Implementing FMask and s2cloudless in Python: those are machine-learning cloud detectors you run yourself, whereas SCL is a categorical bitmask the ESA processor already computed. When the L2A product exists, SCL is the zero-cost first line of defence.

Why This Arises in Remote Sensing Workflows

Any multi-temporal Sentinel-2 analysis — vegetation trends, water extent, change detection — collapses the moment cloudy pixels enter the arithmetic. A cloud reads as a bright, spectrally flat surface, so an NDVI computed over one silently reports bare soil, and a median composite pulls toward whichever dates happened to be clear. Masking has to happen before any index or aggregation.

The SCL band exists precisely so you do not need to run a detector on every scene. During Level-2A processing, the ESA Sen2Cor processor classifies every 20 m pixel into one of twelve classes — vegetation, bare soil, water, cloud at two probability levels, cirrus, shadow, snow, and a few housekeeping categories — and stores the result alongside the reflectance bands. Reading it and selecting a handful of codes gives you a usable mask in milliseconds, with no model weights to download and no per-pixel inference to run.

The two refinements that separate a naive SCL mask from a production one are dilation and resolution alignment. Dilation grows the mask outward to catch the semi-transparent fringe the classifier under-flags: a cloud does not have a crisp edge, and the reflectance in the halo just outside a flagged cloud is still contaminated even though SCL labelled it clear. Resolution alignment matters because SCL is delivered at 20 m while the visible and near-infrared bands you are masking are at 10 m, so the mask has to be upsampled — and it must be upsampled with nearest-neighbour, because the class codes are categorical labels that averaging would turn into nonsense.

This page sits in Cloud and Shadow Masking Strategies under Satellite Processing Workflows & Index Pipelines. A clean SCL mask is the input that feeds Spectral Index Calculation Pipelines: mask first, then compute indices only over valid ground.

SCL Classes and the Masking Decision

The diagram shows the SCL class codes and which ones a typical cloud-and-shadow mask removes.

Class 2 (dark/topographic shadow) is often kept because masking it removes valid terrain-shaded ground; class 11 (snow) is masked only when snow is not your subject. The set is a study decision, not a fixed rule. It is worth being deliberate about the borderline classes: SCL’s cloud-shadow class (3) is geometrically projected from detected clouds and can miss shadows whose parent cloud was itself missed, while the medium-probability cloud class (8) sometimes over-flags bright bare soil or rooftops. For a conservative composite you accept a little over-masking and keep class 8 in the contaminated set; for maximum data retention over an arid site you might treat class 8 more carefully. The registry-style approach of defining the contaminated set as a named constant, rather than scattering literal codes through the code, makes that policy explicit and easy to audit.

Environment & Setup

pip install "rasterio>=1.3.0" "numpy>=1.24" "scipy>=1.10"

Complete Working Example

The function reads SCL, builds and dilates the contamination mask, resamples it from 20 m to the 10 m reflectance grid with nearest-neighbour, and applies it to a stack of reflectance bands.

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

# SCL codes considered contaminated for a cloud + shadow mask
CONTAMINATED = (0, 1, 3, 8, 9, 10, 11)

def scl_cloud_mask(
    scl_path: str,
    reflectance_paths: list[str],
    dilate_px: int = 2,
) -> np.ndarray:
    """
    Return a reflectance stack (bands, rows, cols) as float32 with
    cloud/shadow pixels set to NaN, using the SCL band as the source of truth.
    """
    # 1. Read SCL and build the contamination mask (20 m grid)
    with rasterio.open(scl_path) as scl_src:
        scl = scl_src.read(1)
    mask_20 = np.isin(scl, CONTAMINATED)

    # 2. Dilate to remove the semi-transparent cloud fringe SCL under-flags
    if dilate_px > 0:
        mask_20 = binary_dilation(mask_20, iterations=dilate_px)

    # 3. Read reflectance and resample the mask to its (10 m) shape with nearest
    bands = []
    ref_shape = None
    for path in reflectance_paths:
        with rasterio.open(path) as src:
            ref_shape = (src.height, src.width)
            bands.append(src.read(1).astype("float32"))
    stack = np.stack(bands)

    # nearest-neighbour upsample preserves categorical mask labels
    with rasterio.open(scl_path) as scl_src:
        mask_10 = scl_src.read(
            1,
            out_shape=ref_shape,
            resampling=Resampling.nearest,
        )
    mask_10 = np.isin(mask_10, CONTAMINATED)
    if dilate_px > 0:
        mask_10 = binary_dilation(mask_10, iterations=dilate_px)

    # 4. Apply — contaminated pixels become NaN across every band
    stack[:, mask_10] = np.nan
    return stack


if __name__ == "__main__":
    masked = scl_cloud_mask(
        scl_path="SCL_20m.tif",
        reflectance_paths=["B04_10m.tif", "B08_10m.tif"],
        dilate_px=2,
    )
    valid_fraction = np.isfinite(masked[0]).mean()
    print(f"Valid pixels after masking: {valid_fraction:.1%}")

Note that the mask is resampled to 10 m with Resampling.nearest. SCL holds categorical codes, so a bilinear or cubic resample would produce fractional class values — a meaningless 8.4 between cloud and vegetation — whereas nearest keeps every label intact. The dilation is applied on both the 20 m and the upsampled grid so the fringe is grown consistently regardless of which resolution you inspect. The returned valid_fraction is a cheap quality signal you can log per scene: a value near 1.0 means an almost cloud-free acquisition, while anything below roughly 0.4 usually means the scene is too contaminated to contribute much to a composite.

One subtlety worth internalising is that masking to NaN, rather than dropping pixels or filling them, is what lets the mask compose with everything downstream. A NaN pixel is ignored by an NaN-aware index calculation, skipped by a NaN-aware temporal reducer, and excluded from statistics — all without any per-step bookkeeping. The mask does its job once, at the reflectance stage, and the choice propagates for free.

Variant Patterns

1. Masking an xarray NDVI DataArray

If you compute indices with xarray, apply the mask with .where. NaN is the natural sentinel and propagates cleanly into later temporal aggregation.

import xarray as xr

scl = xr.open_dataarray("SCL_10m.tif")
contaminated = scl.isin([0, 1, 3, 8, 9, 10, 11])
ndvi_clean = ndvi.where(~contaminated)   # contaminated pixels -> NaN

The masked index is the direct input to Calculating NDVI directly from xarray DataArrays.

2. Recording a per-scene valid fraction

Before adding a scene to a composite, reject frames that are mostly cloud. The valid fraction from the mask is a cheap gate.

import numpy as np

valid_fraction = 1.0 - np.isin(scl, CONTAMINATED).mean()
if valid_fraction < 0.4:
    print("Scene >60% contaminated — skip for compositing")

3. Keeping snow when snow is the subject

For cryosphere work, drop 11 from the contaminated set so snow/ice pixels survive while clouds and shadows are still removed.

SNOW_SAFE = (0, 1, 3, 8, 9, 10)   # 11 (snow) retained
mask = np.isin(scl, SNOW_SAFE)

Common Errors

The mask has the wrong shape for the reflectance bands

SCL is delivered at 20 m and the visible/NIR bands at 10 m. Resample the SCL mask to the reflectance shape with Resampling.nearest before applying it, or you will hit a broadcast error or silently misalign the mask.

Thin cloud edges still contaminate the index

The raw SCL classes miss semi-transparent cloud halos. Dilate the boolean mask with scipy.ndimage.binary_dilation by two or three pixels so the fringe around each flagged cloud is removed as well.

Fractional class codes appear after resampling

The SCL band was resampled with bilinear or cubic interpolation, which averages categorical labels into meaningless fractions. Always use nearest-neighbour for class rasters like SCL and QA bands.

Related

• Cloud and Shadow Masking Strategies — the parent section covering QA bands, ML detectors, and mask post-processing.
• Implementing FMask and s2cloudless in Python — run your own cloud detectors when SCL is unavailable or insufficient.
• Spectral Index Calculation Pipelines — where a clean mask feeds NDVI, EVI, and NDWI computation.
• Calculating NDVI directly from xarray DataArrays — apply the mask with .where before index arithmetic.