How to Mask Clouds in Sentinel-2 Imagery with Python

Problem statement

Masking clouds is not a step you add if the scene looks cloudy. It is the step that decides whether every number downstream means anything.

Measured over a 10 km window around Snowdon across the 127 Sentinel-2 scenes from June to December 2025, the median scene had 0.2% of its pixels usable. Building a median NDVI composite from that stack without masking gives a median NDVI of 0.044. With masking, the same 127 scenes give 0.667 (Example 4).

That is not a small correction. It is the difference between "this landscape has no vegetation" and "this landscape is mostly vegetated".

Quick answer

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

# Sentinel-2 L2A scene classification codes to discard
DROP = (0,    # no data
        1,    # saturated or defective
        3,    # cloud shadow
        8,    # cloud, medium probability
        9,    # cloud, high probability
        10)   # thin cirrus

with rasterio.open("B04.tif") as ds:
    shape = ds.shape
with rasterio.open("SCL.tif") as ds:
    scl = ds.read(1, out_shape=shape, resampling=Resampling.nearest)

bad = binary_dilation(np.isin(scl, DROP), iterations=2)
usable = ~bad
print(f"{usable.mean():.1%} usable after a 2 px buffer")

39.7% usable after a 2 px buffer

That is on the 10 m grid, where two pixels are 20 m; the same two pixels on the 20 m SCL grid reach 40 m and leave less (Example 1).

Classification on your grid, named classes, dilation

1. Read the classification band onto the grid you are working on

The SCL band is 20 m; the visible and near-infrared bands are 10 m. Resample it up with nearest neighbour:

scl = ds.read(1, out_shape=shape, resampling=Resampling.nearest)

Bilinear interpolation of class codes invents classes. Halfway between 4 (vegetation) and 6 (water) is 5 (bare soil), which was never observed there.

2. Choose the classes to drop, and name the choice

DROP = (0, 1, 3, 8, 9, 10)

Three classes are decisions rather than facts:

• 2, dark area — often genuinely dark ground, often shadow the detector missed. Test both ways on your scene.
• 7, unclassified — the classifier declined. Usually cloud edge.
• 11, snow or ice — not cloud. Dropping it removes winter mountains from a time series.

Put the tuple at module level with a comment. It is the most consequential parameter in the pipeline and it should not be buried inside a function call.

3. Dilate the mask

Cloud edges are mixed pixels: part cloud, part ground. The classifier resolves them to whichever dominates, leaving a ring of half-cloud pixels labelled as ground.

Measured on the 22 September scene, on the SCL band's own 20 m grid:

with rasterio.open("SCL.tif") as ds:
    scl20 = ds.read(1)
for px in (0, 1, 2, 3):
    bad = np.isin(scl20, DROP)
    bad = binary_dilation(bad, iterations=px) if px else bad
    print(f"{px} px   {(~bad).mean():.1%} usable")

0 px   42.7% usable
1 px   40.0% usable
2 px   37.7% usable
3 px   35.7% usable

Two pixels at 20 m cost five percentage points of coverage. Skipping it keeps pixels that are systematically too bright, and a bright pixel drags a mean upward without ever looking impossible.

4. Apply the mask before any arithmetic

red_raw = rasterio.open("B04.tif").read(1).astype("float32") / 10_000   # cast before any arithmetic
nir_raw = rasterio.open("B08.tif").read(1).astype("float32") / 10_000
red = np.where(usable, red_raw, np.nan)
nir = np.where(usable, nir_raw, np.nan)

Not after. An index is a ratio, and a ratio largely survives the darkening a shadow causes — cloud shadow in this scene had a median NDVI of 0.623, inside the healthy-vegetation range (cloud masking explained measures it). Once the index is computed and averaged, the fact that the pixel was unusable is unrecoverable.

5. Check what survived

brightness = (red + nir) / 2
suspicious = usable & (brightness > 0.35)
print(f"{suspicious.sum():,} kept pixels are brighter than 0.35")

5,363 kept pixels are brighter than 0.35

