Cloud Masking Explained: QA Bands, Scene Classification and What They Miss

Problem statement

Every optical satellite analysis has the same first problem: most pixels are not of the ground. They are of cloud, of cloud shadow, of thin cirrus that looks like haze, or of the edge of a swath.

The scale of it surprises people. 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 the window usable and only one was at least 90% clear (step 1). And the metadata people filter on does not identify it: the correlation between a scene's advertised eo:cloud_cover and how clear it actually was over this window was only −0.656.

Quick answer

Use the per-pixel classification band, not the scene-level cloud percentage:

import numpy as np
import rasterio

# Sentinel-2 L2A scene classification
UNUSABLE = {0,   # no data
            1,   # saturated or defective
            3,   # cloud shadow
            8,   # cloud, medium probability
            9,   # cloud, high probability
            10}  # thin cirrus

with rasterio.open("SCL.tif") as ds:
    scl = ds.read(1)

usable = ~np.isin(scl, list(UNUSABLE))
print(f"{usable.mean():.1%} of this window is usable")

42.7% of this window is usable

That scene, 22 September 2025, advertised 18.8% cloud cover.

Beyond scene cloud cover, the classification band

1. Stop filtering on scene-level cloud cover alone

eo:cloud_cover is one number for a whole 110 km tile. It answers "is this scene worth storing", not "is my field visible".

Measured over the 127-scene stack:

GOOD = [4, 5, 6, 7, 11]
clear = np.array([np.isin(scl, GOOD).mean() for scl in stack_scl])
print(f"median usable fraction of the window per scene: {np.median(clear):6.1%}")
print(f"mean usable fraction:                           {clear.mean():6.1%}")
print(f"scenes at least 90% clear over the window:      {(clear >= 0.9).sum()} of {len(clear)}")
print(f"correlation with eo:cloud_cover:                {np.corrcoef(cloud_cover, clear)[0, 1]:.3f}")
keep = cloud_cover < 20
print(f"filter eo:cloud_cover < 20  keeps {keep.sum():3d} scenes, {(clear[keep] >= 0.9).sum()} of them >=90% clear")
print(f"it discards {(~keep).sum()} scenes, {(clear[~keep] >= 0.9).sum()} of them >=90% clear")

median usable fraction of the window per scene:   0.2%
mean usable fraction:                            10.3%
scenes at least 90% clear over the window:      1 of 127
correlation with eo:cloud_cover:                -0.656
filter eo:cloud_cover < 20  keeps   4 scenes, 1 of them >=90% clear
it discards 123 scenes, 0 of them >=90% clear

Here the filter did not throw away a perfect scene — but it came close. The second-clearest scene over this window, 13 October, was 88% clear and carried a scene cloud cover of 55.9%. Any threshold below 56% discards it.

Use the scene-level number as a cheap pre-filter with a generous threshold, then decide properly per pixel.

2. Read the classification band and know its classes

Sentinel-2 L2A ships a 20 m scene classification layer. Over the AOI on 2025-09-22 (Example 1 prints this):

 4  vegetation                   36.78%
 3  cloud shadow                 23.24%
 9  cloud, high probability      17.21%
 8  cloud, medium probability    13.99%
10  thin cirrus                   2.86%
 5  bare soil                     2.56%
 2  dark area                     1.33%
 6  water                         1.20%
 7  unclassified                  0.83%

Two things stand out. Cloud shadow is the largest unusable class here, bigger than either cloud class — at the scene's 36.7° sun elevation a cloud's shadow falls well away from it, onto ground the cloud itself does not cover. And "dark area" and "unclassified" are ambiguous: real ground the classifier could not label, or missed shadow.

3. Decide what "usable" means and write it down

The classes are not simply good or bad; the boundary is a decision:

Whatever you decide, put the class list in the code as a named constant and in the output metadata. A composite whose mask rule is undocumented cannot be reproduced or compared.

4. Buffer the mask

Cloud masks are optimistic at the edges. A cloud edge pixel is a mixture of cloud and ground, and the classifier resolves it to whichever dominates — leaving a ring of half-cloud pixels labelled as ground.

from scipy.ndimage import binary_dilation

cloudy = np.isin(scl, [0, 1, 3, 8, 9, 10])
for px in (0, 1, 2, 3):                          # 2 px at 20 m = 40 m
    grown = binary_dilation(cloudy, iterations=px) if px else cloudy
    print(f"dilated {px} px: {(~grown).mean():.1%} usable")

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

