Storing Multi-Temporal Index Stacks as Zarr

One-sentence answer: stack the season into a (time, y, x) float32 array, chunk time coarsely and space at 256 pixels, write with consolidated=True and a Zstandard compressor with byte shuffle — then read a pixel’s whole series in one chunk fetch instead of thirty file opens.

Context

Thirty separate index rasters answer “show me 12 June” perfectly and “compute the seasonal integral per pixel” terribly, because the second question opens thirty remote files for every spatial block it touches. A Zarr store inverts that trade: the values for one pixel across all dates sit contiguously, so a per-pixel reduction reads a handful of objects.

Both layouts have their place, which is why the parent topic recommends keeping the per-date cloud-optimized GeoTIFFs and a Zarr stack. This guide is the second one, and its inputs are the accepted windows from satellite imagery APIs and STAC catalogues.

Chunking decides which question is cheap Left: a stack chunked one date per chunk, where reading one pixel's thirty-date series touches thirty separate objects. Right: the same stack with the whole time axis in one chunk, where the same series is a single object fetch, while reading a single date now pulls the whole depth of that spatial block. chunks = (1, 256, 256) — one date per chunk … 30 chunks one pixel's season = 30 object fetches one date's map = 1 fetch chunks = (30, 256, 256) — whole season per chunk one chunk 30 dates × 256 × 256 × 4 B ≈ 7.9 MB one pixel's season = 1 object fetch one date's map = the same fetch, 30× the bytes Chunk for the question you ask most; keep the per-date COGs for the other one.

Prerequisites

Beyond the parent topic’s stack: xarray 2024., zarr 2.18., numcodecs 0.13.* and rioxarray 0.15.*. Every date must already be on one grid — same CRS, transform and shape — which is what the aligned reads in windowed COG reads over HTTP with rasterio produce.

Step-by-step

1. Stack the dates into a (time, y, x) array, sorted by time.

Chunk shape is the whole design decision A table of four chunk shapes for a float32 index stack, with the resulting chunk size in megabytes and the query each favours. One date per chunk suits single-date display, a whole season per chunk suits per-pixel time series, ten dates is a compromise, and very large spatial blocks serve neither. Chunk shape (time, y, x) Chunk size The question it makes cheap (1, 256, 256) one date per chunk 0.26 MB Show me one date (10, 256, 256) the default here 2.6 MB Both, at a small cost to each (30, 256, 256) a whole season per chunk 7.9 MB Per-pixel seasonal statistics (30, 1024, 1024) large spatial blocks 126 MB Neither — a small field reads mostly nothing

2. Chunk with time coarse and space at 256.

3. Write with consolidated metadata.

4. Append new dates as the season progresses.

5. Verify a per-pixel read.

PYTHON 
import numpy as np
import pandas as pd
import xarray as xr
from numcodecs import Blosc

TIME_CHUNK = 10          # dates per chunk: a compromise between series reads and cheap appends
SPACE_CHUNK = 256


def build_stack(arrays: dict[pd.Timestamp, np.ndarray], transform, crs, name: str = "ndvi"):
    """Assemble per-date arrays into one (time, y, x) DataArray on a common grid."""
    dates = sorted(arrays)
    shapes = {a.shape for a in arrays.values()}
    assert len(shapes) == 1, f"dates are not on one grid: {shapes}"
    assert all(arrays[d].dtype == np.float32 for d in dates), "store indices as float32"

    ny, nx = shapes.pop()
    ys = transform.f + transform.e * (np.arange(ny) + 0.5)
    xs = transform.c + transform.a * (np.arange(nx) + 0.5)

    da = xr.DataArray(
        np.stack([arrays[d] for d in dates]),
        dims=("time", "y", "x"),
        coords={"time": pd.DatetimeIndex(dates), "y": ys, "x": xs},
        name=name,
    )
    da.rio.write_crs(crs, inplace=True)
    assert da["time"].to_index().is_monotonic_increasing, "time axis is not sorted"
    return da


def write_stack(da: xr.DataArray, store: str) -> None:
    """Write the stack with a compressor and chunking suited to per-pixel time series."""
    chunked = da.chunk({"time": TIME_CHUNK, "y": SPACE_CHUNK, "x": SPACE_CHUNK})
    encoding = {da.name: {
        "compressor": Blosc(cname="zstd", clevel=5, shuffle=Blosc.SHUFFLE),
        "chunks": (TIME_CHUNK, SPACE_CHUNK, SPACE_CHUNK),
        "_FillValue": np.float32("nan"),
        "dtype": "float32",
    }}
    chunked.to_dataset().to_zarr(store, mode="w", encoding=encoding, consolidated=True)


