Loading STAC Items into an xarray Cube with stackstac

Declare the output grid explicitly and the rest is lazy:

import stackstac

cube = stackstac.stack(
    items,
    assets=["red", "nir", "scl"],
    epsg=32636, resolution=10,
    bounds=(760_000, 9_940_000, 800_000, 9_980_000),   # in the target CRS
    chunksize=2048,
)
# (time, band, y, x) — nothing has been read yet

Letting stackstac infer the grid from the items is the most common way to build a cube far larger than intended. This page belongs to querying STAC catalogs programmatically in Core Raster Fundamentals & STAC Mapping.

From Search Results to a Cube

Many footprints, one grid Each STAC item has its own footprint, CRS and native resolution. stackstac projects every asset onto a single declared grid — one CRS, one resolution, one bounding box — producing a four-dimensional DataArray of time, band, y and x. Areas an item does not cover are filled with NaN for that time step. Items in, aligned cube out item footprints, varying CRS and resolution stackstac.stack epsg · resolution · bounds one grid for every date and band (time, band, y, x), lazy

The declared grid is what makes the cube usable. Every date shares one transform, so a pixel at (y, x) is the same ground location at every time step, and arithmetic across time or bands needs no further alignment. The trade-offs between this library and its main alternative are set out in stackstac vs odc-stac for STAC-to-array.

Environment & Setup

Package Version pin Used for
stackstac >=0.5 Building the lazy cube from items
pystac-client >=0.7 Searching the catalog
dask[array] >=2024.1 Lazy chunked execution
rioxarray >=0.15 Writing results back out with CRS 
pip install "stackstac>=0.5" "pystac-client>=0.7" "dask[array]>=2024.1" "rioxarray>=0.15"

Complete Working Example

import numpy as np
import pystac_client
import stackstac
import xarray as xr

API = "https://earth-search.aws.element84.com/v1"
SCL_KEEP = [4, 5, 6, 11]            # vegetation, not-vegetated, water, snow


def search(bbox_wgs84, period: str, max_cloud: float = 40):
    client = pystac_client.Client.open(API)
    return client.search(collections=["sentinel-2-l2a"], bbox=bbox_wgs84,
                         datetime=period,
                         query={"eo:cloud_cover": {"lt": max_cloud}}).item_collection()


def ndvi_cube(items, *, epsg: int, bounds: tuple[float, float, float, float],
              resolution: int = 10, chunksize: int = 2048) -> xr.DataArray:
    cube = stackstac.stack(
        items,
        assets=["red", "nir", "scl"],
        epsg=epsg, resolution=resolution, bounds=bounds,
        chunksize=chunksize,
        rescale=False,                 # apply scale/offset explicitly, per scene
        fill_value=np.nan, dtype="float32",
    )

    scl = cube.sel(band="scl")
    clear = scl.isin(SCL_KEEP)

    # Per-scene radiometry: newer baselines carry an offset of -0.1
    offset = xr.where(cube["s2:processing_baseline"].astype(str) >= "04.00", -0.1, 0.0)
    red = cube.sel(band="red") * 1e-4 + offset
    nir = cube.sel(band="nir") * 1e-4 + offset

    ndvi = ((nir - red) / (nir + red)).where(clear)
    ndvi.name = "ndvi"
    return ndvi


if __name__ == "__main__":
    items = search([34.6, 0.25, 34.9, 0.55], "2026-04-01/2026-09-30")
    ndvi = ndvi_cube(items, epsg=32636,
                     bounds=(760_000, 9_940_000, 800_000, 9_980_000))
    print(ndvi)                                            # lazy
    seasonal = ndvi.median("time", skipna=True).compute()  # now reads

Turning off rescale and applying the radiometry explicitly is a deliberate choice for Sentinel-2. The processing-baseline offset differs between scenes in the same search, and a single global rescale cannot represent that; the per-scene expression above can, because stackstac exposes item properties as coordinates on the time dimension. The reasoning is set out in converting int16 reflectance to float safely.

How Big Is the Cube, Really?

