⚙️ How-To Satellite imagery Raster & rasterio Satellite Imagery · step 12 of 16 Updated

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 120 Sentinel-2 scenes covering a 10 km window in Snowdonia between June and December 2025:

median usable fraction per scene   0.4%
best single scene                100.0%
scenes at least 50% usable          8 of 120

One scene in fifteen 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 120 scenes gave every pixel between 4 and 30 clear observations, median 14. 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

# stack: (time, 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 4, median 14, max 30

The masking is what makes this work. Measured on the same stack:

composite NDVI median
mean, unmasked 0.063
median, unmasked 0.042
mean, masked 0.506
median, masked 0.658
Four composite methods over the same 120 scenes producing NDVI medians of 0.063, 0.042, 0.506 and 0.658.
Masking moves the answer by 0.6 NDVI. The choice between mean and median moves it by another 0.15.

Step-by-step solution

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_fraction = np.isin(scl, [4, 5, 6, 7, 11]).mean()

Scene-level eo:cloud_cover correlated only −0.654 with this over the 120-scene stack, 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 

mean, masked     NDVI median 0.506   p5 0.150   p95 0.719
median, masked   NDVI median 0.658   p5 0.239   p95 0.816

The masked mean sits 0.15 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 4 observations and a pixel with 30 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 gave an NDVI median of 0.761; the June-to-December composite gave 0.658. 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.

Distribution of clear observations per pixel across 120 scenes: minimum 4, median 14, maximum 30, with no pixel unobserved.
The stack is not uniformly deep. A composite pixel built from 4 looks is not the same measurement as one built from 30.

Code examples

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
  median of 120 scenes, observations 4-30 (median 14), 0 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

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
  summer: 48 scenes, median NDVI 0.734,  7.4% below min_obs
  autumn: 54 scenes, median NDVI 0.535,  0.5% below min_obs
  winter: 18 scenes, median NDVI 0.420, 59.5% below min_obs

Splitting by season recovers the signal that a six-month composite averages away: 0.734 in summer against 0.420 in winter, where the single all-period composite reported 0.658.

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

Explanation

Why the median beats the mean here

Both were computed over the same masked observations, and they differ by 0.15 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 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.4% per scene, you might expect large permanently hidden areas. The measurement says otherwise: minimum 4 observations, no unobserved pixel.

Cloud moves. Over six months, the probability that a given pixel is under cloud on every one of 120 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.

A dated scene against a composite that summarises a period, listing what each can and cannot answer.
The output file looks identical. Only the metadata records which one you have.

Edge cases or notes

  • 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. 120 scenes of one 509 × 543 band as float32 is 133 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 six-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.

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.15 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. Over Snowdonia, 120 scenes across six months gave every pixel between 4 and 30 — but the median individual scene was only 0.4% 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.761 in July; the six-month composite gave 0.658, 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.