Building a Mosaic from a STAC Search

Search, load onto a shared grid, mask, then reduce along time:

import odc.stac
from pystac_client import Client

items = Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=bbox, datetime="2026-06-01/2026-08-31",
    query={"eo:cloud_cover": {"lt": 40}}).item_collection()
cube = odc.stac.load(items, bands=["red", "green", "blue", "scl"], crs="EPSG:32633",
                     resolution=10, bbox=bbox, groupby="solar_day", chunks={})
clear = ~cube.scl.isin([3, 8, 9, 10])
mosaic = cube[["red", "green", "blue"]].where(clear).median("time")

A STAC search returns a list of overlapping scenes in different tiles and sometimes different projections; a mosaic needs one grid and one value per pixel. This page belongs to seamless mosaicking and edge blending in Satellite Processing Workflows & Index Pipelines.

From Search Results to One Grid

Four steps from item list to mosaic A STAC search returns items from several tiles and dates. Filtering and ordering selects the useful ones. odc-stac loads them all onto one shared grid as a time stack. Each time slice is cloud masked, and a reducer collapses the stack to a single mosaic. Search → filter → load → mask and reduce search items, many tiles filter and order cloud cover, date load onto grid one CRS, one res mask and reduce first-valid or median groupby solar_day merges tiles captured on the same pass before reduction.

The groupby="solar_day" argument is the detail that most often gets missed. Sentinel-2 delivers one item per 100 km tile, so a single overpass across a large AOI produces several items with the same date. Without grouping, each becomes its own time slice mostly filled with nodata; with grouping, tiles from the same pass are merged into one slice first, which halves memory use and makes first-valid ordering meaningful.

Environment & Setup

Package Version pin Used for
pystac-client >=0.7 Searching the catalogue
odc-stac >=0.3.9 Loading items onto a shared grid
xarray >=2023.1 Masking and reducing
dask >=2023.1 Lazy loading of large stacks
rioxarray >=0.15 Writing the COG 
pip install "pystac-client>=0.7" "odc-stac>=0.3.9" "xarray>=2023.1" "dask>=2023.1" "rioxarray>=0.15"

Complete Working Example

import odc.stac
import rioxarray  # noqa: F401
import xarray as xr
from pystac_client import Client

SCL_BAD = [0, 1, 3, 8, 9, 10]


def search_items(bbox, dates, max_cloud=50):
    items = Client.open("https://earth-search.aws.element84.com/v1").search(
        collections=["sentinel-2-l2a"], bbox=bbox, datetime=dates,
        query={"eo:cloud_cover": {"lt": max_cloud}}).item_collection()
    # least cloudy first, so first-valid prefers clean scenes
    return sorted(items, key=lambda it: it.properties["eo:cloud_cover"])


def build_mosaic(bbox, dates, crs, res=10, rule="median", out="mosaic.tif"):
    items = search_items(bbox, dates)
    cube = odc.stac.load(items, bands=["red", "green", "blue", "nir", "scl"], crs=crs,
                         resolution=res, bbox=bbox, groupby="solar_day",
                         chunks={"x": 2048, "y": 2048})
    clear = ~cube.scl.isin(SCL_BAD)
    bands = cube[["red", "green", "blue", "nir"]].where(clear)
    if rule == "median":
        mos = bands.median("time", skipna=True)
    else:  # first-valid in item order
        mos = bands.bfill("time").isel(time=0)
    count = clear.sum("time").astype("uint8").rename("clear_count")
    result = xr.merge([mos, count])
    result.to_array("band").rio.to_raster(out, driver="COG", compress="DEFLATE", tiled=True)
    return result

The clear-observation count is written alongside the bands because a mosaic that looks clean can rest on one observation in some places and thirty in others; downstream analysis should know which is which. Loading with odc-stac is covered in more depth in loading STAC items into an xarray cube with stackstac, which compares the two loaders.

Ordering Items for First-Valid Mosaics

