stackstac vs odc-stac for STAC-to-Array

Once a STAC search returns an ItemCollection, you still have to turn a list of JSON items pointing at Cloud-Optimized GeoTIFFs into a single analysis-ready array. Two libraries own this step: stackstac, which collapses everything into one 4D DataArray, and odc-stac, which builds an xarray Dataset with per-band control. They solve the same problem with different philosophies — uniform-and-convenient versus explicit-and-flexible — and the right pick depends on whether your assets share a grid. This page compares them head to head on defaults, reprojection control, mixed-resolution handling, grouping, and speed. For where cube-building sits in the wider raster and STAC picture, see Core Raster Fundamentals & STAC Mapping.

Both produce lazy Dask-backed arrays, so neither reads a pixel until you compute; the difference is in how much you have to spell out and how gracefully each handles a messy collection.

---

Decision framework

The comparison matrix below maps the collection you have onto the loader that fits it. Read across your row before writing the load call.

The pivotal question is homogeneity. If every item shares a CRS and resolution — a single Sentinel-2 tile ID across a season — stackstac’s one-line stack is the least friction. The moment you cross UTM zones or need different resampling for a categorical mask than for reflectance, odc-stac’s explicit per-band model stops fighting you.

---

Prerequisites

Install both loaders alongside a STAC client. odc-stac pulls in odc-geo for its grid model; stackstac depends on rasterio and xarray directly:

pip install "pystac-client>=0.7" "stackstac>=0.5" "odc-stac>=0.3" \
            "xarray>=2023.1" "rioxarray>=0.15" "dask[array]>=2023.1"

Conceptual prerequisites: you should already be able to build a filtered search with pystac-client and know how pixel resolution and scaling forces a target grid choice when bands differ. The examples assume the STAC assets are COGs so range reads work.

---

Head-to-head

Neither is strictly better. stackstac optimises for the case where the collection is already tidy and you want a cube in one call; odc-stac optimises for the case where the collection is not tidy and you need to say exactly how each band lands on the output grid.

---

When to choose stackstac

Choose stackstac when your items come from one tile ID (or otherwise share a CRS and native resolution) and you want a single 4D array to slice by time and band. Its mental model is minimal: give it the items, an EPSG code, and a resolution, and it returns one DataArray you can reduce over the time dimension immediately.

The benchmark below stacks a season of Sentinel-2 items, requests the visible and NIR bands at 10 m in a fixed CRS, and times a median composite over a small window. Setting epsg explicitly is the key habit — letting it infer across mixed zones is where stackstac surprises people.

import time
import stackstac
from pystac_client import Client

catalog = Client.open("https://earth-search.aws.element84.com/v1")
search = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[11.0, 46.9, 11.3, 47.1],
    datetime="2024-06-01/2024-08-31",
    query={"eo:cloud_cover": {"lt": 20}},
)
items = search.item_collection()

# One 4D DataArray: pin the projection and grid explicitly
cube = stackstac.stack(
    items,
    assets=["blue", "green", "red", "nir"],
    epsg=32632,             # do not let it infer across UTM zones
    resolution=10,          # single resolution for every band
    chunksize=2048,
    fill_value=0,
)
print(cube.dims, cube.shape)   # ('time', 'band', 'y', 'x')

t0 = time.perf_counter()
median = cube.sel(band="nir").median(dim="time")   # still lazy
result = median.compute()
print(f"stackstac median: {time.perf_counter() - t0:.2f}s, "
      f"shape={result.shape}")

The appeal is the flatness: one array, one indexing model, sel(band=...) and median(dim="time") fall straight out of xarray. The constraint is symmetrical — because every band lives on one grid with one resampling, stackstac has no way to nearest-sample a categorical mask while bilinear-resampling reflectance. For a homogeneous reflectance stack that limitation never bites.

---

When to choose odc-stac

