How to Turn a STAC Search into an xarray Data Cube
On this page
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, 93 of them under 80% scene cloud. Five bands of those 93 at 10 m are 2.06 GB held eagerly; built lazily, a reduction reads one 22.1 MB chunk at a time (Example 2).
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]},
)What you get — from 1. Search with a generous filter
93 items
Search loosely, one reference grid, one CRS, mask in the cube
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())
print(len(items), "items")93 items
Scene-level cloud cover correlated only −0.754 with clarity over this window (Example 3), 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.
Building, ranking and lazily loading a Snowdon cube
The data used on this page
Every Sentinel-2 L2A item on Earth Search for 2025 over a 10 km window on Snowdon, 4.16–4.00° W, 53.03–53.12° N, with scene cloud cover under 80% — the search in step 1. The setup only opens the catalogue; step 1's code runs the search and every example reads from its items.
from pystac_client import Client
client = Client.open("https://earth-search.aws.element84.com/v1")
aoi = AOI = (-4.16, 53.03, -4.00, 53.12)
print(client.id, "- window", AOI)earth-search-aws - window (-4.16, 53.03, -4.0, 53.12)
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
cube = build_cube(items, AOI, reference="scl") # the 20 m grid 93 items (0 unreadable), (509, 543), 514 MB
usable fraction: 38.7%
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 cube
lazy = lazy_cube(items, AOI, ("blue", "green", "red", "nir", "scl")) # the 10 m grid, nothing read
print(f" one chunk: {np.prod(lazy.data.chunksize) * 4 / 1e6:.1f} MB") lazy cube (93, 5, 1017, 1087), chunks (1, 5, 1017, 1087), 2.06 GB if fully loaded
one chunk: 22.1 MB
Nothing 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 rows
import time
started = time.perf_counter()
ranked = rank_items(items, AOI)
clear = np.array([r["clear"] for r in ranked])
scene = np.array([r["scene_cloud"] for r in ranked])
print(f" {time.perf_counter() - started:.0f} s for all of them; median clear {np.median(clear):.1%}; "
f"{(clear > 0.9).sum()} more than 90% clear, {(clear > 0.5).sum()} more than half clear")
print(f" correlation of scene cloud cover with clear share over the window: {np.corrcoef(scene, clear)[0, 1]:+.3f}") 93 items ranked; best 99.4% clear (2025-07-12, scene cloud 0%)
29 s for all of them; median clear 20.3%; 15 more than 90% clear, 29 more than half clear
correlation of scene cloud cover with clear share over the window: -0.754
The classification band is 20 m and single-byte, so one window is a few kilobytes. Ranking the 93 items took 29 s here, less than loading one full scene, and it tells you which ones are worth the bandwidth.
Why a cube is not a stack of scenes
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 514 MB for one 10 km window at 20 m over a year. The same window at 10 m is four times that, 2.06 GB; 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, the 93 items gave a median clear share of 20.3% over the window; 15 were more than 90% clear and 29 more than half clear.
Reading all 93 at full resolution to keep a few dozen 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.
Mixed CRSs, resampling choice and item IDs
- 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 514 MB eagerly and 2.06 GB at 10 m; laziness keeps peak memory at one 22.1 MB 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.