How to Build a Cloud-Free Composite from Many Satellite Scenes

Problem statement

Over a cloudy region, no single satellite date shows you the ground. Across the 127 Sentinel-2 scenes that cover a 10 km window around Snowdon between June and December 2025, the median scene had 0.2% of the window usable, the best had 99.5%, and 7 of the 127 were at least half usable (step 1). One scene in eighteen was half usable. A composite is not an optimisation for these conditions — it is the only way to get a complete picture at all.

The good news is on the other axis. Stacked, those same 127 scenes gave a median pixel 14 clear observations and the best-observed 31; 904 pixels out of a million — 0.09% — were never seen clear. There is plenty of data; it is just never all on the same day.

Quick answer

Mask, then take a per-pixel median over time:

import numpy as np

# red_scenes, nir_scenes: (y, x) reflectance arrays, already masked to NaN
red = np.stack(red_scenes)
nir = np.stack(nir_scenes)

with np.errstate(all="ignore"):
    red_c = np.nanmedian(red, axis=0)
    nir_c = np.nanmedian(nir, axis=0)
    total = nir_c + red_c
    ndvi = np.where(np.abs(total) < 1e-6, np.nan, (nir_c - red_c) / total)

n_obs = np.isfinite(red).sum(axis=0)
print(f"observations per pixel: min {n_obs.min()}, "
      f"median {np.median(n_obs):.0f}, max {n_obs.max()}")

observations per pixel: min 0, median 14, max 31

The masking is what makes this work. Measured on the same stack (step 4):

Select by local clarity, one grid, a masked median

1. Select scenes on clarity over your area, not on scene cloud cover

Read only the 20 m classification band for each candidate — a few kilobytes each — and rank on the fraction of your window that is usable:

clear = np.array([np.isin(scl, [4, 5, 6, 7, 11]).mean() for scl in scl_scenes])
order = np.argsort(-clear)
print(f"median usable fraction per scene {np.median(clear):6.1%}")
print(f"best single scene                {clear.max():6.1%}")
print(f"scenes at least 50% usable       {(clear >= 0.5).sum():3d} of {len(clear)}")
print(f"correlation with eo:cloud_cover  {np.corrcoef(cloud_cover, clear)[0, 1]:6.3f}")
for i in order[:3]:
    print(f"  {dates[i]}  {clear[i]:6.1%} usable, scene says {cloud_cover[i]:4.1f}% cloud")

median usable fraction per scene   0.2%
best single scene                 99.5%
scenes at least 50% usable         7 of 127
correlation with eo:cloud_cover  -0.656
  2025-07-12   99.5% usable, scene says  0.1% cloud
  2025-10-13   88.1% usable, scene says 55.9% cloud
  2025-08-16   62.7% usable, scene says 30.6% cloud

Scene-level eo:cloud_cover correlated only −0.66 with this, and the second-clearest scene over the window advertised 55.9% cloud.

Do not throw away partly cloudy scenes. A scene that is 30% usable over your window still contributes those pixels, and in a stack that is exactly what you need.

2. Put every scene on one grid

Every scene must share a CRS, shape and transform before stacking. Read each band with out_shape against a single reference, and assert it:

assert all(m["transform"] == reference["transform"] for m in metas)

A scene from a neighbouring UTM zone will have the same date and a different CRS. Stacking it silently produces a smear.

3. Mask every scene before stacking

usable = ~np.isin(scl, (0, 1, 3, 8, 9, 10))
red = np.where(usable, red, np.nan)

This is the step that produces the 0.6 NDVI difference in the table above. Without it, the typical observation of a pixel is cloud, and a median over all observations returns cloud.

4. Reduce over time, and prefer the median

def ndvi(red, nir):
    with np.errstate(all="ignore"):
        total = nir + red
        return np.where(np.abs(total) < 1e-6, np.nan, (nir - red) / total)


raw_red, raw_nir = np.stack(raw_red_scenes), np.stack(raw_nir_scenes)
for label, r, n in [("unmasked", raw_red, raw_nir), ("masked", red, nir)]:
    for name, reduce in [("mean", np.nanmean), ("median", np.nanmedian)]:
        with np.errstate(all="ignore"):
            v = ndvi(reduce(r, axis=0), reduce(n, axis=0))
        p5, p50, p95 = np.nanpercentile(v, [5, 50, 95])
        print(f"{name + ', ' + label:17} NDVI median {p50:.3f}   p5 {p5:.3f}   p95 {p95:.3f}")