Choose odc-stac when the collection is heterogeneous — multiple UTM zones, bands at 10/20/60 m, or a categorical mask that must not be interpolated — or when you want grouping and fusing handled at load time. It returns a Dataset with one variable per band, and every knob is per band.

The benchmark below loads reflectance bands plus the Sentinel-2 SCL mask, resamples reflectance bilinearly but the SCL band with nearest, groups sub-scenes into one observation per solar day, and times the same style of composite. None of the per-band control below is expressible in a single stackstac call.

import time
import odc.stac
from pystac_client import Client

catalog = Client.open("https://earth-search.aws.element84.com/v1")
items = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[11.0, 46.9, 11.3, 47.1],
    datetime="2024-06-01/2024-08-31",
    query={"eo:cloud_cover": {"lt": 20}},
).item_collection()

# Dataset with per-band resolution, resampling, and solar-day grouping
ds = odc.stac.load(
    items,
    bands=["red", "nir", "scl"],
    crs="EPSG:32632",
    resolution=10,
    groupby="solar_day",            # fuse split scenes into one date
    resampling={"scl": "nearest", "*": "bilinear"},  # mask must not interpolate
    chunks={"x": 2048, "y": 2048},
)
print(list(ds.data_vars))          # ['red', 'nir', 'scl']

t0 = time.perf_counter()
median = ds["nir"].median(dim="time")   # lazy until compute
result = median.compute()
print(f"odc-stac median: {time.perf_counter() - t0:.2f}s, "
      f"shape={result.shape}")

The per-band resampling dict is the feature that earns odc-stac its place: nearest on scl keeps the class codes intact while bilinear smooths reflectance, all in one load. groupby="solar_day" folds the along-track scene splits that plague Sentinel-2 into a single date, which otherwise requires a manual xarray groupby. The cost is a slightly larger conceptual surface — you are declaring a grid rather than accepting one.

---

Chunking, grouping, and performance

The two loaders converge on Dask but expose chunking differently, and that difference shapes throughput more than the codec inside the COGs does. stackstac takes a single chunksize that applies square tiles across y and x, keeping every band and date on the same spatial grid; odc-stac takes a chunks dict so you can make the time axis a chunk dimension too, which helps when a temporal reduction dominates the graph. For a spatial reduction — clipping to a watershed and compositing — square spatial chunks aligned to the COG’s internal 512-pixel tiles minimise redundant range reads on both loaders, because a Dask task then maps cleanly onto whole GeoTIFF tiles rather than straddling tile boundaries.

Grouping is where the runtime cost can diverge sharply. odc-stac’s groupby="solar_day" fuses along-track scene splits during the load, so a 24-date search that STAC reports as 40 items collapses to 24 observations before any reduction runs. Doing the same with stackstac means all 40 slices enter the cube and you fold the duplicates yourself with an xarray groupby after the fact, which materialises more intermediate chunks. Neither is wrong, but if your workload is dominated by temporal composites over split scenes, odc-stac’s built-in fusion removes a stage from the graph.

A fair benchmark therefore measures the whole path — search, cube construction, and one representative reduction over a fixed spatial window — not just the constructor call, because both constructors are near-instant (they only build a lazy graph). Time the .compute(), watch peak worker memory, and hold the target CRS, resolution, and chunk size identical across both loaders so the only variable is the library. On a homogeneous single-tile collection the two land within a small margin; the gap opens on messy collections where odc-stac’s per-band control avoids the extra reprojection and grouping passes stackstac would need.

---

Verification & testing

Whatever you build, verify the cube’s geometry and laziness before computing over hundreds of scenes. A cube with the wrong CRS or an accidentally eager read wastes an entire cluster run:

import xarray as xr

