Cloud Masking in Satellite Imagery Before Analysis

A cloud is the brightest thing in most optical scenes, so an unmasked cloud dominates any index, any classifier and any difference you compute. The masking is a handful of spectral thresholds, and knowing which threshold failed is the difference between fixing a pipeline and re-running it.

Why a mask and not a filter

The instinct is to discard cloudy scenes. Scene-level cloud percentage is the wrong granularity: a scene reported as 40% cloudy still has 60% of its pixels usable, and in a maritime or tropical region discarding everything above a scene threshold discards nearly the whole archive. The unit of decision has to be the pixel.

It also has to be conservative in a specific direction. A cloud pixel mistaken for ground corrupts whatever you compute; a ground pixel mistaken for cloud only costs you an observation, and there is usually another date. Cloud masking is therefore tuned for recall on cloud, not for balanced accuracy, which is why the masks that ship with products are visibly over-inclusive at cloud edges.

The spectral tests, worked

Cloud detection rests on two physical facts: clouds are bright across the visible spectrum, and they are approximately grey — reflecting blue, green and red about equally — where most bright ground surfaces are not.

The whiteness test formalises the second. For visible-band reflectances, compute the mean and the summed absolute deviation from it, normalised:

whiteness = sum over visible bands of |rho_i - rho_mean| / rho_mean

cloud pixel:   blue 0.31  green 0.33  red 0.30
  mean = (0.31 + 0.33 + 0.30) / 3 = 0.3133
  sum  = |0.31-0.3133| + |0.33-0.3133| + |0.30-0.3133|
       =  0.0033 + 0.0167 + 0.0133 = 0.0333
  whiteness = 0.0333 / 0.3133 = 0.106

bright bare soil: blue 0.22  green 0.28  red 0.36
  mean = 0.2867
  sum  = 0.0667 + 0.0067 + 0.0733 = 0.1467
  whiteness = 0.1467 / 0.2867 = 0.512

Both pixels are bright, so a brightness test alone flags both. The whiteness index separates them by a factor of five, and a threshold around 0.7 admits the cloud and rejects the soil with a wide margin. The Fmask algorithm of Zhu and Woodcock, published in Remote Sensing of Environment in 2012, is the reference formulation of this family of tests and is the ancestor of the masks shipped with Landsat products.

Snow defeats both tests: it is bright and it is white. The separator is shortwave infrared, where snow absorbs strongly and cloud does not. The normalised difference snow index does the work:

NDSI = (green - swir1) / (green + swir1)

cloud:  green 0.33  swir1 0.28
  NDSI = (0.33 - 0.28) / (0.33 + 0.28) = 0.05 / 0.61 = 0.082

snow:   green 0.75  swir1 0.06
  NDSI = (0.75 - 0.06) / (0.75 + 0.06) = 0.69 / 0.81 = 0.852

A threshold near 0.4 separates them cleanly. Without a shortwave infrared band there is no reliable snow-cloud separation from reflectance alone, which is a hard constraint on sensors that lack one. Where a thermal band exists, as on Landsat, it adds a further strong test: cloud tops are colder than the surface beneath them, by an amount that grows with cloud height.

The cirrus band

Thin high cirrus is the case the visible tests miss. It is semi-transparent, so a cirrus-contaminated pixel is neither bright nor white — it is ground reflectance with a haze added, which shifts every band slightly and passes every threshold above.

The solution is a band placed inside a strong water-vapour absorption feature. Copernicus documents Sentinel-2 band B10 at a central wavelength of 1373.5 nm, 60 m resolution, dedicated to cirrus detection. Almost all radiation at that wavelength is absorbed by water vapour in the lower troposphere, so the surface and any low cloud are invisible; only cloud high enough to sit above most of the atmospheric water vapour reflects anything back. Any signal at all in that band is high thin cloud.

That makes it an unusually clean single-band test — no ratio, no context, just a threshold on one band. It also explains a detail that otherwise looks like an omission: B10 is not delivered in the Level-2A surface reflectance product, because it carries no surface information by design. If you are looking for it in an L2A folder and cannot find it, that is the reason.

Finding the shadow

A cloud shadow is as damaging as the cloud and much harder to detect, because dark pixels are also water, burn scars, and terrain shadow. The reliable approach is geometric: the shadow of a cloud falls at a predictable place given the cloud’s height and the sun’s position.

