Chunked array storage

Deep dive

How it works

Parts:

• chunk grid — shape S cut into boxes of shape C, ceil(S/C) per axis; edge boxes padded
• codec pipeline — per chunk: optional filters (shuffle, delta), then a compressor (LZ4, Zstandard, deflate)
• store — key → bytes: files, object keys or byte ranges in one file; Zarr v3 key of chunk (1, 0) is c/1/0
• metadata — shape, chunk shape, dtype, codecs, fill_value (Zarr: zarr.json)
• chunk index — coordinates → location; implicit in the key name, or a B-tree or array inside the file (HDF5)

Operations:

1. Create: fix shape, chunk shape, dtype, codecs, fill value in the metadata.
2. Write: per touched chunk: filter, compress, store under its key. A partial overwrite reads, merges and rewrites it.
3. Read slice: per axis, chunk range = lo/C ..= (hi-1)/C. For each chunk in the product: fetch, decode, copy the overlap.
4. Missing chunk: no key means every cell is the fill value (Zarr v2 spec), so all-fill chunks need not be stored.
5. Parallelism: chunks are independent: reads fan out as concurrent requests; writers to different chunks need no coordination, writers to one chunk need external locking (Zarr docs).
6. Shard (optional): pack many chunks into one object followed (or preceded) by an index.

 6x6 array, 3x3 chunks         read rows 1-4, column 2:
 ┌─────┬─────┐                 touches chunks (0,0) and (1,0), decodes 18 cells for 4
 │ 0,0 │ 0,1 │   store: c/0/0  c/0/1  c/1/0  c/1/1   (compressed bytes each)
 ├─────┼─────┤   sharded: [chunk][chunk]...[index: (offset, nbytes) per chunk]
 │ 1,0 │ 1,1 │
 └─────┴─────┘

Time-major chunks (long in time, small in space) make a point's time series one chunk and a map 128; map-major is the reverse; a balanced shape pays moderately for both.

Code

use std::collections::HashMap;

type Ix = [usize; 3]; // (t, y, x); row-major inside a chunk

/// Calls `f` for every index in the box lo..hi, last axis fastest.
fn walk(lo: Ix, hi: Ix, mut f: impl FnMut(Ix)) {
    for a in lo[0]..hi[0] { for b in lo[1]..hi[1] { for c in lo[2]..hi[2] { f([a, b, c]); } } }
}
fn lin(i: Ix, shape: Ix) -> usize { (i[0] * shape[1] + i[1]) * shape[2] + i[2] }

/// Delta filter, then run-length pairs (count, value).
fn encode(raw: &[u8]) -> Vec<u8> {
    let (mut prev, mut out) = (0u8, Vec::new());
    let delta: Vec<u8> = raw.iter().map(|&b| { let d = b.wrapping_sub(prev); prev = b; d }).collect();
    let mut i = 0;
    while i < delta.len() {
        let n = delta[i..].iter().take(255).take_while(|&&d| d == delta[i]).count();
        out.extend([n as u8, delta[i]]);
        i += n;
    }
    out
}
fn decode(enc: &[u8]) -> Vec<u8> {
    let (mut prev, mut out) = (0u8, Vec::new());
    for p in enc.chunks(2) { for _ in 0..p[0] { prev = prev.wrapping_add(p[1]); out.push(prev); } }
    out
}

#[derive(Debug, Default, PartialEq)]
struct Stats { chunks: usize, stored: usize, decoded: usize } // chunks overlapped, bytes fetched, bytes decoded

struct Array { shape: Ix, chunk: Ix, fill: u8, store: HashMap<Ix, Vec<u8>> } // store = the object store: key -> bytes

impl Array {
    fn from_dense(dense: &[u8], shape: Ix, chunk: Ix, fill: u8) -> Array {
        let mut store = HashMap::new();
        let grid = [0, 1, 2].map(|d| shape[d].div_ceil(chunk[d]));
        walk([0; 3], grid, |g| {
            let mut buf = vec![fill; chunk.iter().product()]; // edge chunks are padded with the fill value
            let hi = [0, 1, 2].map(|d| ((g[d] + 1) * chunk[d]).min(shape[d]));
            walk([0, 1, 2].map(|d| g[d] * chunk[d]), hi, |i| {
                buf[lin([0, 1, 2].map(|d| i[d] - g[d] * chunk[d]), chunk)] = dense[lin(i, shape)];
            });
            if buf.iter().any(|&b| b != fill) { store.insert(g, encode(&buf)); } // all-fill chunks are never stored
        });
        Array { shape, chunk, fill, store }
    }

