Chunked Arrays and Zarr Explained

Problem statement

A satellite time series is naturally a four-dimensional array: time, band, y, x. Nothing in the raster world stores that well. A GeoTIFF is two dimensions plus bands; a directory of GeoTIFFs makes time an afterthought and every read a file open.

Zarr stores the array as it is — n-dimensional, chunked into blocks, with each chunk a separate object and the metadata in a small JSON file.

The consequence that matters: the chunk shape decides which queries are fast. A cube chunked for maps is slow for time series and vice versa, and the difference can be two orders of magnitude in the number of objects read.

Quick answer

import xarray as xr

cube = xr.open_zarr("cube.zarr")               # "s3://bucket/cube.zarr" works the same way
print(dict(cube.sizes))
print(cube["reflectance"].encoding["chunks"])

{'time': 127, 'band': 2, 'y': 1000, 'x': 1000}
(1, 1, 500, 500)

That chunking — one date, one band and a quarter of the scene per chunk — makes "give me one date's red band" four chunk reads and "give me one pixel's history" 127 chunk reads per band, each of which discards all but one value.

Chunk shape from the query, size, metadata, dask

1. Understand what a chunk is

A chunk is the unit of storage and the unit of transfer. Each is compressed independently and stored as one object with a name encoding its position — 0.3.1.2 for the chunk at index (0, 3, 1, 2).

Reading any element requires fetching, decompressing and discarding the rest of its chunk. That is the same principle as a COG tile, extended to n dimensions.

2. Choose the chunk shape from the query, not the data

For the Snowdon cube, shape (time=127, band=2, y=1000, x=1000) (Example 2):

Neither is wrong; they answer different questions. If both matter, either store two chunkings, or compromise — (32, 2, 256, 256) costs 16 reads for a map and 4 for a pixel history, rather than 1 and 127 or 256 and 1. Timed on a local disk, the time-chunked store returned a map in 7.8 ms and a pixel's history in 268.3 ms; the pixel-chunked one took 312.8 ms and 3.1 ms; the compromise 72.3 ms and 22.5 ms (Example 4). On object storage, where every chunk is a request, the gaps are wider.

3. Size the chunks sensibly

Aim for a few megabytes uncompressed. Too small and the metadata and per-object overhead dominate — a million tiny objects is slow on any store and expensive on some. Too large and every read pulls data you discard.

chunk_bytes = np.prod(chunk_shape) * dtype_size

(1, 2, 1000, 1000) of float32 is 8.0 MB — a reasonable single chunk. (127, 2, 64, 64) is 4.2 MB, also fine.

4. Understand that the metadata is tiny and separate

Zarr stores array shape, dtype, chunk shape, compression and attributes in small JSON documents. A client reads those first — a few kilobytes — and then knows exactly which chunk objects to request.

That is why opening a Zarr store is fast regardless of its size, and why xr.open_zarr returns instantly on a terabyte: it has read the metadata and nothing else.

5. Combine with dask for larger-than-memory work

cube = xr.open_zarr("cube.zarr", chunks={"time": 30})
median = cube["reflectance"].isel(y=slice(0, 500), x=slice(0, 500)).median("time").compute()
print(dict(median.sizes), f"{float(median.sel(band='nir').mean()):.3f}")

{'band': 2, 'y': 500, 'x': 500} 0.617

chunks= makes the arrays dask arrays, so operations build a graph and execute in chunk-sized pieces. Aligning the dask chunks with the storage chunks avoids reading each stored chunk more than once.

One Snowdon cube under two chunkings

The data used on this page

cube.zarr is built from every Sentinel-2 L2A scene in Earth Search over a 10 km window around Snowdon (EPSG:32630, 422,880–432,880 E, 5,875,440–5,885,440 N) from 1 June to 31 December 2025 — 127 scenes, the same stack as Lazy loading explained. It holds reflectance (time, band, y, x) with the red and near-infrared bands, and the scene classification scl. The examples rewrite it with two deliberate chunkings.

import os

import numpy as np
import pystac_client
import rasterio
import xarray as xr
from pyproj import Transformer
from rasterio.windows import from_bounds