mean, unmasked    NDVI median 0.064   p5 0.027   p95 0.090
median, unmasked  NDVI median 0.044   p5 0.016   p95 0.083
mean, masked      NDVI median 0.543   p5 0.194   p95 0.734
median, masked    NDVI median 0.667   p5 0.274   p95 0.821

The masked mean sits 0.12 below the masked median. Both are computed over the same clear observations, so the gap is residual contamination — cloud edges and unflagged thin cirrus that survived the mask. The median ignores them; the mean averages them in.

Use the mean only when you have already established that the mask is essentially perfect, which is rarely.

5. Return the observation count with the composite

n_obs = np.isfinite(red).sum(axis=0)
composite = np.where(n_obs >= 3, composite, np.nan)

A pixel with 3 observations and a pixel with 31 look identical in the output raster. They are not the same measurement, and any user of the composite needs to know which they have.

6. Understand what a composite is not

The best single scene here, 12 July, gave an NDVI median of 0.765; the June-to-December composite gave 0.667. Neither is wrong — they answer different questions. The single scene is one moment in July at peak growth; the composite is a summary across a season that includes senescence.

A composite has no date. It is a statistic over a period, and it must be labelled with that period.

A composite with its bookkeeping

The data used on this page

Every Sentinel-2 L2A scene in Earth Search that covers a 10 km window around Snowdon (EPSG:32630, 422,880–432,880 E, 5,875,440–5,885,440 N) between 1 June and 31 December 2025: 127 scenes of tile 30UVD, from Sentinel-2A, 2B and 2C. Only the red and near-infrared bands and the scene classification (SCL) are read, one window each, and cached per scene. The archive's COGs already carry the reflectance offset, so reflectance is DN / 10,000.

import os

import numpy as np
import pystac_client
import rasterio
from pyproj import Transformer
from rasterio.windows import from_bounds

BOUNDS = (422_880, 5_875_440, 432_880, 5_885_440)
to_wgs84 = Transformer.from_crs(32630, 4326, always_xy=True)
bbox = [*to_wgs84.transform(*BOUNDS[:2]), *to_wgs84.transform(*BOUNDS[2:])]
items = sorted(pystac_client.Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=bbox, datetime="2025-06-01/2025-12-31").items(),
    key=lambda i: i.id)


def window(href):
    with rasterio.open(href) as src:
        return src.read(1, window=from_bounds(*BOUNDS, transform=src.transform))


os.makedirs("scenes", exist_ok=True)
for item in items:
    path = f"scenes/{item.id}.npz"
    if not os.path.exists(path):
        scl = window(item.assets["scl"].href).repeat(2, axis=0).repeat(2, axis=1)   # 20 m onto the 10 m grid
        np.savez_compressed(path, red=window(item.assets["red"].href), nir=window(item.assets["nir"].href),
                            scl=scl, date=str(item.datetime.date()),
                            cloud_cover=item.properties["eo:cloud_cover"])

USABLE = [4, 5, 6, 7, 11]              # vegetation, bare, water, unclassified, snow
scl_scenes, dates, cloud_cover = [], [], []
raw_red_scenes, raw_nir_scenes, red_scenes, nir_scenes = [], [], [], []
for item in items:
    z = np.load(f"scenes/{item.id}.npz")
    red, nir = z["red"].astype("float32") / 10_000, z["nir"].astype("float32") / 10_000
    usable = np.isin(z["scl"], USABLE)
    scl_scenes.append(z["scl"])
    dates.append(str(z["date"]))
    cloud_cover.append(float(z["cloud_cover"]))
    raw_red_scenes.append(red)
    raw_nir_scenes.append(nir)
    red_scenes.append(np.where(usable, red, np.nan))
    nir_scenes.append(np.where(usable, nir, np.nan))
print(len(items), "scenes,", red_scenes[0].shape, "pixels each,", dates[0], "to", dates[-1])

127 scenes, (1000, 1000) pixels each, 2025-06-04 to 2025-12-29

Example 1 — a composite with the bookkeeping included

import numpy as np


def composite(scenes, method="median", min_obs=3):
    """Per-pixel reduction over time with the observation count returned.

    scenes: list of dicts of masked float arrays, all on the same grid.
    """
    bands = sorted(set().union(*(s.keys() for s in scenes)))
    stacks = {b: np.stack([s[b] for s in scenes]) for b in bands}

    reference = stacks[bands[0]]
    n_obs = np.isfinite(reference).sum(axis=0)

    reducer = {"median": np.nanmedian, "mean": np.nanmean,
               "max": np.nanmax, "min": np.nanmin}[method]

    out = {}
    with np.errstate(all="ignore"):
        for band, stack in stacks.items():
            values = reducer(stack, axis=0)
            out[band] = np.where(n_obs >= min_obs, values, np.nan)

    meta = {
        "method": method, "n_scenes": len(scenes), "min_obs": min_obs,
        "obs_min": int(n_obs.min()), "obs_median": float(np.median(n_obs)),
        "obs_max": int(n_obs.max()),
        "pixels_below_min_obs": int((n_obs < min_obs).sum()),
    }
    print(f"  {method} of {len(scenes)} scenes, "
          f"observations {meta['obs_min']}-{meta['obs_max']} "
          f"(median {meta['obs_median']:.0f}), "
          f"{meta['pixels_below_min_obs']:,} pixels below the minimum")
    return out, n_obs, meta