    /// Reads the box lo..hi (hi exclusive): fetches and decodes only the chunks it overlaps.
    fn read_slice(&self, lo: Ix, hi: Ix) -> (Vec<u8>, Stats) {
        let out_shape = [0, 1, 2].map(|d| hi[d] - lo[d]);
        let mut out = vec![self.fill; out_shape.iter().product()];
        let mut st = Stats::default();
        let (g0, g1) = ([0, 1, 2].map(|d| lo[d] / self.chunk[d]), [0, 1, 2].map(|d| (hi[d] - 1) / self.chunk[d] + 1));
        walk(g0, g1, |g| {
            st.chunks += 1;
            let Some(enc) = self.store.get(&g) else { return }; // missing chunk = fill value, no I/O
            let raw = decode(enc);
            (st.stored, st.decoded) = (st.stored + enc.len(), st.decoded + raw.len());
            let a = [0, 1, 2].map(|d| lo[d].max(g[d] * self.chunk[d]));
            let b = [0, 1, 2].map(|d| hi[d].min((g[d] + 1) * self.chunk[d]).min(self.shape[d]));
            walk(a, b, |i| { // copy only the overlap
                out[lin([0, 1, 2].map(|d| i[d] - lo[d]), out_shape)] = raw[lin([0, 1, 2].map(|d| i[d] % self.chunk[d]), self.chunk)];
            });
        });
        (out, st)
    }
}

fn main() {
    let shape = [128, 128, 128]; // (time, y, x), one byte per cell
    let mut dense = vec![0u8; 128 * 128 * 128];
    walk([0; 3], shape, |i| { // blocky field, no data after t = 96
        if i[0] < 96 { dense[lin(i, shape)] = 1 + (i[0] / 16 * 7 + i[1] / 8 * 3 + i[2] / 8) as u8; }
    });
    let want = |lo: Ix, hi: Ix| { let mut v = vec![]; walk(lo, hi, |i| v.push(dense[lin(i, shape)])); v };
    let series = ([0, 40, 70], [128, 41, 71]); // one cell over all time
    let map = ([50, 0, 0], [51, 128, 128]); // the whole map at one time
    for (name, chunk) in [("time-major", [128, 16, 8]), ("map-major", [1, 128, 128]), ("balanced", [16, 32, 32])] {
        let a = Array::from_dense(&dense, shape, chunk, 0);
        for (what, (lo, hi)) in [("series", series), ("map", map)] {
            let (got, st) = a.read_slice(lo, hi);
            assert_eq!(got, want(lo, hi)); // identical to slicing the dense array
            println!("{name} {what:6}: {st:?}");
        }
        println!("{name} stored chunks: {} of {}", a.store.len(), shape.iter().zip(chunk).map(|(s, c)| s.div_ceil(c)).product::<usize>());
    }
    // ragged edges
    let (s, d): (Ix, Vec<u8>) = ([5, 7, 9], (0..315).map(|i| (i % 251) as u8 + 1).collect());
    let a = Array::from_dense(&d, s, [2, 3, 4], 0);
    assert_eq!(a.read_slice([0; 3], s).0, d);
    assert_eq!(a.read_slice([1, 2, 3], [5, 7, 9]).1.chunks, 3 * 3 * 3);
}

Exact results (series wants 128 cells, map 16,384; stored = compressed bytes fetched):

Map-major stores 96 of 128 chunks (no data after t = 96). Read amplification, decoded ÷ wanted cells: 1 for the map read on map-major, 128 on time-major, 12,288 for the series on map-major.

Complexity and constants