BOUNDS = (422_880, 5_875_440, 432_880, 5_885_440)
if not os.path.exists("cube.zarr"):
    to_wgs84 = Transformer.from_crs(32630, 4326, always_xy=True)
    bbox = [*to_wgs84.transform(*BOUNDS[:2]), *to_wgs84.transform(*BOUNDS[2:])]
    items = sorted(pystac_client.Client.open("https://earth-search.aws.element84.com/v1").search(
        collections=["sentinel-2-l2a"], bbox=bbox, datetime="2025-06-01/2025-12-31").items(),
        key=lambda i: i.id)

    def window(href):
        with rasterio.open(href) as src:
            return src.read(1, window=from_bounds(*BOUNDS, transform=src.transform))

    os.makedirs("scenes", exist_ok=True)
    for n, item in enumerate(items):
        path = f"scenes/{item.id}.npz"
        if not os.path.exists(path):
            scl = window(item.assets["scl"].href).repeat(2, axis=0).repeat(2, axis=1)
            np.savez_compressed(path, red=window(item.assets["red"].href), nir=window(item.assets["nir"].href),
                                scl=scl, date=str(item.datetime.date()))
        z = np.load(path)
        ds = xr.Dataset(
            {"reflectance": (("time", "band", "y", "x"), np.stack([z["red"], z["nir"]])[None] / np.float32(10_000)),
             "scl": (("time", "y", "x"), z["scl"][None].astype("uint8"))},
            coords={"time": [np.datetime64(str(z["date"]))], "band": ["red", "nir"],
                    "y": 5_885_435 - 10 * np.arange(1000), "x": 422_885 + 10 * np.arange(1000)})
        ds.reflectance.encoding = {"scale_factor": np.float32(1e-4), "dtype": "uint16", "_FillValue": 0}
        ds.to_zarr("cube.zarr", mode="w" if n == 0 else "a", zarr_format=2,
                  **({} if n == 0 else {"append_dim": "time"}))
print(dict(xr.open_zarr("cube.zarr").sizes))

{'time': 127, 'band': 2, 'y': 1000, 'x': 1000}

Example 1 — writing a cube with a deliberate chunking

import numpy as np
import xarray as xr


def write_cube(arrays, times, bands, transform, crs, path,
               time_chunk=1, spatial_chunk=None):
    """Write a time-band-y-x cube with an explicit chunk shape."""
    stack = np.stack(arrays)                      # (time, band, y, x)
    n_time, n_band, height, width = stack.shape

    a, b, c, d, e, f = transform[:6]
    xs = c + (np.arange(width) + 0.5) * a
    ys = f + (np.arange(height) + 0.5) * e

    cube = xr.Dataset(
        {"reflectance": (("time", "band", "y", "x"), stack.astype("float32"))},
        coords={"time": np.array(times, dtype="datetime64[s]"),
                "band": list(bands), "y": ys, "x": xs},
        attrs={"crs": str(crs)},
    )

    chunks = (time_chunk, n_band,
              spatial_chunk or height, spatial_chunk or width)
    bytes_per_chunk = np.prod(chunks) * 4
    print(f"  cube {stack.shape}, chunks {chunks} "
          f"({bytes_per_chunk / 1e6:.1f} MB each, "
          f"{np.prod([int(np.ceil(s / c)) for s, c in zip(stack.shape, chunks)])} chunks)")
    if bytes_per_chunk < 1e6:
        print("  ! chunks under 1 MB — per-object overhead will dominate")
    if bytes_per_chunk > 100e6:
        print("  ! chunks over 100 MB — every read pulls a lot you discard")

    cube["reflectance"].encoding["chunks"] = chunks
    cube.to_zarr(path, mode="w", zarr_format=2)
    return cube


source = xr.open_zarr("cube.zarr")
arrays = list(source["reflectance"].values)                  # one (band, y, x) array per scene
transform = (10, 0, 422_880, 0, -10, 5_885_440)
write_cube(arrays, source.time.values, source.band.values, transform, "EPSG:32630", "by_time.zarr")
write_cube(arrays, source.time.values, source.band.values, transform, "EPSG:32630", "by_pixel.zarr",
           time_chunk=len(arrays), spatial_chunk=64)

  cube (127, 2, 1000, 1000), chunks (1, 2, 1000, 1000) (8.0 MB each, 127 chunks)
  cube (127, 2, 1000, 1000), chunks (127, 2, 64, 64) (4.2 MB each, 256 chunks)

Printing the chunk size and count at write time is the cheapest possible guard. A store with a million 4 kB chunks is unusable and looks fine until someone tries to read it.

Example 2 — measuring which access pattern your chunking suits

import numpy as np