out, n_obs, meta = composite([{"red": r, "nir": n} for r, n in zip(red_scenes, nir_scenes)])

  median of 127 scenes, observations 0-31 (median 14), 11,460 pixels below the minimum

Example 2 — a "best pixel" composite instead of a statistic

import numpy as np


def best_pixel_composite(stacks, score, bands):
    """Pick one whole observation per pixel rather than mixing them.

    score: (time, y, x) array, higher is better — e.g. NDVI, or negative
    cloud distance. Keeps the bands spectrally consistent per pixel.
    """
    score = np.where(np.isfinite(score), score, -np.inf)
    best = np.argmax(score, axis=0)                     # (y, x)
    valid = np.isfinite(score).any(axis=0)

    ys, xs = np.indices(best.shape)
    out = {b: np.where(valid, stacks[b][best, ys, xs], np.nan) for b in bands}

    chosen, counts = np.unique(best[valid], return_counts=True)
    print(f"  {len(chosen)} of {score.shape[0]} scenes contribute pixels")
    print(f"  largest single contribution: {counts.max() / valid.sum():.1%}")
    return out, best


stacks = {"red": red, "nir": nir}
best_out, best = best_pixel_composite(stacks, ndvi(red, nir), ["red", "nir"])

  61 of 127 scenes contribute pixels
  largest single contribution: 29.8%

A per-band median mixes bands from different dates in one pixel, which can produce a spectrum no real surface ever had — a red from July and a near-infrared from October. A best-pixel composite keeps each pixel spectrally coherent, at the cost of visible seams where the source date changes.

Choose per-band median for index maps and best-pixel for anything that feeds a classifier trained on real spectra.

Example 3 — a seasonal composite, which is usually what you meant

import numpy as np
from collections import defaultdict


def seasonal_composites(scenes, dates, season_of, min_obs=3):
    """One composite per season rather than one for the whole archive."""
    groups = defaultdict(list)
    for scene, date in zip(scenes, dates):
        groups[season_of(date)].append(scene)

    out = {}
    for season, members in sorted(groups.items()):
        stack = np.stack([m["ndvi"] for m in members])
        n_obs = np.isfinite(stack).sum(axis=0)
        with np.errstate(all="ignore"):
            values = np.nanmedian(stack, axis=0)
        out[season] = np.where(n_obs >= min_obs, values, np.nan)
        print(f"  {season}: {len(members):3d} scenes, "
              f"median NDVI {np.nanmedian(out[season]):.3f}, "
              f"{(n_obs < min_obs).mean():.1%} below min_obs")
    return out


SEASON = {6: "summer", 7: "summer", 8: "summer", 9: "autumn", 10: "autumn", 11: "autumn", 12: "winter"}
by_season = seasonal_composites([{"ndvi": ndvi(r, n)} for r, n in zip(red_scenes, nir_scenes)],
                                dates, lambda d: SEASON[int(d[5:7])])

  autumn:  54 scenes, median NDVI 0.556, 10.7% below min_obs
  summer:  55 scenes, median NDVI 0.725, 5.4% below min_obs
  winter:  18 scenes, median NDVI 0.475, 72.1% below min_obs

Splitting by season recovers the signal that a seven-month composite averages away: 0.725 in summer against 0.475 in December, where the single all-period composite reported 0.667.

It also shows the trade-off immediately. Winter — December alone here — has 18 scenes and a median of two clear looks per pixel, so 72.1% of pixels fall below a three-observation minimum. The winter number is real for the 28% of the map that has support and absent for the rest — which is exactly what you want a composite to admit.

Median beats mean, and composites have no date

Why the median beats the mean here

Both were computed over the same masked observations, and they differ by 0.12 NDVI. That gap is contamination: cloud-edge and thin-cirrus pixels that the mask did not catch.

Those survivors are not symmetric. Cloud is bright in both red and near-infrared, and brighter in red relative to a vegetated surface, so contaminated pixels pull NDVI down. A mean incorporates them proportionally; a median ignores anything that is not near the middle.