Dilating costs you real pixels: 42.7% usable falls to 37.7% at two pixels and 35.7% at three. Not dilating costs you wrong pixels, which is worse, because a half-cloud pixel is bright and drags a mean upward without ever looking impossible.

5. Remember the index cannot see what the mask missed

A normalised difference cancels proportional changes, and shadow suppresses red and near-infrared by a similar factor. Compare the shadowed pixels with the vegetation pixels in this scene:

from rasterio.enums import Resampling

with rasterio.open("SCL.tif") as ds:
    scl10 = ds.read(1, out_shape=(1000, 1000), resampling=Resampling.nearest)
red = rasterio.open("B04.tif").read(1).astype("float32") / 10_000
nir = rasterio.open("B08.tif").read(1).astype("float32") / 10_000
with np.errstate(all="ignore"):
    ndvi = (nir - red) / (nir + red)
for code, name in [(4, "vegetation"), (3, "cloud shadow"), (5, "bare soil"), (9, "cloud, high probability")]:
    sel = scl10 == code
    print(f"{name:24} NDVI median {np.nanmedian(ndvi[sel]):.3f}   "
          f"NIR {np.median(nir[sel]):.3f}   red {np.median(red[sel]):.3f}")

vegetation               NDVI median 0.726   NIR 0.362   red 0.057
cloud shadow             NDVI median 0.623   NIR 0.133   red 0.031
bare soil                NDVI median 0.362   NIR 0.285   red 0.126
cloud, high probability  NDVI median 0.085   NIR 0.597   red 0.500

Shadow cut the near-infrared by more than half and the NDVI far less.

That is why masking has to happen before the index, and why "the NDVI looks reasonable" is not evidence the mask worked.

A mask function with explicit decisions

The data used on this page

The same 10 km window around Snowdon (EPSG:32630, 422,880–432,880 E, 5,875,440–5,885,440 N) in every Sentinel-2 L2A scene in Earth Search from 1 June to 31 December 2025 — 127 scenes — plus the red, near-infrared and SCL bands of one of them, 22 September 2025, as GeoTIFFs. The stack cache is shared with the composite guide; only its SCL arrays are used here.

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:])]
catalog = pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
items = sorted(catalog.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"])
stack_scl = [np.load(f"scenes/{item.id}.npz")["scl"] for item in items]
cloud_cover = np.array([item.properties["eo:cloud_cover"] for item in items])

scene = next(i for i in items if i.id == "S2A_30UVD_20250922_0_L2A")
for asset, name in (("red", "B04"), ("nir", "B08"), ("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;", scene.id, "advertises", round(scene.properties["eo:cloud_cover"], 1), "% cloud")

127 scenes; S2A_30UVD_20250922_0_L2A advertises 18.8 % cloud

Example 1 — a mask function with the decisions made explicit

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"}


def usable_mask(scl_path, shape=None, drop=(0, 1, 3, 8, 9, 10),
                dilate_px=2, report=True):
    """Boolean mask of pixels worth analysing, with the rule stated up front."""
    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))
    if dilate_px:
        bad = binary_dilation(bad, iterations=dilate_px)
    usable = ~bad

    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]):
            mark = "drop" if code in drop else "keep"
            print(f"  {code:2d} {SCL_NAMES.get(code, '?'):24} "
                  f"{count / scl.size:6.2%}  {mark}")
        print(f"  usable after {dilate_px} px dilation: {usable.mean():.2%}")

    return usable, {"drop_classes": list(drop), "dilate_px": dilate_px,
                    "usable_fraction": float(usable.mean())}


usable20, rule = usable_mask("SCL.tif")

   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 after 2 px dilation: 37.67%

Dilation cost 5.0 percentage points here — from 42.7% to 37.7%, and 7 points at three pixels. That is the price of not averaging cloud edges into your answer.

Example 2 — ranking scenes by clarity over your area

import numpy as np
import rasterio
from rasterio.warp import transform_bounds


