How to Mask Clouds in Sentinel-2 with the Scene Classification Layer

This guide has 2 routes. Pick one, or read both. Each route is complete on its own.

Clouds ruin an index. A cloud is bright in every band, so NDVI over a cloud comes out low and looks like bare ground, and a map or a mean that includes cloud cells is wrong without showing any error. Every Level-2A Sentinel-2 scene comes with a layer that marks the cloud for you, the scene classification layer, SCL. This guide shows how to use it.

It assumes you have the red band, the near infrared band and the SCL of a scene, as in How to download Sentinel-2. Route A applies the mask in QGIS with the Raster Calculator. Route B does it in Python, where it is easier to compare masks and see what each choice costs.

The scene classification layer

SCL is a raster with a cell of 20 m. Each cell holds a code from 0 to 11 that names what the processing software thinks the cell is.

Code	Class	Hide it?
0	No data	Yes
1	Saturated or defective	Yes
2	Dark area pixels, including terrain shadow	Usually not
3	Cloud shadow	Yes
4	Vegetation	No
5	Not vegetated	No
6	Water	No
7	Unclassified, which includes low probability cloud	Optional
8	Cloud, medium probability	Yes
9	Cloud, high probability	Yes
10	Thin cirrus	Yes for most work
11	Snow or ice	Yes for vegetation work

The classes are the work of an automatic algorithm and it makes mistakes, as the last section explains. It is a good first filter, and it is not a perfect one.

Route A: QGIS

The mask is a division. The Raster Calculator turns a division by zero into NoData, so an expression that divides the NDVI by a term that is 1 for cells to keep and 0 for cells to hide returns NDVI where the cell is clear and NoData where it is not.

Add the three layers with Layer > Add Layer > Add Raster Layer: the red band B04 and the near infrared band B08 from IMG_DATA/R10m, and the scene classification layer SCL from IMG_DATA/R20m. Rename them B04, B08 and SCL.
The SCL has 20 m cells and the bands have 10 m cells. Make the SCL 10 m first. Open Raster > Projections > Warp (Reproject). Set Input layer to SCL, Target CRS to the CRS of B04, and Resampling method to use to Nearest Neighbour. Set Output file resolution in target georeferenced units to 10. In the advanced parameters, set Georeferenced extents of output file to be created to the extent of B04: choose Calculate from Layer and pick B04. Save the result as SCL_10m.tif.
Open Raster > Raster Calculator. In the expression box, enter the expression below. The layers appear in the list as B04@1, B08@1 and SCL_10m@1, where @1 is band 1.
Set Output layer to ndvi_masked.tif. Set the output extent from the layer B04, so that the result has the same grid. Click OK.

(("B08@1" * 0.0001 - 0.1) - ("B04@1" * 0.0001 - 0.1)) /
(("B08@1" * 0.0001 - 0.1) + ("B04@1" * 0.0001 - 0.1)) /
( ("SCL_10m@1" != 0) AND ("SCL_10m@1" != 1) AND ("SCL_10m@1" != 3) AND ("SCL_10m@1" != 8)
  AND ("SCL_10m@1" != 9) AND ("SCL_10m@1" != 10) AND ("SCL_10m@1" != 11) AND ("B04@1" != 0) )

The first two lines are NDVI, with the reflectance conversion inside it: each band is multiplied by 0.0001 and has 0.1 taken off. That is the scale and offset for scenes processed with baseline 04.00 and later, which is every scene since early 2022. For an older scene, leave out the - 0.1. Enter the expression as one line in the dialog. The last part is the mask: each comparison is 1 when the cell is not in that class, and the AND between them is 1 only when the cell is in none of them. Change the list of classes to change the mask.

Run on the scene used in the Python route below, this gives the same 681,120 clear cells and the same mean NDVI, 0.388, as Python, and the values agree cell by cell.

Route B: Python

You need Python with requests, numpy, pandas, rasterio and scipy. The scene below is a cloudy one from July 2023 over Jagraon in Punjab, which is exactly what makes it useful here.

import requests
import numpy as np
import rasterio
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds
from rasterio.enums import Resampling

# Step 1: the area, as west, south, east, north in longitude and latitude
study_area = [75.40, 30.74, 75.55, 30.84]

# Step 2: find the Sentinel-2 scenes of one day over the area, and take the one from tile T43REQ
search = {"collections": ["sentinel-2-c1-l2a"], "bbox": study_area, "datetime": "2023-07-21T00:00:00Z/2023-07-21T23:59:59Z", "limit": 5}
url = "https://earth-search.aws.element84.com/v1/search"
found = requests.post(url, json=search, timeout=60).json()["features"]
scene = [s for s in found if "T43REQ" in s["id"]][0]
print(scene["id"], "| cloud over the whole tile:", round(scene["properties"]["eo:cloud_cover"], 1), "percent")