Given that no cloud mask is perfect, the median is the safer default. The cost is that it discards information: with only 3 or 4 clear observations, a median is a much noisier estimator than a mean would be.

Why a composite has no date

This is the most commonly misused property of composites. The output raster looks exactly like a scene — same grid, same bands, same file format — and it is not one. Each pixel is a statistic over a different subset of dates, and neighbouring pixels can be summarising different months.

Consequences worth stating in the metadata:

• It cannot date an event. A change visible in a composite happened somewhere in the window.
• Neighbouring pixels are not synchronous. A boundary in the output can be a boundary in observation dates rather than on the ground.
• Two composites over different windows are not directly comparable unless both windows are stated and similar.

Why every pixel got at least four looks

With a median usable fraction of 0.2% per scene, you might expect large permanently hidden areas. The measurement says otherwise: 904 pixels in a million were never seen clear, and 98.9% had at least three clear looks. The 904 are scattered, not a block — every one of their 127 classifications was cloud, cirrus, cloud shadow or dark area.

Cloud moves. Over seven months, the probability that a given pixel is under cloud on every one of 127 overpasses is small even when each individual overpass is mostly cloudy. That is the whole statistical argument for compositing, and it is why "how cloudy is this region" is a much less useful question than "how many overpasses does this region get".

Why to keep the observation count

The composite and the count are one product. A user who wants a conservative map masks the composite where n_obs is low; a user who wants full coverage keeps it and knows which parts are thin.

Distributing the composite without the count discards the only information about its own reliability.

All-NaN columns, memory and mixed dates

• np.nanmedian over an all-NaN column warns and returns NaN. Wrap stack reductions in np.errstate(all="ignore") and use the observation count to mask.
• Memory grows with the stack. 127 scenes of one 1000 × 1000 band as float32 is 508 MB; at 10 m over a full tile it is hundreds of gigabytes. Composite in tiles or chunk with dask.
• A per-band median mixes dates within a pixel. Fine for indices, wrong for anything needing a real spectrum.
• Composite the bands, then compute the index — not the other way round — if you want the output to be a plausible spectrum.
• Seasonality dominates long windows. A June-to-December composite is not "the landscape"; it is a seven-month average that includes senescence.
• Sensor differences matter. Sentinel-2A, 2B and 2C have slightly different spectral responses; mixing them in a composite adds a small step.
• Do not fill remaining gaps by interpolation without labelling it. An interpolated composite pixel is a model output, not an observation.
• Record the mask rule, the window and the reducer with the output. Without them the composite is uninterpretable.

Internal links

• How to mask clouds in Sentinel-2 imagery with Python — the step that makes this work
• Cloud masking explained — why the median scene is 0.2% usable
• My composite has holes, stripes or grey patches — diagnosing a composite that went wrong
• How to calculate NDVI from Sentinel-2 in Python — the index this composite feeds
• How to extract a vegetation index time series for a polygon — the alternative to compositing away time
• How to load Sentinel-2 bands into Python as an analysis-ready array — getting every scene on one grid
• How to turn a STAC search into an xarray data cube — the labelled-array version of the stack
• How to detect change between two satellite images in Python — why a dateless composite cannot do this
• How to make a Sentinel-1 composite over a period — the radar equivalent, where cloud is not the problem

FAQ

How do I make a cloud-free composite in Python?

Mask each scene with its classification band, put every scene on one grid, stack them, and take a per-pixel np.nanmedian over time. Return the per-pixel observation count with it.

Median or mean?

Median. On the stack measured here the masked mean sat 0.12 NDVI below the masked median, because residual cloud edges pull the mean down and the median ignores them.

How many scenes do I need?

Enough that every pixel gets several clear looks. Around Snowdon, 127 scenes across seven months gave the median pixel 14 and 98.9% of pixels at least three — but the median individual scene was only 0.2% usable.

Should I filter out cloudy scenes before compositing?

Only extremely cloudy ones. A 30%-usable scene still contributes those pixels, and in a cloudy region you cannot afford to discard them.

Why is my composite darker or greener than a single date?

Because it summarises a period. The best single scene here gave an NDVI median of 0.765 on 12 July; the June-to-December composite gave 0.667, including autumn senescence.

Can I use a composite to detect change?

Not to date it. A composite has no date — each pixel summarises a different subset of the window. Use two composites over two clearly stated periods, and say what the periods were.

What is a best-pixel composite?

One that picks a single whole observation per pixel by some score, instead of mixing bands from different dates. It keeps each pixel spectrally coherent at the cost of visible seams.