def aoi_clarity(item, aoi_wgs84, good=(4, 5, 6, 7, 11)):
    """How clear is this scene over the area you care about?"""
    with rasterio.open(item.assets["scl"].href) as ds:
        bounds = transform_bounds("EPSG:4326", ds.crs, *aoi_wgs84)
        window = rasterio.windows.from_bounds(
            *bounds, transform=ds.transform).round_offsets().round_lengths()
        scl = ds.read(1, window=window)

    return {
        "id": item.id,
        "date": item.properties["datetime"][:10],
        "scene_cloud": round(item.properties["eo:cloud_cover"], 1),
        "aoi_clear": round(float(np.isin(scl, list(good)).mean()), 4),
    }


ranked = sorted((aoi_clarity(item, aoi) for item in items), key=lambda r: -r["aoi_clear"])
print(f"  {'id':25} {'date':11} {'scene_cloud':>11}  aoi_clear")
for r in ranked[:5]:
    print(f"  {r['id']:25} {r['date']:11} {r['scene_cloud']:11.1f}  {r['aoi_clear']:9.4f}")

  id                        date        scene_cloud  aoi_clear
  S2C_30UVD_20250712_0_L2A  2025-07-12          0.1     0.9950
  S2C_30UVD_20251013_0_L2A  2025-10-13         55.9     0.8827
  S2A_30UVD_20250816_1_L2A  2025-08-16         30.6     0.6305
  S2B_30UVD_20250816_0_L2A  2025-08-16         33.9     0.6153
  S2B_30UVD_20251227_0_L2A  2025-12-27         39.8     0.5840

Only the SCL band is read — one small window per scene, not the imagery. Across 127 scenes it took six minutes here, and it saves downloading scenes you would throw away.

Example 3 — auditing what your mask left behind

import numpy as np


def audit_mask(usable, red, nir, scl):
    """Did anything obviously cloud-like survive the mask?"""
    brightness = (red + nir) / 2
    kept = brightness[usable & np.isfinite(brightness)]
    p99 = float(np.percentile(kept, 99))

    suspicious = usable & (brightness > 0.35)     # bright for a land surface
    print(f"  kept p99 brightness      {p99:.3f}")
    print(f"  kept pixels above 0.35   {int(suspicious.sum()):,} "
          f"({suspicious.mean():.2%})")

    if suspicious.sum():
        codes, counts = np.unique(scl[suspicious], return_counts=True)
        for code, count in zip(codes.tolist(), counts.tolist()):
            print(f"     class {code}: {count:,}")
    return suspicious


usable10, _ = usable_mask("SCL.tif", shape=(1000, 1000), report=False)
suspicious = audit_mask(usable10, red, nir, scl10)

  kept p99 brightness      0.365
  kept pixels above 0.35   5,363 (0.54%)
     class 2: 1,308
     class 4: 2,418
     class 5: 1,068
     class 6: 14
     class 7: 555

Bright pixels that survived the mask are almost always cloud edge or unflagged thin cirrus. Reporting which class they were labelled as tells you which class to look at next. In the 22 September scene, 5,363 kept pixels (0.54%) were brighter than 0.35: 2,418 labelled vegetation, 1,308 dark area, 1,068 bare soil and 555 unclassified. A bright pixel labelled "dark area" is a contradiction worth a look; bright vegetation next to a cloud is the cloud-edge mixture that dilation is for.

Why scene percentages fail and shadow is hard

Why a scene-level percentage cannot work

A Sentinel-2 tile is 110 × 110 km. Cloud is not uniform over that area — it organises into fronts, convective cells and orographic caps, all of which are tens of kilometres across.

So a 30% cloudy tile might be completely clear over your 10 km window, or completely covered. The measured correlation of −0.656 between scene cloud cover and window clarity is exactly what that looks like: real information, nowhere near enough to select on.

The asymmetry matters too. Filtering hard on scene cloud cover throws away scenes that were clear over your area, and those are the scarce resource. Around Snowdon there was one scene in 127 that was more than 90% clear over the window; a filter that discards it has cost you the whole period.

Why cloud shadow is the harder half

Cloud detection has a strong physical signal: clouds are bright, cold, and high. Shadow detection has none of that. A shadow is just darker ground, and dark ground exists — water, wet rock, dense conifer, north-facing slopes.

Most shadow detectors work geometrically: take the detected clouds, take the solar and viewing geometry, project each cloud to where its shadow must fall, and look for darkening there. That works when cloud height is known, and cloud height is estimated.

The consequence is that shadow masks are systematically less reliable than cloud masks, and shadow is often a large class — 23.2% against 31.2% of cloud in the scene above, and larger than either individual cloud class.

Why thin cirrus is its own problem