def verify_cube(cube, expected_epsg: int, min_dates: int = 1) -> None:
    """Assert a STAC-derived cube is lazy, projected, and non-empty."""
    # stackstac returns a DataArray; odc-stac a Dataset — normalise
    sample = cube if isinstance(cube, xr.DataArray) else cube[list(cube.data_vars)[0]]

    # Must be Dask-backed (lazy), not already in memory
    assert sample.chunks is not None, "cube is not lazy — reads were eager"

    # CRS must match the requested projection
    epsg = sample.rio.crs.to_epsg() if sample.rio.crs else None
    assert epsg == expected_epsg, f"cube CRS EPSG:{epsg}, expected {expected_epsg}"

    # Time dimension must carry at least the expected observations
    n_time = sample.sizes.get("time", 0)
    assert n_time >= min_dates, f"only {n_time} dates in cube"
    print(f"cube OK: EPSG:{epsg}, {n_time} dates, chunks={sample.chunks[0][:1]}")

Call it right after stackstac.stack or odc.stac.load and before any .compute(). The chunks is not None assertion is the important one: it catches the mistake of accidentally materialising the whole cube, which turns a windowed benchmark into an out-of-memory crash.

---

Pitfalls & troubleshooting

stackstac silently spans UTM zones — if items straddle two zones and you omit epsg=, stackstac picks one and reprojects the rest, which can balloon the output grid and blur data. Always pass epsg= explicitly for multi-tile searches, or filter the search to one tile first.

odc-stac categorical mask interpolated — loading the SCL band without a per-band resampling override defaults to the collection’s method, which may be bilinear and corrupts integer class codes into fractional nonsense. Pass resampling={"scl": "nearest", "*": "bilinear"} so the mask stays discrete.

Duplicate timestamps from scene splits — Sentinel-2 scenes split along orbit tracks appear as two items with the same datetime. In stackstac they become two time slices you must reduce; in odc-stac set groupby="solar_day" to fuse them at load. Forgetting this double-counts pixels in a temporal median.

All-zero or all-nodata cube — a mismatched fill_value (stackstac) or an unread asset nodata (odc-stac) makes valid pixels indistinguishable from fill. Confirm the asset’s declared nodata and set fill_value to match, then mask before compositing.

Runaway memory on compute — the chunk size controls how many pixels each Dask task pulls from the COGs. A 2048 chunk over many bands and dates can exceed worker memory; drop to 1024 or subset the cube in time and space before .compute(). See scaling raster processing with Dask for chunk tuning, and paginating large STAC searches when the item list itself is huge.

---

FAQ

What is the main difference between stackstac and odc-stac?

stackstac produces a single 4D DataArray (time, band, y, x) and leans on stackstac.stack with an explicit epsg and resolution, which is convenient for uniform collections like a single Sentinel-2 tile. odc-stac produces an xarray Dataset with one variable per band and gives finer control over per-band resolution, resampling, and grouping through odc.stac.load, which handles mixed-resolution and multi-CRS collections more gracefully.

How do I handle Sentinel-2 bands at 10m, 20m, and 60m in one cube?

Both loaders resample every asset onto one output grid, so you must pick a target resolution. With stackstac, set resolution=10 and the reader resamples 20 m and 60 m bands up to 10 m. With odc-stac, pass resolution=10 and a per-band resampling dict so you can nearest-sample the SCL mask while bilinear-resampling reflectance bands, which stackstac cannot express per band.

Do stackstac and odc-stac read pixels immediately?

No. Both return a lazy, Dask-backed xarray object whose values are only fetched when you call .compute() or write to disk. This is what lets you build a cube spanning hundreds of scenes and then subset it in time or space before any bytes are read from the COGs, so the choice of chunk size drives how much each worker pulls per task.

---

Related

• Core Raster Fundamentals & STAC Mapping — the parent section covering STAC, COGs, and the raster grid model both loaders assume
• Querying STAC Catalogs Programmatically — build the filtered ItemCollection you feed to either loader
• Using pystac-client to Filter Sentinel-2 Imagery by Date — the search that precedes cube construction
• Paginating Large STAC Searches with pystac-client — keep the item list manageable before stacking
• Scaling Raster Processing with Dask — tune chunking and clusters for the cube both loaders return