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 in Snowdonia, across 120 Sentinel-2 scenes from June to December 2025:
median usable fraction of the window per scene: 0.2%
mean usable fraction: 10.2%
scenes at least 90% clear over the window: 1 of 120
One scene in 120. 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 is only −0.654.
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")48.8% of this window is usable
That scene advertised 18.8% cloud cover.
Step-by-step solution
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 120-scene stack:
filter eo:cloud_cover < 20 keeps 4 scenes
of those, >=90% clear over the window: 1
it discards 116 scenes
of those, >=90% clear over the window: 0
Here the filter did not throw away a perfect scene — but it came close. The second-clearest scene over this window was 88.1% 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:
4 vegetation 42.26%
3 cloud shadow 21.49%
9 cloud, high probability 15.75%
8 cloud, medium probability 12.37%
5 bare soil 3.09%
6 water 1.62%
10 thin cirrus 1.52%
2 dark area 1.30%
7 unclassified 0.60%
Two things stand out. Cloud shadow is the largest unusable class here, bigger than either cloud class — shadows are longer than the clouds that cast them at a 40° sun elevation. 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:
| class | keep for a vegetation index? | keep for a water mask? |
|---|---|---|
| 4 vegetation, 5 bare soil | yes | yes |
| 6 water | yes | yes |
| 11 snow / ice | usually no | no |
| 7 unclassified | your call | your call |
| 2 dark area | risky — often missed shadow | risky |
| 3 cloud shadow | no | no |
| 8, 9 cloud | no | no |
| 10 thin cirrus | no — see below | sometimes |
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, [3, 8, 9, 10])
cloudy = binary_dilation(cloudy, iterations=2) # 2 px at 20 m = 40 m
usable = ~cloudyDilating costs you real pixels — 48.9% usable falls to 44.0% at two pixels and 41.9% 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
Cloud shadow over vegetation in this scene has a median NDVI of 0.622, against 0.730 for unshadowed vegetation. Shadow suppresses red and near-infrared roughly proportionally, and a normalised difference cancels proportional changes.
That is why masking has to happen before the index, and why "the NDVI looks reasonable" is not evidence the mask worked.
Code examples
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())} 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 after 2 px dilation: 43.95%
Dilation cost 4.9 percentage points here — from 48.9% to 44.0%, 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),
} id date scene_cloud aoi_clear
S2C_30UVD_20250712_0_L2A 2025-07-12 0.1 0.9938
S2C_30UVD_20251013_0_L2A 2025-10-13 55.9 0.8805
S2B_30UVD_20250816_0_L2A 2025-08-16 33.9 0.6285
S2A_30UVD_20250816_1_L2A 2025-08-16 30.6 0.6173
S2B_30UVD_20251227_0_L2A 2025-12-27 39.8 0.5696
Only the SCL band is read — one small window per scene, not the imagery. Across 120 scenes that is a few minutes of work that 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 suspiciousBright 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 add to drop next time — usually 2 (dark area) for shadow, or 7 (unclassified) for cloud edge.
Explanation
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.654 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. Over Snowdonia there was one scene in 120 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 the larger class — 21.5% against 28.1% 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.
Edge cases or notes
- 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_covercan 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
droplists 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_covercomes 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
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 120 scenes, its correlation with actual clarity over a 10 km window was −0.654, 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 21.5% 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.622, well inside the vegetation range.
How many scenes do I need for a clear composite?
It depends entirely on the climate. Over Snowdonia the median scene was 0.2% usable and 120 scenes over six months gave every pixel several clear looks. Over a dry region, a handful of scenes is enough.