Estimating Per-Scene Cloud Cover in Python

The number that matters is cloud cover over your study area, not over the tile. Read the quality layer for the AOI window and compute it directly:

import numpy as np
import rasterio

CLOUDY = (3, 8, 9, 10)      # shadow, cloud medium, cloud high, thin cirrus

with rasterio.open(scl_href) as src:
    scl = src.read(1, window=aoi_window, out_shape=(1, 512, 512))[0]

valid = scl != 0
cloud_fraction = float(np.isin(scl, CLOUDY)[valid].mean())
print(f"AOI cloud cover: {cloud_fraction:.1%}")

This is the scene-selection step that precedes everything in Cloud and Shadow Masking Strategies.

Why This Arises in Remote Sensing Workflows

Every STAC item advertises a cloud cover percentage, and every pipeline that filters on it eventually discovers the same two problems.

The first is scale mismatch. eo:cloud_cover is computed over the whole tile — 110 by 110 km for Sentinel-2. A study area covering one percent of that tile has its own weather, and the two numbers are only loosely related. Filtering at 20 percent throws away scenes that are perfectly clear over the site, and keeps scenes whose entire cloud mass sits on it.

The second is that the tile-level number says nothing about where the cloud is or how much of the site was observed at all. A scene at the edge of a swath may cover only part of the AOI, and the unobserved part is not cloudy — it is absent, which is a different problem with a different remedy.

Computing the fraction yourself over the AOI fixes both, and it is cheap: the quality layer is one band, it can be read through an overview, and the result is a single number per scene that makes the rest of the pipeline’s decisions obvious.

Tile cloud cover is a poor proxy for site cloud cover Four scenes with tile cloud cover between 8 and 62 percent. Over a small study area the ordering is completely different: the 62 percent tile is clear over the site, and the 8 percent tile is entirely clouded there. Filtering on the tile value would keep the wrong scenes and discard the right ones. Same four scenes, two different numbers scene tile eo:cloud_cover over the 12 km study area 2023-06-05 8% 97% — unusable 2023-06-15 34% 2% — excellent 2023-06-25 62% 0% — clear, and would have been filtered out 2023-07-05 19% 41% — partly usable Use the tile value as a cheap pre-filter; make the keep-or-drop decision on the AOI number.

Environment & Setup

Package Version Why
rasterio ≥1.3.0 Windowed and decimated reads of the quality layer
pystac-client ≥0.7 Resolving quality-layer assets per item
numpy ≥1.23 Class counting
pandas ≥2.0 Ranking the candidate scenes 
pip install "rasterio>=1.3.0" "pystac-client>=0.7" "pandas>=2.0"

Complete Working Example

This function reads only the quality layer, only over the AOI, and only at a decimated resolution — then reports the three fractions a scene-selection decision needs.

Selection outcomes from the three computed fractions Coverage and cloud combine into four practical outcomes: process the scene, process it knowing part of the site is missing, defer it, or reject it. Splitting the decision this way keeps partially covered scenes from being confused with cloudy ones. Turning three numbers into one decision covered cloud over the AOI decision > 0.95 < 0.1 process — the ideal case > 0.95 0.1 – 0.4 process, record the retained fraction 0.5 – 0.95 any process only if the AOI part observed is the part you need < 0.5 any reject — the swath barely reached the site Coverage and cloud are different failures: one is geometry, the other is weather. 
import numpy as np
import pandas as pd
import rasterio
from rasterio.warp import transform_bounds
from rasterio.windows import Window, from_bounds

SCL_NODATA = 0
SCL_SHADOW = (3,)
SCL_CLOUD = (8, 9, 10)         # medium, high, thin cirrus
SCL_SNOW = (11,)

GDAL_OPTS = {"GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
             "GDAL_INGESTED_BYTES_AT_OPEN": "65536"}


def aoi_quality(scl_href: str, aoi_bounds_wgs84, *, sample: int = 512) -> dict:
    """Cloud, shadow and no-observation fractions over an AOI, from a decimated read."""
    with rasterio.Env(**GDAL_OPTS):
        with rasterio.open(scl_href) as src:
            left, bottom, right, top = transform_bounds(
                "EPSG:4326", src.crs, *aoi_bounds_wgs84, densify_pts=21)
            win = from_bounds(left, bottom, right, top, transform=src.transform)
            win = win.round_offsets(op="floor").round_lengths(op="ceil")
            win = win.intersection(Window(0, 0, src.width, src.height))
            if win.width < 1 or win.height < 1:
                return {"covered": 0.0, "cloud": None, "shadow": None, "usable": 0.0}

            # Decimate: cloud fraction is a statistic, full resolution buys nothing
            out_h = min(int(win.height), sample)
            out_w = min(int(win.width), sample)
            scl = src.read(1, window=win, out_shape=(1, out_h, out_w))[0]

    total = scl.size
    observed = scl != SCL_NODATA
    n_observed = int(observed.sum())
    if n_observed == 0:
        return {"covered": 0.0, "cloud": None, "shadow": None, "usable": 0.0}

    cloud = np.isin(scl, SCL_CLOUD)[observed].mean()
    shadow = np.isin(scl, SCL_SHADOW)[observed].mean()
    snow = np.isin(scl, SCL_SNOW)[observed].mean()

    return {
        "covered": n_observed / total,              # how much of the AOI the swath reached
        "cloud": float(cloud),
        "shadow": float(shadow),
        "snow": float(snow),
        "usable": float((1 - cloud - shadow - snow) * (n_observed / total)),
    }


