Virtual stores: the archive as one array

The COGs are the storage. A virtual store is an index over them that lets the whole archive be opened as a single array (time, y, x) for each variable, without copying or converting any data. This page covers how that works, what you can do today, including the GDAL multidimensional VRT, and the indexes that are coming.

A COG is already chunked
What you can do today
GDAL multidimensional VRT
Coming
Things the indexes have to get right
Why not convert to Zarr?

A COG is already chunked

A tiled GeoTIFF is a chunked array with its chunk manifest in the header. Each day’s analysed_sst COG is a 17999 by 36000 array cut into 512 by 512 tiles; its header lists the byte offset and length of every tile, then the same again for each overview level. A reader that has the header can fetch any tile with one HTTP range request.

You can see this for one file on the playproj tiles and overviews page: it reads only the header, draws the tiles and overview levels, and gives the byte range of any tile you click.

So the archive is already a chunked cube: about 8,900 days, each 36 by 71 tiles per variable, plus overviews. What’s missing is a single description of the cube: the time axis, which file holds which day, and where each tile is, in a form that array tools can open directly.

What you can do today

Build the cube from the catalogue at read time. Each tool reads the STAC index, opens each day’s COG header, and assembles a lazy array:

Tool	How
Python, odc-stac	`odc.stac.load(items, chunks={})` gives a dask-backed xarray cube (example)
Python, stackstac	`stackstac.stack(items)`, the same idea
R, terra	`rast(urls)` with `time()` set, one layer per day (example)
GDAL	`gdal raster stack` over a list of URLs, one band per day (example)
GDAL multidimensional	The published VRT, every day as one array per variable (below)

The first four build the cube in each session: they read the index and then the header of every day they use, a few hundred requests for a year, before any pixels.

GDAL multidimensional VRT

ghrsst-mur-v2.vrt describes every COG as one GDAL multidimensional dataset. It’s rebuilt from the STAC items each day, with the STAC Parquet file.

Arrays: analysed_sst, analysis_error, mask, sea_ice_fraction and sst_anomaly, each (time, y, x), plus the time, y and x coordinates.
Time: one step per day from 2002-06-01, in days since 1970-01-01.
Grid: the shared 0.01 degree grid, with y and x at pixel centres.
Values: scale, offset, nodata and unit as in the COGs, so unscaled analysed_sst is in degrees Celsius. Days without a variable (sst_anomaly before 2014) read as nodata.

The VRT holds no data. Each time step points at that day’s COG, and GDAL reads the tiles as it would for any GeoTIFF. Opening it reads one file (about 17 MB, growing with each day) instead of the STAC index, but a read still opens the header of each day it touches, so a long time series at one point costs a request or two per day. Set GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR so GDAL doesn’t also list each day’s folder: it roughly halves the time.

export GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR
gdal mdim info /vsicurl/https://data.source.coop/ausantarctic/ghrsst-mur-v2/ghrsst-mur-v2.vrt

Python, with GDAL’s bindings:

from osgeo import gdal

gdal.UseExceptions()
gdal.SetConfigOption("GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR")
ds = gdal.OpenEx("/vsicurl/https://data.source.coop/ausantarctic/ghrsst-mur-v2/ghrsst-mur-v2.vrt", gdal.OF_MULTIDIM_RASTER)
root = ds.GetRootGroup()
sst = root.OpenMDArray("analysed_sst")

# The last 7 days at 44S 150E, finding the row and column from the coordinates
y = root.OpenMDArray("y").ReadAsArray()
x = root.OpenMDArray("x").ReadAsArray()
row, col = int(abs(y + 44).argmin()), int(abs(x - 150).argmin())
n = sst.GetDimensions()[0].GetSize()
series = sst.GetUnscaled().ReadAsArray(
    array_start_idx=[n - 7, row, col], count=[7, 1, 1]
).ravel()  # degrees Celsius

Other readers of GDAL’s multidimensional API can open it too:

GDAL: gdal mdim convert to subset or copy out a region and time range
R: terra::rast(x, md = TRUE), or gdalraster’s multidimensional functions
Python: xarray through a GDAL backend such as gdalxarray

Coming

Three more indexes over the same COGs, to be published alongside them on Source Cooperative. Each records where every tile is, so unlike the VRT, opening one doesn’t read any COG headers.

A flat tabular index

One row per tile: date, variable, overview level, tile column and row, the COG’s URL, the tile’s byte offset and length, and the COG’s checksum. This is the next level down from the STAC GeoParquet, which has one row per day. As a Parquet table it can be queried with DuckDB, arrow or Polars to plan exactly which byte ranges a region and date range need, and the other indexes can be generated from it.

kerchunk Parquet references

A kerchunk reference store in Parquet: a Zarr (version 2) view of the archive in which each chunk is a reference to a tile’s byte range in a COG. It opens with fsspec’s reference file system and xarray’s Zarr engine, in any environment that has them, with no GDAL.

Icechunk

An Icechunk store whose chunks are virtual references into the COGs: a Zarr (version 3) view, with Icechunk’s transactions and version history. When a preliminary day is replaced by its final, the update is a new commit, so earlier analyses can still name the exact version they used.

Things the indexes have to get right

Codecs. The tiles are ZSTD-compressed with horizontal differencing (TIFF predictor 2). GDAL handles that natively. Zarr-based stores have to declare matching codecs for every chunk.
Edges. The grid is 17999 by 36000, which isn’t a whole number of 512-pixel tiles. TIFF stores the edge tiles at full size, padded, which is what a regular Zarr chunk grid expects.
Overviews. Each overview level is its own array, with its own tiles: a set of multiscales rather than one array.
Days without sst_anomaly. Older source files don’t have it, so that array has missing time steps.
Change. Preliminary days are rewritten at the same URLs when their final arrives, so references to them go stale. The VRT points at files, not byte ranges, so it isn’t affected. The other indexes will record each COG’s file:checksum and be rebuilt as days are finalised. final days don’t change.

Why not convert to Zarr?

A Zarr copy would duplicate several terabytes that are already chunked, compressed and public, and it would need its own process to keep it in step as days are finalised. The virtual stores give array tools the same view with only an index to maintain, and GDAL-based tools keep reading the COGs directly.