⚙️ How-To Satellite imagery Raster & rasterio Satellite Imagery · step 11 of 16 Updated

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 in Snowdonia across 120 Sentinel-2 scenes, the median scene had 0.4% of its pixels usable. Building a median NDVI composite from that stack without masking gives a median NDVI of 0.042. With masking, the same 120 scenes give 0.658.

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")
44.0% usable after a 2 px buffer
Cloud masking in four steps: read the classification band at the reference resolution, select classes to drop, dilate the mask, apply it before any index.
Four steps, and the dilation is the one people skip.

Step-by-step solution

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 one scene:

no dilation   48.9% usable
1 px          46.3%
2 px          44.0%
3 px          41.9%

Two pixels at 20 m costs about 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 = np.where(usable, red, np.nan)
nir = np.where(usable, nir, 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.622, inside the healthy-vegetation range. 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 any land surface")

Bright pixels that survived the mask are cloud edge or unflagged cirrus. Whatever class they carry is the class to add to DROP next time.

Usable fraction falling from 48.9 percent to 41.9 percent as the cloud mask is dilated from zero to three pixels.
Dilation costs coverage you have plenty of and removes error you cannot detect later.

Code examples

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
   4 vegetation             42.26%  keep
   3 cloud shadow           21.49%  drop
   9 cloud high             15.75%  drop
   8 cloud medium           12.37%  drop
   5 bare soil               3.09%  keep
   6 water                   1.62%  keep
  10 thin cirrus             1.52%  drop
   2 dark area               1.30%  keep
   7 unclassified            0.60%  keep
  usable 48.9% -> 44.0% after 2 px dilation

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

This finds cloud and misses shadow entirely, which is the larger 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

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.

Explanation

Why the median composite collapses without masking

The measurement at the top is worth sitting with: an unmasked median NDVI composite over 120 scenes gave 0.042, and the masked one gave 0.658. Mean absolute difference per pixel: 0.575.

The reason is that clouds are not rare outliers here. With a median usable fraction of 0.4%, 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 Snowdonia is a large majority.

Masking changes what the median is taken over — from "all 120 observations" to "the 4 to 30 clear ones" — 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 21.5% against 28.1% for both cloud classes combined — and at a 40° sun elevation a shadow is longer than the cloud that casts it.

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 4 to 30 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.

Five checks on the output of a cloud mask, including brightness of surviving pixels and whether the class list was recorded.
The brightest surviving pixels name the class you should have dropped.

Edge cases or notes

  • 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 reaches NDVI 0.622 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.

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 120 scenes in Snowdonia: an unmasked median NDVI composite gave 0.042 and a masked one gave 0.658, a mean absolute difference of 0.575 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.