Byte-reference stores

A multidim VRT mosaic over netCDF, the byte table behind it, and kerchunk-Parquet Zarr.

Many archives are thousands of netCDF or HDF5 files. Two ways to read such a collection as one array, without copying it:

• a multidimensional VRT (gdal mdim mosaic) stitches the files along a dimension, and GDAL reads each file through its own driver;
• a byte-reference table records, for every chunk of every array, which file holds it and at what byte offset and size; any reader that knows the codec can then fetch chunks directly. kerchunk’s Parquet references are one form of this; the table below is another.

Both are small text or tabular artifacts that describe data rather than holding it.

A mosaic of three days

Three consecutive daily files of NOAA OISST v2.1, as NOAA publishes them, and a multidim VRT that gdal mdim mosaic wrote over them:

Data: NOAA 0.25-degree Daily Optimum Interpolation Sea Surface Temperature (OISST), Version 2.1; doi:10.25921/RE9P-PT57; licence: CC0-1.0; days 20250101, 20250102, 20250103. Built by scripts/fixtures/make_oisst.py; details in docs/fixtures/oisst/provenance.json.

vrt = data("oisst/oisst.vrt")
mos = xr.open_dataset(vrt, engine="gdalxarray")
mos

<xarray.Dataset> Size: 100MB
Dimensions:  (time: 3, zlev: 1, lat: 720, lon: 1440)
Coordinates:
  * time     (time) float64 24B 1.717e+04 1.717e+04 1.717e+04
  * zlev     (zlev) float64 8B 0.0
  * lat      (lat) float64 6kB -89.88 -89.62 -89.38 -89.12 ... 89.38 89.62 89.88
  * lon      (lon) float64 12kB 0.125 0.375 0.625 0.875 ... 359.4 359.6 359.9
Data variables:
    sst      (time, zlev, lat, lon) float64 25MB ...
    anom     (time, zlev, lat, lon) float64 25MB ...
    err      (time, zlev, lat, lon) float64 25MB ...
    ice      (time, zlev, lat, lon) float64 25MB ...

xarray.Dataset

• Dimensions:
  • time: 3
  • zlev: 1
  • lat: 720
  • lon: 1440
• Coordinates: (4)
  • time
    (time)
    float64
    1.717e+04 1.717e+04 1.717e+04

    long_name :

    Center time of the day

    axis :

    T

    array([17167., 17168., 17169.])

  • zlev
    (zlev)
    float64
    0.0

    actual_range :

    0, 0

    long_name :

    Sea surface height

    positive :

    down

    axis :

    Z

    array([0.])

  • lat
    (lat)
    float64
    -89.88 -89.62 ... 89.62 89.88

    grids :

    Uniform grid from -89.875 to 89.875 by 0.25

    long_name :

    Latitude

    axis :

    Y

    array([-89.875, -89.625, -89.375, ...,  89.375,  89.625,  89.875], shape=(720,))

  • lon
    (lon)
    float64
    0.125 0.375 0.625 ... 359.6 359.9

    grids :

    Uniform grid from 0.125 to 359.875 by 0.25

    long_name :

    Longitude

    axis :

    X

    array([1.25000e-01, 3.75000e-01, 6.25000e-01, ..., 3.59375e+02, 3.59625e+02,
           3.59875e+02], shape=(1440,))

• Data variables: (4)
  • sst
    (time, zlev, lat, lon)
    float64
    ...

    long_name :

    Daily sea surface temperature

    valid_max :

    4500

    valid_min :

    -300

    units :

    Celsius

    [3110400 values with dtype=float64]

  • anom
    (time, zlev, lat, lon)
    float64
    ...

    long_name :

    Daily sea surface temperature anomalies

    valid_max :

    1200

    valid_min :

    -1200

    units :

    Celsius

    [3110400 values with dtype=float64]

  • err
    (time, zlev, lat, lon)
    float64
    ...

    long_name :

    Estimated error standard deviation of analysed_sst

    valid_max :

    1000

    valid_min :

    0

    units :

    Celsius

    [3110400 values with dtype=float64]

  • ice
    (time, zlev, lat, lon)
    float64
    ...

    long_name :

    Sea ice concentration

    valid_max :

    100

    valid_min :

    0

    units :

    %

    [3110400 values with dtype=float64]

fig, axs = plt.subplots(1, 3, figsize=(13, 2.8), sharey=True, layout="constrained")
for i, ax in enumerate(axs):
    mos["sst"].isel(time=i, zlev=0).plot(ax=ax, cmap="RdYlBu_r", vmin=-2, vmax=32,
                                          add_colorbar=(i == 2))
    ax.set_facecolor("0.85")
    ax.set_title(f"time index {i}")
plt.show()

Sea surface temperature on each day of the mosaic, each day from its own netCDF file. Grey is land and sea ice.

time arrives as plain numbers here because gdal mdim mosaic carries the values over but not their units, so xarray has nothing to decode them with. Opening one of the source files directly gives decoded times; the Virtual paths page does that.

The byte table

The same three files, described chunk by chunk. The table was read from the files with GDAL (MDArray.GetRawBlockInfo) when the fixture was built, and GDAL’s Parquet driver reads it back:

def read_table(path):
    ds = gdal.OpenEx(path, gdal.OF_VECTOR)   # keep the dataset alive while
    lyr = ds.GetLayer(0)                      # its layer is in use
    defn = lyr.GetLayerDefn()
    cols = [defn.GetFieldDefn(i).GetName() for i in range(defn.GetFieldCount())]
    return pd.DataFrame([{c: f.GetField(c) for c in cols} for f in lyr])

refs = read_table(data("oisst/refs_sst.parquet"))
refs

	time	zlev	lat	lon	present	path	offset	size	info
0	0	0	0	0	True	nc/oisst-avhrr-v02r01.20250101.nc	58888	667949	COMPRESSION=DEFLATE; FILTER=SHUFFLE; ENDIANNES...
1	1	0	0	0	True	nc/oisst-avhrr-v02r01.20250102.nc	58888	667794	COMPRESSION=DEFLATE; FILTER=SHUFFLE; ENDIANNES...
2	2	0	0	0	True	nc/oisst-avhrr-v02r01.20250103.nc	58888	668444	COMPRESSION=DEFLATE; FILTER=SHUFFLE; ENDIANNES...

One row per chunk: the chunk’s index along each dimension, the file it lives in, its byte offset and size, and how it is encoded, as GDAL reports it. OISST stores each day’s field as a single chunk:

tables = {t["file"]: t for t in data_json("oisst/provenance.json")["products"]["byte_tables"]}
chunks = tables["refs_sst.parquet"]["chunks"]
dict(zip(tables["refs_sst.parquet"]["dims"], chunks))

{'time': 1, 'zlev': 1, 'lat': 720, 'lon': 1440}

Decoding a chunk by hand

With a row of the table, a chunk is a range read and whatever the info column lists: here a zlib inflate and a byte unshuffle. No netCDF library is involved:

row = refs.iloc[1]
info = set(x.strip() for x in row["info"].split(";"))
print(sorted(info))

f = gdal.VSIFOpenL(data("oisst/" + row["path"]), "rb")
gdal.VSIFSeekL(f, int(row["offset"]), 0)
raw = gdal.VSIFReadL(1, int(row["size"]), f)
gdal.VSIFCloseL(f)

buf = zlib.decompress(raw) if "COMPRESSION=DEFLATE" in info else raw
itemsize = 2                                              # int16
if "FILTER=SHUFFLE" in info:                               # bytes grouped by position
    n = len(buf) // itemsize
    buf = np.frombuffer(buf, np.uint8).reshape(itemsize, n).T.copy().tobytes()
order = "<" if "ENDIANNESS=LITTLE" in info else ">"
chunk = np.frombuffer(buf, f"{order}i2").reshape(chunks)
f"{len(raw):,} bytes on disk, {chunk.nbytes:,} bytes decoded, shape {chunk.shape}"

['COMPRESSION=DEFLATE', 'ENDIANNESS=LITTLE', 'FILTER=SHUFFLE']

'667,794 bytes on disk, 2,073,600 bytes decoded, shape (1, 1, 720, 1440)'

and it matches what gdalxarray reads through the VRT, value for value:

via_gdal = xr.open_dataset(vrt, engine="gdalxarray", mask_and_scale=False)["sst"].isel(time=1).values
assert np.array_equal(chunk[0], via_gdal)
bool(np.array_equal(chunk[0], via_gdal))

True

That equality is what makes a byte table a faithful replacement for the files’ own metadata: the files can stay where they are, and the table travels.

kerchunk-Parquet Zarr

GDAL’s Zarr driver reads kerchunk references stored as Parquet, so a reference store opens like any Zarr. This one is a pinned commit of the BRAN2023 ocean reanalysis temp variable: about 13 TB of netCDF on NCI’s THREDDS server, described by a few megabytes of Parquet on GitHub.

NoteLive: BRAN2023 ocean temp as kerchunk-Parquet Zarr (GitHub, pinned commit)

• Address: ZARR:"/vsicurl/https://raw.githubusercontent.com/mdsumner/virtualized/d0b9c23317d2de07c3851eb9dfc968be1761d5b8/zarr/ocean_temp.zarr"
• Provider: mdsumner/virtualized at d0b9c23; chunks on NCI THREDDS
• Pinned: a fixed commit or immutable object, will not change.
• Status: answered the store check at 2026-11-01 09:15 UTC.

if live("bran-virtualized-pinned"):
    bran = xr.open_dataset(store("bran-virtualized-pinned")["dsn"], engine="gdalxarray",
                           drop_variables=["Time_bnds"])
    print(f"{bran['temp'].nbytes / 1e12:.1f} TB described; dims {dict(bran.sizes)}")
    print(bran["temp"].encoding.get("preferred_chunks"))

12.9 TB described; dims {'Time': 5844, 'st_ocean': 51, 'yt_ocean': 1500, 'xt_ocean': 3600, 'nv': 2}
{'Time': 1, 'st_ocean': 1, 'yt_ocean': 300, 'xt_ocean': 300}

Reading values fetches the referenced byte ranges from THREDDS, which rate-limits parallel requests; the BRAN2023 gallery page shows that read and the settings it needs.

Rendered with GDAL 3.13.3, xarray 2026.4.1.dev32+g71cfb13ac, numpy 2.5.4, gdalxarray 0.5.1, rasterix 0.2.4.