S2A_T43REQ_20230721T054756_L2A | cloud over the whole tile: 58.1 percent

The tile is 58 percent cloudy. Read the red band, the near infrared band and SCL for the study area. The bands have a 10 m cell and SCL has a 20 m cell, so SCL is read with the shape of the red band. Resampling.nearest copies the value of the closest 20 m cell into each 10 m cell. This is the only correct choice for a layer of codes, because averaging codes would invent classes that do not exist.

# Step 1: a function that reads one window of one band. out_shape makes SCL come out on the grid of the red band.
def read(asset, out_shape=None, resampling=Resampling.bilinear):
    with rasterio.open(asset["href"]) as src:
        left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *study_area)
        window = from_bounds(left, bottom, right, top, src.transform).round_offsets().round_lengths()
        return src.read(1, window=window, out_shape=out_shape, resampling=resampling)

# Step 2: read the bands and the SCL
assets = scene["assets"]
red_numbers = read(assets["red"])
nir_numbers = read(assets["nir"])
scl = read(assets["scl"], out_shape=red_numbers.shape, resampling=Resampling.nearest)
print("window:", red_numbers.shape[1], "x", red_numbers.shape[0], "cells | SCL has the same shape:", scl.shape == red_numbers.shape)

# Step 3: whole numbers to reflectance, using the scale and offset from the metadata. 0 means no data.
def to_reflectance(numbers, asset):
    info = asset["raster:bands"][0]
    return np.where(numbers == 0, np.nan, numbers * info["scale"] + info["offset"])

red = to_reflectance(red_numbers, assets["red"])
nir = to_reflectance(nir_numbers, assets["nir"])
ndvi = (nir - red) / (nir + red)

window: 1440 x 1114 cells | SCL has the same shape: True

What is in the scene

# Step 1: the share of the window in each SCL class
names = {0: "no data", 1: "saturated", 2: "dark area", 3: "cloud shadow", 4: "vegetation", 5: "not vegetated",
         6: "water", 7: "unclassified", 8: "cloud medium", 9: "cloud high", 10: "thin cirrus", 11: "snow"}
codes, counts = np.unique(scl, return_counts=True)
for code, count in zip(codes, counts):
    print(f"{code:2d} {names[int(code)]:14s} {count / scl.size * 100:5.1f} percent")

 3 cloud shadow     0.7 percent
 4 vegetation      14.1 percent
 5 not vegetated   23.1 percent
 6 water            0.0 percent
 7 unclassified     5.2 percent
 8 cloud medium    22.2 percent
 9 cloud high       2.0 percent
10 thin cirrus     32.6 percent

The window is mostly cloud and cirrus, with vegetation and bare ground in between. Hiding the right classes is the whole job.

Build the mask

# Step 1: the classes to hide. True means hide.
hide = [0, 1, 3, 8, 9, 10, 11]
bad = np.isin(scl, hide)
print("hidden:", round(float(bad.mean()) * 100, 1), "percent of the window")

# Step 2: apply it. Hidden cells become NaN, "not a number", which numpy functions can skip.
ndvi_masked = np.where(bad, np.nan, ndvi)
print("cells left:", int(np.isfinite(ndvi_masked).sum()))

hidden: 57.5 percent of the window
cells left: 681120

How much does the mask change the answer?

import pandas as pd

# Step 1: four masks. Each one is a list of classes to hide.
masks = {
    "no mask": [],
    "clouds, shadow, no data": [0, 1, 3, 8, 9],
    "plus thin cirrus and snow": [0, 1, 3, 8, 9, 10, 11],
    "plus low probability cloud": [0, 1, 3, 7, 8, 9, 10, 11],
}

# Step 2: for each one, the share hidden and the mean NDVI of the cells that are left
rows = []
for label, classes in masks.items():
    keep = ~np.isin(scl, classes) & np.isfinite(ndvi)
    rows.append({"mask": label, "hidden %": round(100 * (1 - keep.mean()), 1), "mean NDVI": round(float(ndvi[keep].mean()), 3)})
print(pd.DataFrame(rows).to_string(index=False))

# Step 3: for comparison, the mean NDVI of the cells that the SCL calls cloud
under_cloud = np.isin(scl, [8, 9]) & np.isfinite(ndvi)
print("mean NDVI of the cloud cells themselves:", round(float(ndvi[under_cloud].mean()), 3))

                      mask  hidden %  mean NDVI
                   no mask       0.0      0.323
   clouds, shadow, no data      24.9      0.371
 plus thin cirrus and snow      57.5      0.388