def rank_scenes(items, aoi_bounds, *, asset: str = "SCL", min_usable: float = 0.7) -> pd.DataFrame:
    """Score every candidate scene and return them best-first."""
    rows = []
    for item in items:
        q = aoi_quality(item.assets[asset].href, aoi_bounds)
        rows.append({"id": item.id,
                     "date": item.datetime.date().isoformat(),
                     "tile_cloud": item.properties.get("eo:cloud_cover"),
                     **q})
    df = pd.DataFrame(rows).sort_values("usable", ascending=False)
    df["keep"] = df["usable"] >= min_usable
    return df


if __name__ == "__main__":
    ranked = rank_scenes(list(search.items()), (36.70, -1.45, 36.82, -1.33))
    print(ranked[["date", "tile_cloud", "cloud", "covered", "usable", "keep"]].head(10))
    print("scenes worth processing:", int(ranked["keep"].sum()), "of", len(ranked))

The usable fraction is the field to sort on, because it combines both failure modes: a scene can be unusable because it is cloudy or because the swath did not reach the site, and a single cloud percentage cannot express the difference.

Variant Patterns

1. Landsat QA bits instead of SCL

Landsat encodes quality as bit flags rather than class codes, so the classification step differs while everything around it stays the same.

import numpy as np

# QA_PIXEL: bit 3 = cloud, bit 4 = cloud shadow, bit 5 = snow, bit 6 = clear
qa = src.read(1, window=win, out_shape=(1, 512, 512))[0].astype("uint16")

cloud = (qa & (1 << 3)) != 0
shadow = (qa & (1 << 4)) != 0
observed = qa != 1                  # 1 = fill in Collection 2

Normalising both encodings into the same boolean vocabulary early is what lets one selection function serve both missions, as argued in Matching Landsat and Sentinel-2 Grids.

2. Reading through overviews to make the sweep cheap

Decimated reads make scene selection almost free Computing the AOI cloud fraction at full 20 metre resolution reads about 1.4 megabytes per scene. Reading through an overview into a 512 pixel array reads about 90 kilobytes and gives a cloud fraction within a tenth of a percent of the full-resolution answer. Cost of screening 200 candidate scenes full resolution 1.4 MB × 200 = 280 MB · ~3 min out_shape=(512, 512) 90 KB × 200 = 18 MB · ~25 s Cloud fraction from the decimated read: 23.4%. From the full read: 23.5%. A statistic over hundreds of thousands of pixels does not need every pixel — it needs an unbiased sample. Use full resolution only when the mask itself will be applied, not when it is only being summarised.

For a long time series, run the screen as its own pass and store the result, so that reprocessing does not repeat it.

ranked = rank_scenes(items, aoi_bounds)
ranked.to_parquet("scene_quality_2023.parquet", index=False)

# Later runs read the table instead of the archive
keep = pd.read_parquet("scene_quality_2023.parquet").query("usable >= 0.7")

Because the screen depends only on the AOI and the quality layer, it stays valid across pipeline changes, and it makes “why was this date excluded?” a query rather than an investigation.

Interpreting the Numbers

Three patterns in the output recur often enough to name.

A low cloud fraction with a low covered value means the swath clipped the AOI. The observed part may be perfectly clear, but any statistic over the site is computed from a fraction of it, and a composite built from such scenes will have a spatially uneven observation count. Treat covered as a separate acceptance criterion rather than folding it into a single score.

A moderate cloud fraction with a high shadow fraction usually indicates terrain, not weather — steep slopes at low sun elevation cast shadow that the classifier flags. In mountainous areas this can make winter scenes systematically unusable, which is a real property of the data rather than a bug in the mask.

A cloud fraction that jumps between adjacent dates while the tile value stays flat is normal and is exactly why the AOI number exists. If instead it stays flat while the imagery obviously changes, suspect that the quality asset is being read from the wrong item or that the AOI window is landing outside the swath — the failure mode covered in Cropping a STAC Item to an AOI Without Full Download.

Common Errors

The cloud fraction is always zero

The nodata class was counted as clear. Exclude class 0 (or the mission’s fill value) from the denominator before computing any fraction.

The screen is slower than processing the scenes

Full-resolution reads over the whole tile. Pass a window and an out_shape; the statistic is unchanged and the transfer drops by more than an order of magnitude.

Fractions disagree with the visible imagery

The quality layer is at a different resolution from the imagery and was read for a different window. Derive the window from the quality asset’s own transform, not the reflectance band’s.

Frequently Asked Questions

Q: Why not just use eo:cloud_cover from the STAC item? It describes the whole tile. A tile that is 60 percent cloudy can be perfectly clear over a study area occupying one percent of it, and a tile reported at 10 percent can have all of its cloud sitting exactly on your site.

Q: Can I compute this without reading full-resolution data? Yes. Read the quality layer through an overview with out_shape — cloud fraction is a statistic, and a decimated read gives the same answer to within a fraction of a percent for a small fraction of the bytes.

Q: Should shadow count as cloud? For scene selection, yes: a shadowed pixel is as unusable as a cloudy one. Report them separately as well, because a high shadow fraction with low cloud usually indicates terrain rather than weather.