Mean red and near-infrared above 0.35 is bright for this landscape, and a kept pixel that bright is worth a look: cloud edge and unflagged cirrus are the usual reasons. In this scene most of them were labelled vegetation or dark area (cloud masking explained breaks them down) — the class they carry says where the mask leaked, not always what to drop.

The mask function and its rule

The data used on this page

One Sentinel-2 L2A scene from Earth Search, 22 September 2025, cut to a 10 km window around Snowdon (EPSG:32630, 422,880–432,880 E, 5,875,440–5,885,440 N): blue, red, near-infrared, SWIR 1.6 µm and SCL as GeoTIFFs. Plus the red, near-infrared and SCL windows of all 127 scenes over the same window from June to December 2025, cached per scene and shared with the composite guide. Reflectance is DN / 10,000; the archive's COGs already carry the offset.

import os

import numpy as np
import pystac_client
import rasterio
from pyproj import Transformer
from rasterio.windows import from_bounds

BOUNDS = (422_880, 5_875_440, 432_880, 5_885_440)
to_wgs84 = Transformer.from_crs(32630, 4326, always_xy=True)
aoi = [*to_wgs84.transform(*BOUNDS[:2]), *to_wgs84.transform(*BOUNDS[2:])]
items = sorted(pystac_client.Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=aoi, datetime="2025-06-01/2025-12-31").items(), key=lambda i: i.id)


def window(href):
    with rasterio.open(href) as src:
        return src.read(1, window=from_bounds(*BOUNDS, transform=src.transform))


os.makedirs("scenes", exist_ok=True)
for item in items:
    path = f"scenes/{item.id}.npz"
    if not os.path.exists(path):
        scl = window(item.assets["scl"].href).repeat(2, axis=0).repeat(2, axis=1)
        np.savez_compressed(path, red=window(item.assets["red"].href), nir=window(item.assets["nir"].href),
                            scl=scl, date=str(item.datetime.date()),
                            cloud_cover=item.properties["eo:cloud_cover"])

scene = next(i for i in items if i.id == "S2A_30UVD_20250922_0_L2A")
for asset, name in (("blue", "B02"), ("red", "B04"), ("nir", "B08"), ("swir16", "B11"), ("scl", "SCL")):
    if not os.path.exists(f"{name}.tif"):
        with rasterio.open(scene.assets[asset].href) as src:
            win = from_bounds(*BOUNDS, transform=src.transform)
            profile = src.profile | {"width": int(win.width), "height": int(win.height), "driver": "GTiff",
                                     "transform": src.window_transform(win), "compress": "deflate"}
            with rasterio.open(f"{name}.tif", "w", **profile) as dst:
                dst.write(src.read(1, window=win), 1)
print(len(items), "scenes cached;", {n: rasterio.open(f"{n}.tif").shape for n in ("B04", "B11", "SCL")})

127 scenes cached; {'B04': (1000, 1000), 'B11': (500, 500), 'SCL': (500, 500)}

Example 1 — the mask function, with the rule returned alongside

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

SCL_NAMES = {0: "no data", 1: "saturated", 2: "dark area", 3: "cloud shadow",
             4: "vegetation", 5: "bare soil", 6: "water", 7: "unclassified",
             8: "cloud medium", 9: "cloud high", 10: "thin cirrus",
             11: "snow or ice"}
DROP = (0, 1, 3, 8, 9, 10)


def cloud_mask(scl_path, shape=None, drop=DROP, dilate_px=2, report=False):
    """Boolean mask of usable pixels, plus the rule that produced it."""
    with rasterio.open(scl_path) as ds:
        scl = (ds.read(1, out_shape=shape, resampling=Resampling.nearest)
               if shape else ds.read(1))

    bad = np.isin(scl, list(drop))
    raw_usable = float((~bad).mean())
    if dilate_px:
        bad = binary_dilation(bad, iterations=dilate_px)

    usable = ~bad
    rule = {"drop_classes": list(drop), "dilate_px": dilate_px,
            "usable_fraction": round(float(usable.mean()), 4),
            "usable_before_dilation": round(raw_usable, 4)}

    if report:
        codes, counts = np.unique(scl, return_counts=True)
        for code, count in sorted(zip(codes.tolist(), counts.tolist()),
                                  key=lambda x: -x[1]):
            print(f"  {code:2d} {SCL_NAMES.get(code, '?'):22} "
                  f"{count / scl.size:6.2%}  "
                  f"{'drop' if code in drop else 'keep'}")
        print(f"  usable {raw_usable:.1%} -> {usable.mean():.1%} "
              f"after {dilate_px} px dilation")

    return usable, rule