Cirrus is semi-transparent. The ground is visible through it, but with reduced contrast and a bluish bias, because thin ice scatters short wavelengths more.

That leaves a genuine choice. Dropping cirrus loses pixels that carry real, if degraded, signal. Keeping it admits pixels with a systematic bias that will look like a real change when compared against a clear date. For anything comparing across dates, drop it; for a single-date visual product, keeping it may be fine.

Sentinel-2 has a dedicated cirrus band at 1.375 µm, in a water-vapour absorption region where the lower atmosphere is opaque — anything bright there is high and thin.

Why more scenes beat a better mask

Given that the median scene here is 0.2% usable, no mask refinement will produce a clear picture of one date. The only way out is stacking: with enough dates, most pixels are seen clear at least once, and a per-pixel composite fills the map.

That reframes cloud masking. Its job is not to rescue a scene; it is to label pixels honestly so a compositing step can choose among many observations. A mask that is too permissive poisons the composite; a mask that is slightly too aggressive just costs you observations you have plenty of.

SCL resolution, dark areas and snow

• The SCL band is 20 m, the visible bands are 10 m. Resample with nearest neighbour, never bilinear — interpolating class codes invents classes that do not exist.
• Class 2 "dark area" is ambiguous. Often genuinely dark ground, often missed shadow. Test both ways on your scene.
• Class 11 is snow or ice, not cloud. Dropping it silently removes winter mountains from your time series.
• Dilate the mask. Two pixels at 20 m cost about 5 percentage points of coverage here and remove the cloud-edge mixtures.
• Snow and cloud are hard to separate in the visible bands and easy in the shortwave infrared, where ice absorbs and water droplets do not.
• eo:cloud_cover can be missing or wrong. Some providers compute it over the tile including no-data area, which inflates it for edge scenes.
• A mask is per-sensor. Landsat's QA_PIXEL band is a bit field, not a class code, and needs bitwise tests rather than np.isin.
• Record the mask rule with the output. Two composites with different drop lists are not comparable.

Internal links

• How to mask clouds in Sentinel-2 imagery with Python — the implementation
• How to build a cloud-free composite from many satellite scenes — what masking is actually for
• Spectral indices explained — why an index cannot see the shadow the mask missed
• Spectral bands explained — the shortwave-infrared contrast that separates cloud from snow
• How to search and download satellite imagery with STAC in Python — where eo:cloud_cover comes from
• My composite has holes, stripes or grey patches — what an over- or under-aggressive mask produces
• How to extract a vegetation index time series for a polygon — masking inside a time series
• Sentinel-2 or Landsat? Choosing a satellite imagery source — two different mask formats
• How to detect new forest clearing with Sentinel-2 — masks, NDVI drops and validation
• Clouds show as cold spots — what unmasked cloud does to a heat map

FAQ

What is the SCL band?

Sentinel-2 L2A's scene classification layer: a 20 m raster where each pixel carries a class code from 0 to 11, including cloud, cloud shadow, cirrus, vegetation, water and snow.

Which SCL classes should I drop?

Start with 0, 1, 3, 8, 9 and 10 — no data, saturated, cloud shadow, both cloud classes and thin cirrus. Treat 2 (dark area) and 7 (unclassified) as decisions to test on your own scene.

Can I just filter on eo:cloud_cover?

Only as a coarse pre-filter. Measured over 127 scenes, its correlation with actual clarity over a 10 km window was −0.656, and the second-clearest scene over the window advertised 55.9% cloud.

Why is cloud shadow the biggest problem?

Because it has no distinctive physical signature — it is just darker ground — so it is detected geometrically from estimated cloud heights. In the scene measured here it covered 23.2% of the window, more than either cloud class.

Should I dilate the cloud mask?

Yes. Cloud edges are mixed pixels that the classifier resolves to whichever side dominates. Two pixels of dilation cost about 5 percentage points of coverage here and remove the mixtures.

Why does my NDVI look fine over an area I know was cloudy?

Because NDVI is a ratio, and shadow reduces red and near-infrared proportionally. Cloud shadow in this scene had a median NDVI of 0.623, well inside the vegetation range.

How many scenes do I need for a clear composite?

It depends entirely on the climate. Around Snowdon the median scene was 0.2% usable, and 127 scenes over seven months gave the median pixel 14 clear looks. Over a dry region, a handful of scenes is enough.