def append_date(da_new: xr.DataArray, store: str) -> None:
    """Append one acquisition along time. The grid must match exactly."""
    existing = xr.open_zarr(store, consolidated=True)
    var = list(existing.data_vars)[0]
    assert da_new.shape[1:] == existing[var].shape[1:], (
        f"grid mismatch: new {da_new.shape[1:]} vs store {existing[var].shape[1:]}")
    assert da_new["time"].values[0] > existing["time"].values[-1], (
        "appending a date at or before the last one — rewrite the store instead of appending")
    da_new.to_dataset(name=var).to_zarr(store, append_dim="time", consolidated=True)

Inline verification — read one pixel’s whole series and check both the values and the cost:

PYTHON 
import time

ds = xr.open_zarr(STORE, consolidated=True)
start = time.perf_counter()
series = ds["ndvi"].isel(y=120, x=140).values
elapsed = time.perf_counter() - start

print(f"{len(series)} dates in {elapsed * 1000:.0f} ms; "
      f"peak {np.nanmax(series):.3f} on {ds['time'].values[int(np.nanargmax(series))]}")
assert np.isfinite(series).sum() >= 6, "fewer than six finite dates — mostly cloud, or a bad mask"
assert elapsed < 1.0, "per-pixel series read is slow — check the time chunking"
assert np.nanmax(series) <= 1.0 and np.nanmin(series) >= -1.0, "NDVI outside physical range"

Settings worth being deliberate about

Setting Default here Why it matters
Time chunk 10 dates A compromise: a whole-season chunk makes series reads cheapest but rewrites on every append, while one date per chunk makes appends free and series reads thirty times more expensive
Spatial chunk 256 px At 128 the metadata outweighs the data for a single field; at 1024 a small field reads mostly nothing
Compressor Zstandard level 5 with byte shuffle Shuffle groups the exponent bytes of neighbouring floats, which is most of the compression on smoothly varying index data
dtype float32 float64 doubles the store for precision an index does not have; int16 scaling saves more but adds a conversion every consumer must remember
Fill value NaN Propagates correctly through means and maxima without every consumer remembering to mask a sentinel
consolidated True, on write and append One metadata object instead of one per array; on object storage this alone is usually a tenfold difference in open time

Chunk size in bytes is the number to sanity-check: 10 dates × 256 × 256 × 4 bytes is about 2.6 MB, which sits in the range object stores serve efficiently. Chunks far below a megabyte turn a read into a request storm; chunks far above ten waste bandwidth on data the query did not want.

Gotchas and edge cases

  • Appending into a whole-axis time chunk rewrites that chunk. If the store grows one date at a time, a time chunk of ten keeps appends cheap while still making a series read ten times cheaper than thirty separate opens. Chunk the whole axis only when the season is complete and the store is written once.
Appending one acquisition without rewriting the season Five steps for appending a new date: open the store with consolidated metadata, assert the spatial grid matches exactly, assert the new timestamp is later than the last one, append along the time dimension, and re-consolidate the metadata so readers still open the store in one request. Open store consolidated metadata Check grid shape must match exactly Check time strictly increasing Append along the time axis Re-consolidate or readers list every array Forgetting the final step is invisible locally and costs a request per array on object storage.

  • Consolidated metadata is not automatic on append. Passing consolidated=True on the append keeps it current; forgetting it leaves readers to list every array’s metadata object, which on object storage is the slowest part of opening the store.

  • Zarr does not carry geospatial metadata natively. rioxarray’s CRS is stored as an attribute by convention, and a consumer that ignores the convention gets an unprojected array. Always write the CRS, and assert it on read.

  • A NaN fill value costs nothing and prevents an entire class of bug. Using a sentinel like −9999 in a float stack means every consumer must remember to mask it; NaN propagates correctly through means and maxima by default — the same argument made in fixing nodata and NaN propagation in band math.

  • One store per field group, not per field. A store for a single 40-hectare field at 10 m is 400 × 400 pixels — smaller than one spatial chunk, so the chunking does nothing and the metadata outweighs the data. Group neighbouring fields, or accept the per-date COGs alone for small operations.

  • Time must be strictly increasing. Two scenes with the same timestamp — the tile-overlap duplicate discussed in searching Sentinel-2 scenes — break the append assertion, which is the correct outcome: deduplicate before stacking.

This guide is part of Cloud-Optimized Storage for Field Imagery — see there for object-key design, compression and lifecycle rules.

Parent section Cloud-Optimized Storage for Field Imagery Related Converting GeoTIFFs to Cloud-Optimized GeoTIFF with rio-cogeo