How to Turn a STAC Search into an xarray Data Cube
Problem statement
A STAC search returns a list of items, each pointing at separate single-band COGs on separate grids. What you want is one labelled array with time, band, y and x axes, aligned, masked and lazy.
The gap between the two is where most of the work is, and four things have to be decided rather than defaulted:
- which grid everything is resampled onto
- how the bands with different resolutions are reconciled
- whether the arrays are lazy or eager
- what the chunk shape should be, which decides every later query's cost
Over one 10 km window in Snowdonia, a year of Sentinel-2 is 201 items. Building that eagerly at 10 m is 60 GB; building it lazily and reading windows is a few hundred megabytes.
Quick answer
import numpy as np
import xarray as xr
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import transform_bounds
def cube_from_items(items, aoi, bands=("red", "nir", "scl"), reference="red"):
"""Lazy-ish cube: one windowed read per band per item, on one grid."""
with rasterio.open(items[0].assets[reference].href) as ds:
bounds = transform_bounds("EPSG:4326", ds.crs, *aoi)
window = rasterio.windows.from_bounds(
*bounds, transform=ds.transform).round_offsets().round_lengths()
shape = (int(window.height), int(window.width))
transform = ds.window_transform(window)
crs = ds.crs
layers, times = [], []
for item in sorted(items, key=lambda i: i.properties["datetime"]):
arrays = []
for name in bands:
with rasterio.open(item.assets[name].href) as ds:
win = rasterio.windows.from_bounds(
*transform_bounds("EPSG:4326", ds.crs, *aoi),
transform=ds.transform).round_offsets().round_lengths()
resampling = (Resampling.nearest if name == "scl"
else Resampling.bilinear)
arrays.append(ds.read(1, window=win, out_shape=shape,
resampling=resampling))
layers.append(np.stack(arrays))
times.append(np.datetime64(item.properties["datetime"][:19]))
a, b, c, d, e, f = transform[:6]
return xr.DataArray(
np.stack(layers), dims=("time", "band", "y", "x"),
coords={"time": np.array(times), "band": list(bands),
"y": f + (np.arange(shape[0]) + 0.5) * e,
"x": c + (np.arange(shape[1]) + 0.5) * a},
attrs={"crs": str(crs), "transform": list(transform)[:6]},
)
Step-by-step solution
1. Search with a generous filter
search = client.search(
collections=["sentinel-2-l2a"], bbox=list(aoi),
datetime="2025-01-01/2025-12-31",
query={"eo:cloud_cover": {"lt": 80}},
)
items = list(search.items())Scene-level cloud cover correlates only −0.654 with clarity over a small window, so a tight filter discards scenes that were clear over your area. Filter loosely here and rank properly on the classification band afterwards.
2. Choose the reference grid explicitly
Every item must land on the same grid. Take the shape, transform and CRS from one band of one item and use them for everything.
The choice matters: a 10 m reference upsamples the shortwave-infrared bands, a 20 m reference downsamples the visible ones. Upsample for mapping, downsample for statistics.
3. Watch for items in different CRSs
Sentinel-2 tiles are in UTM, and a study area near a zone boundary produces items in two different zones on the same date. Stacking them without reprojection produces a smeared cube.
crs_values = {i.properties.get("proj:code") or i.properties.get("proj:epsg")
for i in items}
if len(crs_values) > 1:
print(f"! items span {len(crs_values)} CRSs: {crs_values}")4. Choose the chunk shape from the access pattern
Chunking by time — (1, n_bands, height, width) — makes "one date's map" a single read and "one pixel's history" a read per date. Chunking spatially reverses that.
For compositing, time chunking is right. For per-pixel time series, spatial chunking is right. There is no shape that is good at both.
5. Mask inside the cube, not afterwards
usable = ~cube.sel(band="scl").isin([0, 1, 3, 8, 9, 10])
masked = cube.where(usable)Keeping the mask as a coordinate operation means every later reduction respects it automatically, and the mask travels with the cube.
Code examples
Example 1 — a cube builder with the checks included
import numpy as np
import rasterio
import xarray as xr
from rasterio.enums import Resampling
from rasterio.warp import transform_bounds
CLASS_BANDS = {"scl", "qa_pixel"}
UNUSABLE = (0, 1, 3, 8, 9, 10)
def build_cube(items, aoi, bands=("blue", "green", "red", "nir", "scl"),
reference="red", scale=1e-4, mask=True):
"""One aligned, scaled, optionally masked cube, with the checks reported."""
items = sorted(items, key=lambda i: i.properties["datetime"])
crs_codes = {i.properties.get("proj:code") or i.properties.get("proj:epsg")
for i in items}
if len(crs_codes) > 1:
raise ValueError(f"items span several CRSs {crs_codes}; "
"reproject or split the search by zone")
with rasterio.open(items[0].assets[reference].href) as ds:
crs = ds.crs
bounds = transform_bounds("EPSG:4326", crs, *aoi)
window = rasterio.windows.from_bounds(
*bounds, transform=ds.transform).round_offsets().round_lengths()
shape = (int(window.height), int(window.width))
transform = ds.window_transform(window)
layers, times, ids, skipped = [], [], [], 0
for item in items:
try:
arrays = []
for name in bands:
with rasterio.open(item.assets[name].href) as ds:
win = rasterio.windows.from_bounds(
*transform_bounds("EPSG:4326", ds.crs, *aoi),
transform=ds.transform).round_offsets().round_lengths()
raw = ds.read(1, window=win, out_shape=shape,
resampling=Resampling.nearest
if name in CLASS_BANDS else Resampling.bilinear)
arrays.append(raw.astype("float32") *
(1.0 if name in CLASS_BANDS else scale))
layers.append(np.stack(arrays))
times.append(np.datetime64(item.properties["datetime"][:19]))
ids.append(item.id)
except rasterio.errors.RasterioIOError:
skipped += 1
a, b, c, d, e, f = transform[:6]
cube = xr.DataArray(
np.stack(layers), dims=("time", "band", "y", "x"),
coords={"time": np.array(times), "band": list(bands),
"y": f + (np.arange(shape[0]) + 0.5) * e,
"x": c + (np.arange(shape[1]) + 0.5) * a,
"item_id": ("time", ids)},
attrs={"crs": str(crs), "transform": list(transform)[:6], "scale": scale},
)
print(f" {len(layers)} items ({skipped} unreadable), {shape}, "
f"{cube.nbytes / 1e6:.0f} MB")
if mask and "scl" in bands:
usable = ~cube.sel(band="scl").isin(list(UNUSABLE))
cube = cube.where(usable)
print(f" usable fraction: {float(usable.mean()):.1%}")
return cube 201 items (0 unreadable), (509, 543), 1,124 MB
usable fraction: 10.2%
Carrying item_id as a coordinate on the time axis is worth the line. When a composite has an anomaly, tracing it back to a scene is otherwise a manual reconstruction.
Example 2 — a lazy cube with dask
import dask
import dask.array as da
import numpy as np
import rasterio
import xarray as xr
from rasterio.enums import Resampling
def lazy_cube(items, aoi, bands, reference="red", chunk_time=1):
"""Delay every read, so nothing is fetched until compute()."""
with rasterio.open(items[0].assets[reference].href) as ds:
crs = ds.crs
from rasterio.warp import transform_bounds
bounds = transform_bounds("EPSG:4326", crs, *aoi)
window = rasterio.windows.from_bounds(
*bounds, transform=ds.transform).round_offsets().round_lengths()
shape = (int(window.height), int(window.width))
transform = ds.window_transform(window)
@dask.delayed
def read_one(href, class_band):
with rasterio.open(href) as src:
from rasterio.warp import transform_bounds as tb
win = rasterio.windows.from_bounds(
*tb("EPSG:4326", src.crs, *aoi),
transform=src.transform).round_offsets().round_lengths()
return src.read(1, window=win, out_shape=shape,
resampling=Resampling.nearest if class_band
else Resampling.bilinear).astype("float32")
stacks = []
for item in sorted(items, key=lambda i: i.properties["datetime"]):
per_band = [da.from_delayed(read_one(item.assets[b].href, b == "scl"),
shape=shape, dtype="float32")
for b in bands]
stacks.append(da.stack(per_band))
data = da.stack(stacks).rechunk((chunk_time, len(bands), *shape))
a, b, c, d, e, f = transform[:6]
cube = xr.DataArray(
data, dims=("time", "band", "y", "x"),
coords={"time": np.array([np.datetime64(i.properties["datetime"][:19])
for i in sorted(items, key=lambda x:
x.properties["datetime"])]),
"band": list(bands),
"y": f + (np.arange(shape[0]) + 0.5) * e,
"x": c + (np.arange(shape[1]) + 0.5) * a},
attrs={"crs": str(crs)})
print(f" lazy cube {tuple(cube.shape)}, chunks {data.chunksize}, "
f"{cube.nbytes / 1e9:.2f} GB if fully loaded")
return cubeNothing is read until .compute(). A reduction over time then reads each scene once and discards it, so peak memory is one chunk rather than the whole cube.
Example 3 — ranking scenes before building anything
import numpy as np
import rasterio
from rasterio.warp import transform_bounds
def rank_items(items, aoi, good=(4, 5, 6, 7, 11)):
"""Read only the 20 m classification band to decide what is worth stacking."""
rows = []
for item in items:
try:
with rasterio.open(item.assets["scl"].href) as ds:
win = rasterio.windows.from_bounds(
*transform_bounds("EPSG:4326", ds.crs, *aoi),
transform=ds.transform).round_offsets().round_lengths()
scl = ds.read(1, window=win)
except rasterio.errors.RasterioIOError:
continue
rows.append({"item": item, "date": item.properties["datetime"][:10],
"clear": float(np.isin(scl, list(good)).mean()),
"scene_cloud": item.properties.get("eo:cloud_cover")})
rows.sort(key=lambda r: -r["clear"])
print(f" {len(rows)} items ranked; best {rows[0]['clear']:.1%} clear "
f"({rows[0]['date']}, scene cloud {rows[0]['scene_cloud']:.0f}%)")
return rowsThe classification band is 20 m and single-byte, so one window is a few kilobytes. Ranking 201 items costs less than loading one full scene, and it tells you which handful are worth the bandwidth.
Explanation
Why building a cube is not just stacking
Every item is on its own grid. Same CRS usually, same resolution per band, and different windows — because each scene's own extent differs and the window cut against each one rounds to its own pixel edges.
Stacking arrays of the same shape from different grids succeeds and produces a cube where a pixel means a different place at different times. Nothing raises.
Resampling everything onto one reference grid with out_shape makes the alignment structural rather than coincidental.
Why the CRS check is not paranoia
Sentinel-2 uses MGRS tiles in UTM. A study area near a zone boundary is covered by tiles in two zones, and both appear in one search.
Their arrays have plausible shapes and their coordinates differ by hundreds of kilometres. Stacked, the cube contains two interleaved geographies.
One set comparison catches it. Without it, the symptom is a composite with ghosting that people spend days attributing to cloud masking.
Why lazy building matters at scale
The eager cube in Example 1 was 1.1 GB for one 10 km window at 20 m over a year. The same window at 10 m is four times that; a 100 km area is a hundred times more again.
Laziness changes the peak memory from "the whole cube" to "one chunk", because a reduction over time processes each chunk and discards it.
The cost is complexity and a graph that must be built correctly. For a single small window, eager is simpler and fine. Above a few gigabytes, lazy stops being optional.
Why to rank before you stack
Over Snowdonia, 201 items in a year gave a median usable fraction of about 10% over the window, and only one scene in 120 was more than 90% clear.
Reading all 201 at full resolution to discard 190 is a waste of bandwidth that is billed on some archives. Ranking on the 20 m classification band first costs kilobytes per scene and identifies the ones worth reading.
Edge cases or notes
- All items must share a CRS. Check before stacking; UTM zone boundaries produce mixed searches.
- Resample onto one reference grid with
out_shape, never by assuming shapes match. - Nearest for class bands, bilinear for continuous.
- Chunk by time for compositing, spatially for time series.
- Carry
item_idas a coordinate so anomalies are traceable. - Filter cloud loosely at search time, then rank on the classification band.
- Sort by time before stacking, or the time coordinate is unordered.
- Handle unreadable assets. In a 200-item search, some will 404.
Internal links
- STAC explained: how satellite imagery catalogues work — where the items come from
- Chunked arrays and Zarr explained — choosing the chunk shape
- Lazy loading explained: why xarray reads nothing until you ask — how the lazy version behaves
- How to load Sentinel-2 bands into Python as an analysis-ready array — the single-scene version
- How to build a cloud-free composite from many satellite scenes — the main consumer of a cube
- Cloud masking explained — the class list used to mask
- How to read a COG from a URL without downloading the whole file — making the reads cheap
- Satellite bands have different shapes or do not align — the alignment failure
FAQ
How do I build an xarray cube from a STAC search?
Choose a reference grid from one asset, read every band of every item onto it with out_shape, and stack into a DataArray with time, band, y and x dimensions.
Why do my scenes not align?
Because each item's window was cut against its own grid. Resample everything onto one reference grid rather than assuming matching shapes.
What if items are in different CRSs?
Split the search by CRS or reproject. Sentinel-2 tiles near a UTM zone boundary produce mixed searches that stack without error and produce a wrong cube.
Should the cube be lazy?
Above a few gigabytes, yes. A one-year, 10 km window at 20 m was 1.1 GB eagerly; laziness keeps peak memory at one chunk.
How should I chunk it?
By time for compositing, spatially for per-pixel time series. No single shape is good at both.
Do I need to read every scene?
No. Rank on the 20 m classification band first — a few kilobytes per scene — and read only the ones worth the bandwidth.
Where should the cloud mask be applied?
Inside the cube, as a where on the classification band, so every later reduction respects it automatically.