horizontal displacement = cloud height x tan(solar zenith angle)
direction = solar azimuth + 180 degrees

cloud base 2,000 m, solar zenith 40 deg:
  2,000 x tan(40) = 2,000 x 0.8391 = 1,678 m
  at 10 m pixels, that is 168 pixels from the cloud

Cloud height is not measured, which is the difficulty. Fmask’s approach is to project each detected cloud object across a plausible range of heights, compute where the shadow would land for each, and select the height whose projection best overlaps the pixels that are actually dark in the near-infrared. The near-infrared is used because vegetation is bright there, so shadowing produces a large relative drop, and because water — the main confuser — is dark there whether shadowed or not and can be excluded separately.

Both cloud and shadow masks then get dilated by a few pixels. The reason is optical rather than algorithmic: the sensor’s point spread function mixes neighbouring ground samples, so the pixels immediately outside a cloud boundary contain some cloud signal without triggering any threshold. The cost is quantifiable — a cloud with a 4 km perimeter dilated by 3 pixels at 10 m loses an extra 4,000 × 30 = 120,000 m², about 0.12 km² — and it is worth paying, because edge pixels are the ones most likely to look like real but subtle change.

The masks that ship with the data

You rarely need to implement any of this. Both major open archives ship a per-pixel classification, and knowing their encodings is most of the practical work.

Sentinel-2 Level-2A carries a Scene Classification layer. Copernicus documents its values as 0 no data, 1 saturated or defective, 2 cast shadows, 3 cloud shadows, 4 vegetation, 5 not vegetated, 6 water, 7 unclassified, 8 cloud medium probability, 9 cloud high probability, 10 thin cirrus, and 11 snow or ice. A usable mask is typically classes 3, 8, 9 and 10, and whether you include 2 and 11 depends on the analysis — snow is a legitimate surface for a land-cover map and a contaminant for a vegetation index.

Landsat Collection 2 packs its flags into the bits of a single QA_PIXEL band. The USGS documents the layout as: bit 0 fill, bit 1 dilated cloud, bit 2 cirrus, bit 3 cloud, bit 4 cloud shadow, bit 5 snow, bit 6 clear, bit 7 water, then confidence values in bit pairs — 8 to 9 cloud, 10 to 11 cloud shadow, 12 to 13 snow or ice, 14 to 15 cirrus — each encoding none, low, medium or high.

import numpy as np, rasterio

with rasterio.open(qa_path) as src:
    qa = src.read(1)

fill     = (qa & (1 << 0)) != 0
dilated  = (qa & (1 << 1)) != 0
cirrus   = (qa & (1 << 2)) != 0
cloud    = (qa & (1 << 3)) != 0
shadow   = (qa & (1 << 4)) != 0

# bits 8-9: cloud confidence, 3 == high
cloud_conf = (qa >> 8) & 0b11

usable = ~(fill | dilated | cirrus | cloud | shadow)

The shift-and-mask pattern is the whole of it, and the mistake to avoid is reading the confidence pairs as single bits. A pixel with cloud confidence 2 (medium) has bit 9 set and bit 8 clear, so a careless single-bit test on bit 8 reports low confidence for it.

What masking does to a composite

Masking leaves holes, and the standard repair is to composite several dates — take the per-pixel median of every unmasked observation in a window. The coverage arithmetic is favourable:

assume each scene is clear at a given pixel with probability 0.4
12 monthly scenes:

P(no clear observation) = 0.6^12 = 0.00218

so about 0.2% of pixels have no valid observation all year

That is the good news and it hides the bad news. The surviving observations are not a random sample of the year: they are the clear days. In a monsoon climate that means the composite is built almost entirely from the dry season, so a “2026 annual composite” is a dry-season image with an annual label. Any comparison between two such composites is partly a comparison of which months happened to be clear in each year.

Two habits address it. Carry a per-pixel count of contributing observations alongside the composite, so downstream work can see where the estimate rests on one scene and where it rests on ten. And keep the acquisition dates rather than collapsing to a label, because when the composite is used for change detection between years the seasonal composition of each side is the first thing to check. The same care applies to any index computed from the composite, including NDVI, where a seasonal shift in which observations survived is indistinguishable from a change in vegetation.