usable20, rule = cloud_mask("SCL.tif", report=True)
print(rule)

   4 vegetation             36.78%  keep
   3 cloud shadow           23.24%  drop
   9 cloud high             17.21%  drop
   8 cloud medium           13.99%  drop
  10 thin cirrus             2.86%  drop
   5 bare soil               2.56%  keep
   2 dark area               1.33%  keep
   6 water                   1.20%  keep
   7 unclassified            0.83%  keep
  usable 42.7% -> 37.7% after 2 px dilation
{'drop_classes': [0, 1, 3, 8, 9, 10], 'dilate_px': 2, 'usable_fraction': 0.3767, 'usable_before_dilation': 0.427}

Example 2 — masking without an SCL band

Not every product ships a classification layer. A physical fallback, in decreasing order of preference:

import numpy as np


def spectral_cloud_mask(blue, red, nir, swir16,
                        blue_threshold=0.20, ndsi_max=0.42):
    """A crude cloud mask from reflectance alone, for products with no QA band."""
    # clouds are bright in blue and bright everywhere
    bright = blue > blue_threshold

    # snow is bright too, but collapses in the shortwave infrared
    with np.errstate(divide="ignore", invalid="ignore"):
        ndsi = (blue - swir16) / (blue + swir16)
    snow = ndsi > ndsi_max

    # vegetation is never this bright
    with np.errstate(divide="ignore", invalid="ignore"):
        ndvi = (nir - red) / (nir + red)
    vegetated = ndvi > 0.3

    cloud = bright & ~snow & ~vegetated
    print(f"  bright {bright.mean():.1%}, of which snow {snow[bright].mean():.1%}, "
          f"vegetated {vegetated[bright].mean():.1%}")
    print(f"  cloud estimate {cloud.mean():.1%}")
    return cloud


def band(name, shape=(1000, 1000)):
    with rasterio.open(f"{name}.tif") as ds:     # B11 is 20 m: read it onto the 10 m grid
        return ds.read(1, out_shape=shape, resampling=Resampling.bilinear).astype("float32") / 10_000


cloud = spectral_cloud_mask(band("B02"), red_raw, nir_raw, band("B11"))
scl10 = scl                                          # the quick answer's SCL on the 10 m grid
print(f"  SCL says cloud (8, 9, 10): {np.isin(scl10, [8, 9, 10]).mean():.1%}, "
      f"shadow (3): {(scl10 == 3).mean():.1%}")
print(f"  of the spectral cloud pixels, SCL calls {np.isin(scl10[cloud], [8, 9, 10]).mean():.1%} cloud")

  bright 22.5%, of which snow 0.0%, vegetated 3.8%
  cloud estimate 21.7%
  SCL says cloud (8, 9, 10): 34.1%, shadow (3): 23.2%
  of the spectral cloud pixels, SCL calls 98.0% cloud

On the 22 September scene it flagged 21.7% as cloud where SCL has 34.1%, and 98% of what it flagged SCL also calls cloud — so it is conservative rather than wrong. It flagged none of the 23.2% that SCL calls shadow, the largest single class. Treat it as a last resort and say so in the output metadata — a spectral mask is not comparable with an SCL mask.

Example 3 — deciding whether class 2 helps or hurts

import numpy as np