def chunk_read_cost(shape, chunks, queries):
    """How many chunks does each query touch?"""
    print(f"  array {shape}, chunks {chunks}")
    for name, selection in queries.items():
        touched = 1
        for dim, (size, chunk) in enumerate(zip(shape, chunks)):
            span = selection.get(dim, size)
            touched *= int(np.ceil(span / chunk)) if span > 1 else 1
        useful = np.prod([selection.get(d, s) for d, s in enumerate(shape)])
        fetched = touched * np.prod(chunks)
        print(f"    {name:28} {touched:6,} chunks, "
              f"{fetched / max(useful, 1):8.1f}x amplification")


queries = {"one date, whole scene": {0: 1},
           "one pixel, all dates": {1: 1, 2: 1, 3: 1},
           "small window, all dates": {1: 1, 2: 64, 3: 64}}
for chunks in ((1, 2, 1000, 1000), (127, 2, 64, 64), (32, 2, 256, 256)):
    chunk_read_cost((127, 2, 1000, 1000), chunks, queries)

  array (127, 2, 1000, 1000), chunks (1, 2, 1000, 1000)
    one date, whole scene             1 chunks,      1.0x amplification
    one pixel, all dates            127 chunks, 2000000.0x amplification
    small window, all dates         127 chunks,    488.3x amplification
  array (127, 2, 1000, 1000), chunks (127, 2, 64, 64)
    one date, whole scene           256 chunks,    133.2x amplification
    one pixel, all dates              1 chunks,   8192.0x amplification
    small window, all dates           1 chunks,      2.0x amplification
  array (127, 2, 1000, 1000), chunks (32, 2, 256, 256)
    one date, whole scene            16 chunks,     33.6x amplification
    one pixel, all dates              4 chunks, 132104.1x amplification
    small window, all dates           4 chunks,     32.3x amplification

The amplification factor — bytes fetched divided by bytes used — is the number to design against. Above about 10 the chunking is wrong for that query.

Example 3 — rechunking an existing store

import xarray as xr


def rechunk(src_path, dst_path, chunks, variable=None):
    """Rewrite a store with a different chunk shape."""
    source = xr.open_zarr(src_path)
    names = [variable] if variable else list(source.data_vars)

    for name in names:
        old = source[name].encoding.get("chunks")
        source[name] = source[name].chunk(chunks)
        source[name].encoding = {}
        source[name].encoding["chunks"] = tuple(
            chunks.get(d, source[name].sizes[d]) for d in source[name].dims)
        print(f"  {name}: {old} -> {source[name].encoding['chunks']}")

    source.to_zarr(dst_path, mode="w", zarr_format=2)
    print(f"  written to {dst_path}")
    return dst_path


rechunk("by_time.zarr", "compromise.zarr", {"time": 32, "y": 256, "x": 256})

  reflectance: (1, 2, 1000, 1000) -> (32, 2, 256, 256)
  written to compromise.zarr

Rechunking reads everything and writes everything, so it is expensive — but it is a one-off, and for a store queried in a pattern its chunking does not suit, it pays back immediately.

For very large stores, a dedicated rechunking tool that stages through an intermediate is more memory-efficient than a straight read-write.

Example 4 — the two chunkings, timed on disk

import os
import time


def timed(fn, repeats=3):
    best = float("inf")
    for _ in range(repeats):
        start = time.perf_counter()
        fn()
        best = min(best, time.perf_counter() - start)
    return best * 1000


for store in ("by_time.zarr", "by_pixel.zarr", "compromise.zarr"):
    files = sum(len(names) for _, _, names in os.walk(os.path.join(store, "reflectance")))
    size = sum(os.path.getsize(os.path.join(d, n)) for d, _, names in os.walk(store) for n in names)
    c = xr.open_zarr(store)["reflectance"]
    one_map = timed(lambda: c.isel(time=60).values)
    one_pixel = timed(lambda: c.isel(y=500, x=500).values)
    print(f"{store:17} {files:5} chunk files, {size / 1e6:6.0f} MB   one date's map {one_map:7.1f} ms   "
          f"one pixel's history {one_pixel:7.1f} ms")

by_time.zarr        129 chunk files,    761 MB   one date's map     7.8 ms   one pixel's history   268.3 ms
by_pixel.zarr       258 chunk files,    774 MB   one date's map   312.8 ms   one pixel's history     3.1 ms
compromise.zarr      66 chunk files,    765 MB   one date's map    72.3 ms   one pixel's history    22.5 ms