Inferred bounds versus declared bounds A search for a 40 kilometre area of interest returns items from four MGRS tiles. With inferred bounds the cube covers all four tiles, about 200 by 200 kilometres, and a seasonal median reads roughly 90 gigabytes. With declared bounds matching the area of interest the cube is 40 by 40 kilometres and the same median reads under 4 gigabytes. Same search, two grids inferred bounds 4 tiles · ~90 GB read for a median declared bounds 40 km AOI · ~4 GB read Check cube.nbytes before computing anything — it is metadata, and it is free. If it is larger than you expected, the grid is wrong, not the data.

cube.nbytes and cube.shape are computed from metadata alone, so they are free to inspect and they answer the question that matters before anything is read. Make a habit of printing both after every stack call; an unexpectedly large number almost always means inferred bounds or a resolution that was not set.

Reducing before computing matters just as much. ndvi.median("time").compute() reads every pixel of every date once and returns one layer; ndvi.compute().median("time") materialises the whole cube in memory first. The lazy graph makes both look identical until one of them runs out of memory.

Grouping Split Acquisitions by Solar Day

A single satellite pass over an area of interest often produces several items — one per MGRS tile it intersects — all with nearly identical timestamps. Stacking them yields separate time steps that each cover part of the area and are NaN elsewhere, which is correct but awkward: a “date” in the cube is then only partly observed, and counting valid observations per pixel double-counts nothing but looks as if it does.

Grouping those items by solar day merges them into one complete observation per pass:

import numpy as np

# One observation per pass: first valid value across items from the same day
daily = ndvi.groupby(ndvi.time.dt.floor("D")).first(skipna=True)
daily = daily.rename({"floor": "time"})

first(skipna=True) takes, at each pixel, the first non-NaN value among that day’s items — which, since the tiles do not overlap much, is simply the value from whichever tile covered that pixel. Do this before any temporal statistics, so that seasonal percentiles and observation counts refer to passes rather than to catalog records. The downstream reductions are covered in temporal aggregation and time series analysis.

Verification

Inspect the metadata before the pixels Before computing, confirm the cube's CRS is the one requested, its size in bytes is what the area of interest implies, its time axis has the expected number of dates, and the requested bands are present. All four checks read only metadata. Four free checks CRS as requested nbytes matches the AOI time steps = items returned bands all requested present None of these touches a pixel; all of them prevent a very expensive mistake. 
cube = stackstac.stack(items, assets=["red", "nir", "scl"], epsg=32636,
                       resolution=10, bounds=(760_000, 9_940_000, 800_000, 9_980_000))
assert cube.rio.crs is None or cube.rio.crs.to_epsg() == 32636 or cube.attrs["crs"].endswith("32636")
print(cube.shape, f"{cube.nbytes / 1e9:.1f} GB if fully read")
assert cube.sizes["time"] == len(items)
assert set(cube.band.values) == {"red", "nir", "scl"}

The time-step count catches one quiet failure: items with identical timestamps are kept as separate steps, so a scene split across two MGRS tiles appears twice. That is usually wanted, and grouping by solar day afterwards merges them — but it should be a deliberate step rather than a surprise in a later median.

Common Errors

The cube is hundreds of gigabytes

Bounds were inferred from items spanning several tiles. Pass bounds for the area of interest.

ValueError: Cannot pick a common CRS

Items span several UTM zones and no epsg was given. Choose one explicitly.

Every value is NaN after masking

The SCL band was resampled bilinearly, so class codes became fractions. Stackstac uses nearest by default; check that nothing forced another method.

Indices jump at a date in the middle of the series

The processing-baseline offset was not applied per scene. Apply it from each item’s properties.

Frequently Asked Questions

Q: Does stackstac download the imagery? Not when you call stack — it returns a lazy Dask-backed DataArray built from the item metadata. Pixels are read only when a computation needs them, and then only the windows that computation touches.

Q: Why pass bounds and epsg explicitly? Because the default grid is derived from the items’ footprints, which for a search spanning several tiles or UTM zones can be very large and in an unexpected projection. Declaring the grid makes the cube’s size and alignment a choice rather than a consequence.

Q: How do I avoid reading more than I need? Select assets, clip bounds to the area of interest, filter by cloud cover in the search, and reduce over time before computing. Each of those cuts the bytes read, and together they are usually the difference between minutes and hours.

Q: Can I write the cube straight to Zarr? Yes, after choosing a chunking that suits how it will be read. Writing a stackstac cube to Zarr materialises it once, which turns repeated catalog reads into fast local-region reads, as described in writing an xarray datacube to Zarr.