def test_class_inclusion(scl, red, nir, code, baseline_drop=DROP):
    """Does keeping this class change the answer? Measure rather than guess."""
    def ndvi_median(drop):
        usable = ~np.isin(scl, list(drop))
        with np.errstate(all="ignore"):
            r = np.where(usable, red, np.nan)
            n = np.where(usable, nir, np.nan)
            total = n + r
            v = np.where(np.abs(total) < 1e-6, np.nan, (n - r) / total)
        return float(np.nanmedian(v)), float(usable.mean())

    keep_med, keep_frac = ndvi_median(baseline_drop)
    drop_med, drop_frac = ndvi_median(tuple(baseline_drop) + (code,))

    print(f"  keeping class {code}: NDVI median {keep_med:.3f}, "
          f"{keep_frac:.1%} of pixels")
    print(f"  dropping it     : NDVI median {drop_med:.3f}, "
          f"{drop_frac:.1%} of pixels")
    print(f"  effect: {drop_med - keep_med:+.3f} NDVI for "
          f"{keep_frac - drop_frac:.1%} fewer pixels")
    return drop_med - keep_med


for code in (2, 7):
    test_class_inclusion(scl10, red_raw, nir_raw, code)

  keeping class 2: NDVI median 0.710, 42.7% of pixels
  dropping it     : NDVI median 0.713, 41.4% of pixels
  effect: +0.003 NDVI for 1.3% fewer pixels
  keeping class 7: NDVI median 0.710, 42.7% of pixels
  dropping it     : NDVI median 0.712, 41.9% of pixels
  effect: +0.003 NDVI for 0.8% fewer pixels

If dropping a class moves the median NDVI noticeably, that class was contributing biased pixels and should go. If it moves nothing, keep the coverage. On the 22 September scene, dropping class 2 moved the median by +0.003 for 1.3% of the pixels and class 7 by +0.003 for 0.8% — keep both, for this scene.

Example 4 — the composite with and without the mask

unmasked_red, unmasked_nir, masked_red, masked_nir = [], [], [], []
for item in items:
    z = np.load(f"scenes/{item.id}.npz")
    r, n = z["red"].astype("float32") / 10_000, z["nir"].astype("float32") / 10_000
    ok = np.isin(z["scl"], [4, 5, 6, 7, 11])      # the composite guide's keep list: drops 2 as well
    unmasked_red.append(r)
    unmasked_nir.append(n)
    masked_red.append(np.where(ok, r, np.nan))
    masked_nir.append(np.where(ok, n, np.nan))


def median_ndvi(reds, nirs):
    with np.errstate(all="ignore"):
        r, n = np.nanmedian(np.stack(reds), axis=0), np.nanmedian(np.stack(nirs), axis=0)
        return (n - r) / (n + r)


clear = np.array([(~np.isnan(m)).mean() for m in masked_red])
without = median_ndvi(unmasked_red, unmasked_nir)
with_mask = median_ndvi(masked_red, masked_nir)
print(f"median usable fraction per scene {np.median(clear):.1%}")
print(f"median NDVI, unmasked composite  {np.nanmedian(without):.3f}")
print(f"median NDVI, masked composite    {np.nanmedian(with_mask):.3f}")
print(f"mean absolute difference         {np.nanmean(np.abs(with_mask - without)):.3f}")
print(f"clear looks per pixel            {np.percentile((~np.isnan(np.stack(masked_red))).sum(0), [0, 50, 100])}")

median usable fraction per scene 0.2%
median NDVI, unmasked composite  0.044
median NDVI, masked composite    0.667
mean absolute difference         0.587
clear looks per pixel            [ 0. 14. 31.]

Why a median collapses without masking

Why the median composite collapses without masking

The measurement at the top is worth sitting with: an unmasked median NDVI composite over 127 scenes gave 0.044, and the masked one gave 0.667. Mean absolute difference per pixel: 0.587.

The reason is that clouds are not rare outliers here. With a median usable fraction of 0.2%, the typical observation of a pixel is a cloud. A median picks the middle observation, and the middle observation is cloud. The median's robustness is a defence against a minority of bad values, and cloud over Snowdon is a large majority.

Masking changes what the median is taken over — from "all 127 observations" to "the clear ones", 14 for the median pixel — and that is why it works.

Why shadow is the biggest class

Cloud detection has physics behind it: clouds are bright, cold and high. Shadow detection has none — a shadow is darker ground, and dark ground exists.