Chunks as the design, and Zarr against NetCDF

Why chunking is the whole design

In a chunked format there is no index scan, no query planner and no server. The client computes which chunks a query touches from the array's shape and chunk shape, and requests exactly those objects.

That makes performance completely predictable — and completely determined by the chunk shape. There is no runtime optimisation to compensate for a bad choice.

The amplification factor above makes the stakes concrete: 2,000,000× for one band of one pixel's history against a time-chunked cube. The client fetches an 8 MB chunk to return four bytes, 127 times over — and on a local disk that took 268.3 ms against 3.1 ms from the pixel-chunked store (Example 4).

Why Zarr and not NetCDF or HDF5

All three store chunked n-dimensional arrays with metadata. The difference is the storage model.

NetCDF and HDF5 put everything in one file with an internal index, which suits a filesystem and needs a library that can seek within it. Reading them from object storage means either downloading the file or a complicated protocol over range requests.

Zarr puts each chunk in a separate object and the metadata in separate small documents. On object storage that is exactly right: every chunk is a key, and fetching one is a plain GET.

The trade is a very large number of objects, which some filesystems handle badly. Zarr v3's sharding addresses that by grouping chunks into larger objects with an internal index.

Why chunk size has a sweet spot

Too small: each chunk carries per-object overhead — an HTTP request, a metadata lookup, a compression header. At 4 kB chunks, the overhead exceeds the data.

Too large: every read fetches the whole chunk regardless of how little you wanted. At 500 MB chunks, reading a single pixel costs 500 MB.

A few megabytes balances the two: large enough that per-object overhead is negligible, small enough that a targeted read is targeted. That is also roughly the size at which HTTP transfer becomes bandwidth-bound rather than latency-bound.

Why to align dask chunks with storage chunks

Dask chunks describe the computation; Zarr chunks describe the storage. When they differ, each dask chunk reads parts of several stored chunks, and each stored chunk is read by several dask tasks.

Aligning them — or making the dask chunks whole multiples of the storage chunks — means each stored chunk is fetched exactly once. Misaligned by one, and a chunk can be fetched twice for no benefit.

Tiny chunks, rechunking, sharding and consolidated metadata

• The chunk shape decides the performance, and there is no runtime compensation.
• Aim for a few megabytes per chunk uncompressed.
• A million tiny chunks is slow everywhere and expensive on object storage.
• Align dask chunks with storage chunks, or read each one twice.
• Rechunking rewrites everything. Choose well the first time.
• Zarr v3 sharding groups chunks into larger objects with an internal index.
• Metadata is separate and tiny, which is why opening a huge store is instant.
• Consolidated metadata (zarr.consolidate_metadata) turns many small metadata reads into one.

Internal links

• Cloud-native geospatial explained — the pattern Zarr belongs to
• Lazy loading explained: why xarray reads nothing until you ask — how a cube stays cheap to open
• How to turn a STAC search into an xarray data cube — building one from imagery
• How to choose chunk and tile sizes that actually help — sizing chunks in practice
• Cloud-optimised GeoTIFF explained — the two-dimensional equivalent
• How to scale a GeoPandas job with dask-geopandas — the same chunking idea for tables
• How to process a very large GeoPackage in chunks with Python — chunking without a chunked format
• How to build a cloud-free composite from many satellite scenes — a natural cube operation
• GRIB, NetCDF or Zarr — choosing a gridded format

FAQ

What is Zarr?

A format for chunked, compressed, n-dimensional arrays where each chunk is a separate object and the metadata is a small JSON document — designed for object storage.

How do I choose a chunk shape?

From the query you will run most. Time-chunked cubes make maps cheap and pixel histories expensive; spatially chunked cubes do the reverse.

How big should a chunk be?

A few megabytes uncompressed. Much smaller and per-object overhead dominates; much larger and every read pulls data you discard.

Why is my pixel time series so slow?

Because the cube is chunked by time. Reading one pixel across 127 dates touches 127 chunks, each of which is a whole scene.

Should I use Zarr or NetCDF?

Zarr for object storage, because each chunk is an independently fetchable object. NetCDF and HDF5 for local filesystems, where their single-file model is an advantage.

Do I need dask to use Zarr?

No — xarray reads Zarr without it. Use dask when the computation is larger than memory, and align the dask chunks with the storage chunks.

Can I change the chunking later?

Yes, by rewriting the store. It reads and writes everything, so it is expensive but one-off.