• Read: chunks overlapped per axis = (hi-1)/C - lo/C + 1, multiplied over axes. Decode work = chunks × chunk volume.
• Index: one entry per chunk; a Zarr v3 shard entry is two uint64 (16 bytes), so 4 chunks plus CRC32C is 68 bytes (shard spec).
• HDF5: a chunk holds at most 2^32-1 elements and 4 GB. netCDF-C default chunk size: 4 MiB.

Efficiency

• Cost: one slice = (chunks overlapped) requests or seeks plus whole-chunk decodes. Bound: request count and latency for small chunks, bytes decoded and bandwidth for large ones.
• Levers:
  • chunk shape matched to the query — 128 chunks → 1 (this page's code) — one dominant access pattern
  • chunks of 1 MB or more uncompressed — fewer requests (Zarr-Python docs) — object stores
  • sharding — 100 GB in 1 MB chunks: 100,000 files unsharded, 100 with 1 GB shards (Zarr-Python docs) — many small chunks
  • skip all-fill chunks — 96 of 128 stored here; Zarr docs: no objects for empty data — sparse data
  • chunk cache sized to a selection — repeats skip I/O (netCDF default: 10 chunks up to 64 MB) — overlapping reads
• Naive → standard → state of the art:
  • one contiguous file · a row is one read, a column one read per row
  • regular chunks, per-chunk compression, a key per chunk · reads proportional to the overlap
  • shards with an index, resolution levels, chunk shape per access pattern · few objects, range reads of inner chunks

Implementation notes

• Shape, then size: shape from the dominant read, then grow chunks past 1 MB uncompressed (netCDF guide: no default suits all uses).
• Codecs: Zarr v3 order: array→array (any), one array→bytes, bytes→bytes (any). Zarr v2's example: delta filter, then Blosc with LZ4 and shuffle.
• Zarr v3 sharding: inner chunk shape must divide the shard shape; index at start or end, recommended little-endian uint64 plus CRC32C; an empty inner chunk has offset and nbytes both 2^64-1 and reads as fill value. Zarr-Python orders inner chunks by Morton code by default (space-filling curves).
• HDF5 index: version 1 B-tree (node type 1 points at raw data chunks); the format also has single-chunk, implicit, fixed-array, extensible-array and version 2 B-tree indexes. Each open dataset has its own chunk cache; slot count: a prime, about 100 times the chunks that fit.
• COG: tiled GeoTIFF with overviews (image pyramids) as extra IFDs. GDAL's COG driver defaults: 512-pixel tiles, LZW. OGC maps tiles to the 2D tile matrix set (map tile pyramids).

Pitfalls

• Wrong shape: the series read on map-major chunks decodes 12,288 cells per cell wanted.
• Tiny chunks: more lookups, more I/O, bigger metadata (HDF5 docs); millions of objects stress stores. Shard.
• Huge chunks: a one-cell read decodes the whole chunk.
• Two writers on one chunk race on read-modify-write: align writes to chunk (or shard) boundaries. Shard writes come in whole-shard steps (Zarr-Python docs).
• Fill value hides absence: a fill of 0 or NaN looks like real 0 or NaN. Zarr v2 fill_value: null leaves absent chunks undefined; edge padding is undefined too.
• Rechunking is a separate full job: Zarr-Python has no in-place change; copy into a new array. Decide the shape early.

Variants and alternatives

In production

• netCDF-4 files use HDF5 as their base format, with chunking and per-variable compression.
• OME-Zarr (OME-NGFF 0.5): each image resolution level is a Zarr array.
• Cloud-Optimized GeoTIFF: an OGC standard (21-026) for rasters read by HTTP range requests, also from object storage.

Implementations

Sources

• Zarr specification v3.1, core. spec
• Zarr specification v2. spec
• Zarr sharding codec, ZEP 0002. spec
• Zarr-Python user guide, Optimizing performance. docs
• The HDF Group, Chunking in HDF5. docs
• The HDF Group, HDF5 File Format Specification Version 3.0. spec
• Unidata, NetCDF Users Guide: Improving Performance with Chunking. docs
• OGC, OGC Cloud Optimized GeoTIFF Standard, 21-026. spec · cogeo.org
• GDAL, COG -- Cloud Optimized GeoTIFF generator. docs
• OME-NGFF specification 0.5. spec