First-valid takes the first clear pixel in stack order Three scenes are stacked in order of increasing cloud cover. For each pixel, the first-valid rule takes the value from the first scene in the stack that is clear there. Scene one fills most of the mosaic, scene two fills holes under scene one's clouds, and scene three fills what remains. Stack order decides the winner scene 1 — 5% cloud scene 2 — 20% cloud scene 3 — 35% cloud first-valid mosaic mostly scene 1, hole filled by scene 2 Sort by cloud cover, by date nearness, or by a quality score — the sort is the policy.

First-valid produces a mosaic built from real, individual observations, which keeps spectral relationships intact and is the right choice for visual base maps and for anything that needs a single acquisition date per pixel. Median smooths residual cloud and noise but mixes dates. Choosing a pixel-selection rule for composites compares the options in detail.

Keeping Memory in Check

A summer of Sentinel-2 over a 100 km square is several hundred gigabytes at full resolution, so the load must stay lazy until the reduction. Chunking along x and y but not time keeps each chunk’s full time series together, which is what median and first-valid need. For very large areas, running the mosaic tile by tile over a fixed grid and writing each tile separately, then assembling with a VRT, is more robust than one enormous Dask graph.

Choosing the Date Window

The date window is the main lever on mosaic quality, and it pulls in two directions. A wider window gives every pixel more chances at a clear observation, which drives nodata towards zero and lets the median reject more residual cloud. It also mixes phenological states: a window from April to September blends bare spring fields with full summer canopy, and the median of those is a pixel that never existed. For a visual base map this is usually acceptable. For anything that feeds a model or an index, keep the window inside one phenological stage — typically one to three months in temperate regions — and accept a slightly higher nodata fraction, filling the remainder from an adjacent window only where necessary.

A practical pattern is a primary window plus a fallback: build the mosaic from the preferred window, then fill any pixel whose clear count is zero from a wider window, and record in a separate band which window supplied it. That keeps most of the mosaic seasonally coherent while still covering the persistently cloudy corners.

Handling Mixed Projections

Large AOIs straddle UTM zones, and items from different zones arrive in different CRSs. odc-stac reprojects everything onto the requested grid during load, so the mosaic is always consistent, but the choice of target CRS matters. For an AOI mostly inside one zone, use that zone; for a continental mosaic, an equal-area projection avoids distorting pixel areas, which matters for any later area statistics.

Verification

Three checks on a finished mosaic The nodata fraction should be near zero inside the AOI. The distribution of clear observation counts shows where the mosaic rests on few scenes. A visual inspection at tile boundaries reveals seams from differing acquisition dates. Check the mosaic before publishing it nodata inside AOI should be ~0% clear count flag pixels with 1–2 obs seams inspect tile boundaries Low clear counts predict where residual cloud and seams will appear. 
import numpy as np

m = build_mosaic(bbox, "2026-06-01/2026-08-31", "EPSG:32633")
nodata = float(m.red.isnull().mean())
counts = m.clear_count.values
print(f"nodata {nodata:.2%}, median clear obs {np.median(counts):.0f}, "
      f"pixels with <=2 obs {np.mean(counts <= 2):.1%}")
assert nodata < 0.01

Common Errors

The mosaic is mostly empty

Items were not grouped by solar day, so each slice is one tile. Add groupby="solar_day".

Memory explodes during load

Chunks are missing or chunk along time. Pass chunks={"x": 2048, "y": 2048}.

Visible seams between tiles

Neighbouring areas come from different dates. Use a median, or feather and match histograms as in histogram matching across scenes.

Values are scaled integers with an offset

Newer processing baselines add an offset to L2A values. Apply the scale and offset before reducing.

Frequently Asked Questions

Q: Why group items by solar day? Because tiled products deliver one item per tile, so a single overpass appears as several items with the same date. Grouping merges them into one time slice, which saves memory and makes first-valid ordering work as intended.

Q: Should I use median or first-valid? Median for analysis where residual cloud and noise matter more than acquisition date; first-valid for visual base maps and anything needing a single real observation per pixel.

Q: How many scenes do I need for a clean mosaic? It depends on cloudiness. In temperate summers five to ten clear-ish scenes are usually enough; in cloudy tropics you may need a whole season. The clear-count layer tells you whether you had enough.

Q: Can I build the mosaic without downloading scenes? Yes. odc-stac reads only the windows needed from cloud-optimised GeoTIFFs over HTTP, so nothing is downloaded in full.