Most shadow masks are geometric: take the detected clouds, take the sun and view angles, project each cloud to where its shadow must fall, and darken-test there. It depends on an estimated cloud height, so it is systematically less reliable than the cloud mask it derives from.

In the scene measured here, shadow covered 23.2% against 31.2% for both cloud classes combined, more than either one alone — at the scene's 36.7° sun elevation each shadow falls well clear of its cloud, onto ground the cloud does not cover.

Why dilation is not optional

A 10 m pixel on a cloud edge contains some cloud and some ground. Its reflectance is a weighted average, so it is brighter than the ground and darker than the cloud, and the classifier assigns it to whichever side is over half.

The pixels labelled ground are therefore systematically contaminated toward cloud. That is a bias, not noise, and it does not average out over many scenes — every scene has cloud edges, and they are always bright.

Two pixels of dilation cost five percentage points of coverage. With a median of 14 clear observations per pixel in the stack, coverage is the thing you have most of.

Why to record the rule

Two composites built with different DROP lists are not comparable, and nothing in the output rasters says so. A composite whose masking rule is undocumented cannot be reproduced or defended.

Write the class list, the dilation radius and the resulting usable fraction into the output file's tags or a sidecar JSON. It costs three lines.

Nearest resampling, snow and dark areas

• Resample SCL with nearest neighbour, never bilinear.
• Class 11 is snow, not cloud. Dropping it deletes winter from an alpine time series.
• Class 2 "dark area" is ambiguous — test it on your own scene rather than following a rule of thumb.
• Mask before the index, not after. Cloud shadow reached NDVI 0.623 here and cannot be detected downstream.
• binary_dilation counts pixels, not metres. Two iterations on a 20 m band is 40 m; on a 10 m band it is 20 m.
• A spectral fallback mask misses shadow. Label products that used one.
• Landsat's QA_PIXEL is a bit field, not class codes — use bitwise tests, not np.isin.
• Masking makes the array float. np.where(usable, band, np.nan) on an integer band raises or silently casts; convert first.

Internal links

• Cloud masking explained — the concepts behind the class list
• How to build a cloud-free composite from many satellite scenes — what the mask is for
• How to calculate NDVI from Sentinel-2 in Python — where the mask goes in the order of operations
• Spectral indices explained — why a ratio cannot see a shadow
• How to load Sentinel-2 bands into Python as an analysis-ready array — reading SCL onto the reference grid
• My composite has holes, stripes or grey patches — when the mask is too aggressive or not aggressive enough
• How to extract a vegetation index time series for a polygon — masking inside a time series
• Spectral bands explained — the shortwave-infrared contrast that separates cloud from snow
• SAR imagery explained — the sensor that ignores cloud entirely

FAQ

How do I mask clouds in Sentinel-2 with Python?

Read the SCL band at your working resolution with nearest-neighbour resampling, drop classes 0, 1, 3, 8, 9 and 10, dilate the result by a couple of pixels, and apply it before computing anything.

How much difference does masking make?

Measured over 127 scenes around Snowdon: an unmasked median NDVI composite gave 0.044 and a masked one gave 0.667, a mean absolute difference of 0.587 per pixel.

Should I dilate the mask?

Yes. Cloud-edge pixels are mixtures that the classifier assigns to whichever side dominates, so the ones labelled ground are biased bright. Two pixels cost about five percentage points of coverage.

What if my product has no classification band?

Use a spectral fallback — bright in blue, not snow by NDSI, not vegetated by NDVI — and accept that it will miss cloud shadow entirely. Record that the product used a fallback mask.

Should I drop class 2, "dark area"?

Test it. Compute your statistic with and without; if dropping it moves the answer, it was contributing biased pixels.

Why is my masked array still showing cloud?

Thin cirrus (class 10) and cloud edge are the usual survivors. Check the class labels of the brightest surviving pixels and extend the drop list.

Can I mask after computing an index?

No. An index is a ratio and largely cancels the darkening a shadow causes, so the information that a pixel was unusable is already gone.