plus low probability cloud      62.7      0.412
mean NDVI of the cloud cells themselves: 0.173

The table shows what the mask is worth. With no mask the mean NDVI is pulled down by the cloud cells, which are near 0.17, and every stricter mask raises the mean as more cloud and cirrus go. The strictest masks hide more than half of the window, so what is left is a smaller, cleaner sample. That is the trade. Hide too little and the cloud biases the answer. Hide too much and not enough is left to use. Thin cirrus shows it well: hiding it costs 33 percentage points of the window and lifts the mean by only 0.017. Cirrus is thin, so some of the ground still shows through. If cells are scarce, you may keep it and accept a small bias.

A reasonable default for vegetation and land cover work is the third row: no data, saturated, cloud shadow, cloud, thin cirrus and snow. This is the mask used for the monthly NDVI series in the time series tutorials. Add class 7 only when you see leftover haze in the result.

Grow the mask around the clouds

The SCL tends to miss the faint edge of a cloud and the thin haze around it. A common fix is to grow the mask by one or two cells, so that the cells next to a cloud are hidden as well.

from scipy.ndimage import binary_dilation

# Step 1: grow the hidden area by 2 cells in every direction (2 cells of 10 m is 20 m)
grown = binary_dilation(bad, iterations=2)
print("hidden before:", round(float(bad.mean()) * 100, 1), "percent | after growing by 2 cells:", round(float(grown.mean()) * 100, 1), "percent")

# Step 2: the mean NDVI with the grown mask
ndvi_grown = np.where(grown, np.nan, ndvi)
print("mean NDVI, plain mask:", round(float(np.nanmean(ndvi_masked)), 3), "| grown mask:", round(float(np.nanmean(ndvi_grown)), 3))

hidden before: 57.5 percent | after growing by 2 cells: 58.7 percent
mean NDVI, plain mask: 0.388 | grown mask: 0.39

Growing hides a little more of the window, 58.7 percent instead of 57.5, and changes the mean only in the third decimal, because this scene has big cloud areas whose edges are already caught. In a scene with many small clouds, growing matters more.

Save the masked NDVI

# Step 1: write the masked NDVI as a GeoTIFF. The transform places it, and NaN marks the hidden cells.
with rasterio.open(assets["red"]["href"]) as src:
    left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *study_area)
    window = from_bounds(left, bottom, right, top, src.transform).round_offsets().round_lengths()
    transform, crs = src.window_transform(window), src.crs

profile = {"driver": "GTiff", "height": ndvi_masked.shape[0], "width": ndvi_masked.shape[1], "count": 1,
           "dtype": "float32", "crs": crs, "transform": transform, "nodata": np.nan, "compress": "deflate"}
with rasterio.open("ndvi_masked_2023-07-21.tif", "w", **profile) as dst:
    dst.write(ndvi_masked.astype("float32"), 1)

# Step 2: read it back to check
with rasterio.open("ndvi_masked_2023-07-21.tif") as check:
    values = check.read(1)
    print("saved:", check.width, "x", check.height, "cells | CRS EPSG", check.crs.to_epsg(), "| valid cells:", int(np.isfinite(values).sum()))

saved: 1440 x 1114 cells | CRS EPSG 32643 | valid cells: 681120

What the mask cannot catch

Small clouds and thin haze. The SCL has a cell of 20 m and it can miss a cloud smaller than a few cells.
Cloud shadow. Class 3 catches part of it. The rest hides in the dark area class or in vegetation, and it lowers reflectance without looking like cloud.
Bright ground that looks like cloud, such as sand, roofs and salt, which can be called cloud by mistake. This hides good data.
Cloud that is just above the cloud probability threshold. Plain thin cloud can pass as vegetation.

So always look at the result. A map of masked NDVI should show no cloud shapes, no dark halos at their edges, and no large holes where there is bare ground. When it does, change the classes or grow the mask.

Common mistakes

Reading SCL with bilinear resampling. Averaged codes mean nothing. Use nearest.
Using the tile cloud percentage in place of the mask. The percentage says nothing about your window.
Hiding class 2, dark area, which is often just a shadowed slope or a dark field and good data.
Not checking how much of the window is hidden. A mask that hides 90 percent leaves too little to analyse.
Masking after taking the mean. Mask first, then summarise.
Hiding cloud in one scene and not in the others in a time series. Apply the same rule to every scene.