NDSI Snow Is Confused with Cloud or Water
On this page
Problem statement
A snow map built from the Normalised Difference Snow Index usually looks right until you check it against a second date or a second product. Lakes in July appear as snow. A cloud-free scene suddenly has a quarter of its snowpack flagged as cloud. A patch of "snow" disappears ten days later in weather that should not have melted it. Snow, water droplet clouds, ice clouds and open water each sit in particular places in reflectance space, and a snow test or a cloud mask that ignores one of those places mislabels it.
This page plots where Sentinel-2's own scene classes fall in NDSI and short-wave infrared on 20 May 2024 over Rocky Mountain National Park, measures how often each passes a standard snow test, and checks the confused pixels against a clear scene ten days later.
Quick answer
Compute NDSI and the snow test, then cross-tabulate against the scene classification layer (SCL) to see which classes the test and the mask disagree on:
g, nir, swir = refl["B03"], refl["B08"], refl["B11"]
ndsi = (g - swir) / (g + swir + 1e-9)
snow_test = (ndsi > 0.4) & (nir > 0.11)
for k, name in ((11, "snow"), (8, "cloud, medium probability"), (9, "cloud, high probability"), (6, "water"), (3, "cloud shadow")):
m = scl == k
print(f"{name:26} {m.sum():9,} px median NDSI {np.median(ndsi[m]):+.2f} SWIR {np.median(swir[m]):.2f} "
f"pass the snow test {snow_test[m].mean():4.0%}")snow 1,559,924 px median NDSI +0.75 SWIR 0.06 pass the snow test 79%
cloud, medium probability 503,456 px median NDSI +0.53 SWIR 0.11 pass the snow test 63%
cloud, high probability 380,024 px median NDSI +0.15 SWIR 0.47 pass the snow test 1%
water 10,336 px median NDSI +0.89 SWIR 0.01 pass the snow test 72%
cloud shadow 134,976 px median NDSI -0.38 SWIR 0.05 pass the snow test 7%
Where each surface sits, then fix water and cloud
1. Know where each surface sits
Snow is bright in green and dark in short-wave infrared: high NDSI, low SWIR. Clouds of water droplets are bright in both: low NDSI, high SWIR. Ice clouds and thin cloud over snow sit in between. Open water is dark in both bands, so its NDSI is unstable and can be high. Plotting NDSI against SWIR separates these where NDSI alone cannot.
2. Water that passes as snow
On 20 May the park's high lakes were still frozen and snow-covered, so "water" pixels with high NDSI were right. By July the lakes were open:
for d in ("20240520", "20240704"):
r, c = load(d)
nd = (r["B03"] - r["B11"]) / (r["B03"] + r["B11"] + 1e-9)
w = c == 6
print(f"{d}: water pixels {w.sum():,}; NDSI > 0.4 {((nd > 0.4) & w).mean() / w.mean():.0%}; "
f"also NIR > 0.11 {((nd > 0.4) & (r['B08'] > 0.11) & w).sum() / w.sum():.0%}")20240520: water pixels 10,336; NDSI > 0.4 95%; also NIR > 0.11 72%
20240704: water pixels 3,464; NDSI > 0.4 8%; also NIR > 0.11 1%
In July a minority of open-water pixels still have NDSI above 0.4; the near-infrared test removes nearly all of them, because snow is bright in the near infrared and water is not.
3. Snow that the mask calls cloud
The more consequential confusion runs the other way. On 20 May, Sen2Cor flagged half a million pixels as medium-probability cloud — and most of them sit squarely in the snow cluster of the NDSI–SWIR plot. High-probability cloud, by contrast, has high SWIR and fails the snow test almost entirely:
med, high = scl == 8, scl == 9
print(f"medium-probability cloud: {med.mean():.0%} of the window, {snow_test[med].mean():.0%} pass the snow test")
print(f"high-probability cloud: {high.mean():.0%} of the window, {snow_test[high].mean():.0%} pass the snow test")medium-probability cloud: 15% of the window, 63% pass the snow test
high-probability cloud: 11% of the window, 1% pass the snow test
4. Check against a clear date
Ten days later, on 30 May, the same window was almost cloud-free. If the pixels Sen2Cor called cloud on 20 May were really snow, most of them should still be snow on 30 May, at the same rate as pixels at the same elevations that were labelled snow:
r30, c30 = load("20240530")
snow30 = (((r30["B03"] - r30["B11"]) / (r30["B03"] + r30["B11"] + 1e-9)) > 0.4) & (r30["B08"] > 0.11)
clear30 = ~np.isin(c30, [0, 3, 8, 9, 10])
ok20 = ~np.isin(scl, [0, 3, 8, 9, 10])
for lo in (3000, 3200, 3400):
band = (DEM >= lo) & (DEM < lo + 200) & clear30
a = band & med & snow_test
b = band & ok20 & snow_test
print(f"{lo}-{lo + 200} m: snow on 30 May - of 'cloud' that passed the snow test {snow30[a].mean():.0%}, "
f"of clear pixels that passed {snow30[b].mean():.0%}")3000-3200 m: snow on 30 May - of 'cloud' that passed the snow test 56%, of clear pixels that passed 61%
3200-3400 m: snow on 30 May - of 'cloud' that passed the snow test 86%, of clear pixels that passed 86%
3400-3600 m: snow on 30 May - of 'cloud' that passed the snow test 88%, of clear pixels that passed 89%
The two rates match within a few points at every elevation. The pixels the mask called medium-probability cloud and the snow test called snow behaved exactly like snow: they were snow. Masking them all removed about a sixth of the snow that the 20 May scene could have shown (Example 1).
5. Decide per class
High-probability cloud and cirrus should be masked — they fail the snow test and are genuinely cloud. Medium-probability cloud over a snowpack deserves a second look: pixels that pass the NDSI and near-infrared tests with low SWIR are far more likely to be snow than cloud. Water needs the near-infrared test whatever the season.
A snow-aware cloud mask in May and April
The data used on this page
Copernicus Sentinel-2 Level-2A, tile 13TDE, 20 May, 30 May and 4 July 2024 (Copernicus Sentinel data 2024, processed by ESA), via Microsoft Planetary Computer; Copernicus DEM GLO-30 (© DLR e.V. 2010–2014 and © Airbus Defence and Space GmbH 2014–2018, provided under COPERNICUS by the European Union and ESA).
import numpy as np
import planetary_computer
import pystac_client
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import reproject, transform_bounds
from rasterio.windows import from_bounds
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace)
BOX = (-105.80, 40.22, -105.58, 40.38)
IDS = {"20240520": "S2B_MSIL2A_20240520T174909_R141_T13TDE_20240521T001144",
"20240530": "S2B_MSIL2A_20240530T174909_R141_T13TDE_20240530T231838",
"20240704": "S2A_MSIL2A_20240704T174911_R141_T13TDE_20240706T011649"}
def load(date):
"""Reflectance for B03, B08, B11 on the 10 m grid and the SCL, for BOX."""
item = catalog.get_collection("sentinel-2-l2a").get_item(IDS[date])
off = 1000 if float(item.properties["s2:processing_baseline"]) >= 4 else 0
out = {}
for band in ("B03", "B08", "B11", "SCL"):
with rasterio.open(item.assets[band].href) as src:
w = from_bounds(*transform_bounds(4326, src.crs, *BOX), src.transform).round_offsets().round_lengths()
sc = 2 if src.res[0] == 20 else 1
out[band] = src.read(1, window=w, out_shape=(int(w.height * sc), int(w.width * sc)),
resampling=Resampling.nearest).astype(float)
if src.res[0] == 10:
load.transform, load.crs = src.window_transform(w), src.crs
r, c = min(v.shape[0] for v in out.values()), min(v.shape[1] for v in out.values())
refl = {b: (out[b][:r, :c] - off) / 1e4 for b in ("B03", "B08", "B11")}
return refl, out["SCL"][:r, :c].astype(int)
refl, scl = load("20240520")
DEM = np.zeros(scl.shape, np.float32)
with rasterio.open("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N40_00_W106_00_DEM/"
"Copernicus_DSM_COG_10_N40_00_W106_00_DEM.tif") as src:
reproject(src.read(1), DEM, src_transform=src.transform, src_crs=src.crs, dst_transform=load.transform,
dst_crs=load.crs, resampling=Resampling.bilinear)
print(scl.shape, "pixels;", {int(k): int(v) for k, v in zip(*np.unique(scl, return_counts=True))})(1790, 1884) pixels; {2: 50984, 3: 134976, 4: 88692, 5: 601772, 6: 10336, 7: 42196, 8: 503456, 9: 380024, 11: 1559924}
Example 1: a snow-aware cloud mask
Keep medium-probability cloud pixels that look like snow, mask the rest:
looks_like_snow = snow_test & (swir < 0.15)
mask_strict = np.isin(scl, [0, 3, 8, 9, 10])
mask_aware = np.isin(scl, [0, 3, 9, 10]) | ((scl == 8) & ~looks_like_snow)
for name, mask in (("strict SCL mask", mask_strict), ("snow-aware mask", mask_aware)):
valid = ~mask
print(f"{name:16} masked {mask.mean():5.1%}; snow {((snow_test & valid).sum() / valid.sum()):5.1%} of the rest; "
f"snow area {(snow_test & valid).sum() * 100 / 1e6:6.1f} km²")strict SCL mask masked 30.2%; snow 55.4% of the rest; snow area 130.4 km²
snow-aware mask masked 22.6%; snow 59.7% of the rest; snow area 155.9 km²
The snow-aware mask keeps 25 km² of snow that the strict mask throws away — a sixth of the snow it finds — and still masks the high-probability cloud.
Example 2: the same check in April
The confusion is not a one-off. On 5 April, a day with scattered cloud, compared with the clear scene of 1 March:
IDS.update({"20240405": "S2A_MSIL2A_20240405T174901_R141_T13TDE_20240406T024518",
"20240301": "S2B_MSIL2A_20240301T175239_R141_T13TDE_20240301T231548"})
r5, c5 = load("20240405")
nd5 = (r5["B03"] - r5["B11"]) / (r5["B03"] + r5["B11"] + 1e-9)
t5 = (nd5 > 0.4) & (r5["B08"] > 0.11)
print(f"5 April: medium-probability cloud {t5[c5 == 8].mean():.0%} pass the snow test, high-probability {t5[c5 == 9].mean():.0%}")5 April: medium-probability cloud 77% pass the snow test, high-probability 1%
Why the cloud mask flags snow
Sen2Cor's cloud detection works partly by brightness, and a sunlit snowpack is the brightest thing in a mountain scene. When the algorithm is unsure — medium probability — it often labels bright snow as possible cloud. Its own snow class and the NDSI are built on the same physics as the snow test, but the cloud mask gets the last word in a standard workflow, so snow labelled as cloud disappears from the map as unknown.
The NDSI–SWIR plane makes the distinction clear. Water-droplet clouds are bright in the short-wave infrared, so they sit at high SWIR and low NDSI. Snow sits at low SWIR and high NDSI. Pixels in the snow cluster that the mask calls cloud are, on the evidence of a clear scene ten days later, snow. Open water is the opposite case: it can sneak into the snow cluster by being dark in both bands, and the near-infrared test removes it.
Cirrus, medium-probability cloud and water
- Plot NDSI against SWIR to see where each class sits before choosing thresholds.
- High-probability cloud and cirrus: mask them.
- Medium-probability cloud over snow: test before masking.
- Water: always add a near-infrared floor.
- Check against a nearby clear date before trusting any mask.
- Thin cirrus over snow can pass every test; the cirrus band (B10, not in L2A) or SCL class 10 helps.
Internal links
- Snow cover and NDSI explained — the index and its physics
- How to map snow cover from Sentinel-2 — the full recipe
- How to make a MODIS snow cover time series — cloud at 500 m
- How to find the snow line elevation — where masking errors bite
- How to measure glacier retreat — water in the snow test
FAQ
Why does my NDSI snow map show lakes as snow?
Water can have a high NDSI because it is dark in both green and short-wave infrared. Require a minimum near-infrared reflectance (about 0.11) as well.
Why does the Sentinel-2 cloud mask remove snow?
Sen2Cor's medium-probability cloud class often includes bright snow. On 20 May 2024 over Rocky Mountain National Park, 63% of those pixels passed the snow test, and a clear scene ten days later showed them snow-covered at the same rate as pixels labelled snow.
How do I tell snow from cloud in Sentinel-2?
Use the short-wave infrared: snow is dark at 1.6 µm, water-droplet clouds are bright. High NDSI with low SWIR is snow; low NDSI with high SWIR is cloud.
Should I mask all cloud classes before mapping snow?
Mask high-probability cloud, cirrus and shadow; test medium-probability cloud pixels against the snow criteria before discarding them.