Earth observation and ocean model archives are mostly not cloud-native, and they will not be rewritten. They do not need to be. Every chunk of every array already sits at a known byte range in a file that someone is already serving. Record those ranges once, in a plain table, beside a concatenation plan that GDAL already understands, and the archive can be published as whatever its readers want: a GDAL VRT for R and GIS users, a Parquet table for query engines and browsers, kerchunk for xarray, an Icechunk repository for versioned pipelines. No copies, no new format, and no single pathway that everyone must adopt.
This post sets out that position, how I arrived at it, and where the work stands as of October 2026. It sits beside three other pieces that each state one stance and argue it throughout: sidebyside (warp and cell: two stances on viewing gridded data), regular-in-what (a grid is regular in a frame) and trianglewins (one input contract, many backends, judged by code that belongs to none of them).
The stance
A virtual store is not a format. It is a table of byte ranges plus a plan, and every format is a projection of that table.
Every array dataset is a set of granules (netCDF, HDF5, GeoTIFF, GRIB) in which each chunk of each array occupies a byte range. Kerchunk, VirtualiZarr, Icechunk, a GDAL multidimensional VRT, DMR++ and a COG’s tile index all reduce to knowing those ranges and pointing a reader at them. The tools disagree about where that knowledge lives: in memory, in a bespoke manifest, in a transactional repository. I put it in one durable, queryable, language-neutral table and treat each format as a rendering of it.
Four corollaries follow, and the rest of the post works them out.
- Separate what VirtualiZarr fuses. The concatenation plan (a
gdal mdim mosaicVRT), the byte references (a flat table) and the file addressing (source plus remap rules) are three things, each held once. - Endpoints are staged and complementary, not competing. The mdim VRT, a flat Parquet table, kerchunk Parquet and Icechunk each serve a different reader and usage model. None is “the” answer and none is derived from another; all come from the table.
- The catalog points at published forms; it never holds bytes. bowerbird knows how to get data that is not FAIR, sds knows the string that opens it now, and the byte-reference table sits between them.
- Scan once, where it is cheap; reference where it is durable. Chunk offsets do not change when a file is copied byte for byte, so a scan of a local HPC mirror can produce references to a public URL. That is what turns a research mirror into a publishing pipeline.
How the idea arrived
It came in pieces, mostly by building something and finding the same table underneath it.
| When | Piece | What it showed |
|---|---|---|
| May 2025 | zarrquet | VirtualiZarr’s kerchunk Parquet is a table of (path, offset, size, raw), and the array’s chunk keys belong in that table too |
| Feb 2026 | idea-nci-thredds | VirtualiZarr to Icechunk working on Pawsey, and the gotchas along the way |
| Mar 2026 | cogcache | COG tiles fetched by range, decoded in Rust, catalogued in DuckDB; warp as target chunk -> source chunk refs -> read -> warp -> write, bit-identical to GDAL. The reference table is already here, as a cache key |
| Mar 2026 | parq2ice | kerchunk Parquet to Icechunk in Rust; hit a ~2 GB manifest limit on a 16.7M-chunk BRAN array |
| Apr 2026 | pawsey-aad | bowerbird sync into Pawsey object storage, COG conversion, VirtualiZarr Parquet and curated listings, all scheduled |
| Jun 2026 | virtualized | BRAN2023 built twice, by GDAL and by VirtualiZarr, and compared |
| Jul 2026 | Icechunk pilots | Icechunk 2.1 built from tables: OISST in three access patterns; 5,800 REMA COGs as ~19.5M refs from TIFF IFDs |
| Jul 2026 | raad.virtualize | one central store of byte references, with published forms as projections; 16,387 OISST files scanned in ~20 minutes into SQLite |
| Jul 2026 | blocklist | the clean split: VRT for the logic, a chunk walk for the physics, kerchunk Parquet out; an access location versus a public one |
| Aug 2026 | aad-filelist | bowerbird listings as Parquet, and a richer reference schema (arrays, inline bytes, affine) |
| Oct 2026 | mdim-refs/0.2 | the two SQLite stores merged, with column names from GDAL’s block info; three readers agree byte for byte on OISST |
| Oct 2026 | wxylib | the table read in Rust and written to Icechunk; then the table as a cache, with held bytes and lineage |
The arc closes on itself. cogcache began as a tile cache keyed by byte range, and the latest cache design ends as a byte-range cache described by the same rows. I take that as evidence for the stance rather than a coincidence.
Why a table
The case for rows over objects:
- Rows compute without shared state. Scan each source in parallel and bind the results. There is no live lazy-array graph to serialise between workers; a Slurm array or a mirai pool is just
map -> rows -> rbind. - Rows are queryable. With explicit chunk keys (
dim_0 .. dim_n) a table answers “which chunks touch this box and these days” by predicate pushdown before any byte is read. kerchunk’s positional index cannot. - Rows carry anything. A per-chunk codec, a checksum, a modification time, inline bytes, a held copy, a spatial extent, a lineage pointer from a derived chunk back to its sources. Adding a column does not change the format.
- Rows are language-neutral. SQLite and Parquet open in R, Python, Rust, GDAL’s vector drivers and DuckDB with no shared library.
- Rows make a cache honest. bowerbird mirrors whole files eagerly; a range cache mirrors only what was looked at. The same rows describe both, with a
heldstate (ref,inline,copy,link).
Five existing schemas, built independently, converge on the same core:
array, chunk key (explicit), location, offset, length, present
+ array metadata held once (dims, shape, chunks, dtype, fill, codec, attrs)
The mdim-refs/0.2 schema settles the remaining choices: explicit keys; chunk rows keyed within their source, so a rescan is an upsert and the global index is a view; locations recorded verbatim at scan time, with remap rules applied at publish; 64-bit byte columns; and one codec per array, so a change of codec ends one array and starts another (an era).
The wider argument is in Where the line falls: tensors, tables and four layers. Every dataset is values over an index space plus a map from index to place, either implicit (a recipe) or explicit (stored values). A tensor is a table whose key is implicit; a table is a tensor whose key is explicit. Byte references are the clearest case of the two needing each other: a table that describes a tensor.
The model
discover bowerbird, blueant, CMR, STAC, THREDDS -> sources (access, public, mtime)
plan gdal mdim mosaic -> mdim VRT
physical GDAL block info | HDF5 chunk walk | DMR++ | TIFF IFDs -> chunk rows per source
table mdim-refs: sources, arrays, chunk rows, placements, remap -> refs_<array> view
emit Icechunk | kerchunk Parquet | flat Parquet | OGR layer
publish object storage (the truth) -> public mirrors
register one sds row per published form
read GDAL, gdalraster, xarray, zaro, wxylib, DuckDB, ...
The table has three layers, each held once:
- Physical. For each source and array, the byte range of every chunk, keyed by its position inside the source. No time axis, no dataset, no remapping.
- Plan. Where each source’s slab sits in the global array, copied from the mosaic VRT so that the table stands alone outside GDAL.
- View. The global chunk index,
refs_<array>: exactly the rows that GDAL’s block info would emit if run on the mosaic as one array.
The view is a simple formula over the key, dim_i(global) = dim_i(source) - src_offset_i / c_i + dest_i / c_i, valid when every placement is aligned to chunk boundaries. If a mosaic has no single regular chunk grid it cannot be virtualised, and the table says so rather than writing a wrong store.
The full specification is mdim-refs/0.2.
The endpoints, staged and complementary
Each endpoint earns its place by serving a reader the others do not.
| Endpoint | Who reads it | Usage | Strength | Limit |
|---|---|---|---|---|
| mdim VRT | GDAL and everything on it (gdalraster, osgeo.gdal, R and Python bindings) | open one string, GDAL does the rest; also the build plan for the table | R access on day one; inspectable XML; no new format | slow to build for many sources; drops some CF attributes; reads go through the source files |
| flat Parquet | DuckDB, arrow, zaro, wxylib, any range fetcher | query the table, fetch only the ranges needed; the natural home of a cache | pushdown, no codec registry, browser-friendly, append-only | not a standard; every reader is ours |
| kerchunk Parquet | xarray and fsspec; GDAL’s Zarr driver via /vsikerchunk_parquet_ref/ |
interoperate with the Python stack, and with GDAL as Zarr | widest existing reader base | positional keys; full rewrite on update |
| Icechunk | icechunk-python, zarrs, icechunk-js | transactional history: bulk load once, append days as commits | versioning, checksums, browser readers | fast-moving format; manifest layout needs care for daily appends |
| SQLite / OGR layer | DBI, GDAL vector drivers | the canonical store, inspected with ogrinfo | one file per dataset, no server | not a reader format |
The VRT and the table are the heart. kerchunk comes first among the emitters because GDAL reads it, so R users get it. Icechunk is a sink where versioning matters. The flat Parquet export is a middle stage rather than a peer of the renders: harvest and collate into it, render from it, and keep the cache there. Its file layout is not fixed yet; a single self-describing file with Parquet key-value metadata is the likely end.
What Icechunk does, as table operations
A suspicion worth testing: much of Icechunk’s machinery is a view, or a simple operation, over a byte-reference table. The people building it frame it in similar terms. Tom Nicholas and Sebastian Galkin describe virtual chunks as playing “a similar role that an ‘external table’ does in more traditional databases”, and Nicholas’s ESIP talk treats the manifest as an abstraction layer in which physical and logical layout change independently. Earthmover’s recent posts set Icechunk beside Iceberg (“Iceberg is tabular, while Icechunk is tensor-based”) and put tables and tensors side by side in one platform.
Checked against the Icechunk 2.x specification, the suspicion holds for the addressing and layout features. It does not hold for the versioning features, which need real storage semantics: atomic, conditional writes on object storage. “View” below means a pure function of rows the table already holds; “layout” means a physical choice about storing those rows; “storage” means it needs atomic writes that no view can supply.
| Icechunk feature | Table equivalent | Kind |
|---|---|---|
| Virtual chunk containers; 2.0 named containers, so assets relocate by changing the container | a remap table (prefix, replacement) joined at read by longest-prefix match | view |
| Container credentials and read-time authorization | a per-prefix access policy, closer to row-level security | policy |
shift_array, reindex_array (regular grids only) |
dim_i + k over the chunk key: the placement formula above, with the same alignment restriction |
view |
move (metadata only) |
update one row of arrays |
one-row update |
| Chunk references (index, location, offset, length), inline, native | the refs_<array> columns |
same rows |
| References sorted by index, binary search | a clustered primary key | layout |
| Location dictionary | dictionary encoding of the path column | layout |
| Manifest splitting by extent | partitioning by floor(dim_0 / n), with min-max statistics for pruning |
layout |
| Checksums by etag or modification time | an etag or mtime column plus a reader-side check | column + reader policy |
| Diff between commits | EXCEPT of two snapshot row sets |
view, given snapshots |
| Snapshots, branches, tags | append-only rows with validity intervals and an “as of” view; publishing that atomically is the hard part | storage |
| Transactions and conflict detection | the conflict test is an intersection of changed chunk keys, a join; the commit itself is a conditional write | storage |
| Expiration and garbage collection | an anti-join of reachable against listed objects | storage operations |
Four readings follow, each skeptical in both directions.
Remap is the clearest case, and Icechunk has now built it. Named containers are a two-column prefix table resolved at read time. The mdim-refs table binds remap rules at publish instead, so relocating data means re-rendering. Writing one named container per remap rule would let an Icechunk repository carry the remap table itself, bound late.
Layout features are not free just because they are views. Splitting, sorting and dictionary coding are where Icechunk spends its engineering. A Parquet table gets the equivalents from Parquet and DuckDB, but only for readers that are ours. The 2 GB BRAN limit was a layout limit (one manifest), not a table limit, and manifest splitting addresses it.
Versioning is not a view. Snapshots-as-views are easy to write and wrong without atomic publication: a reader must see one consistent snapshot across all arrays, and two writers must not both win. Icechunk does this with conditional writes and no server. Among table formats, Lance does the same (a put-if-absent commit of each version manifest, with optimistic retry and no catalog); Delta Lake can, where the object store supports conditional writes; Iceberg and DuckLake need a catalog. If the byte-reference table needs history, it should borrow a table format with transactions rather than grow one.
Most of what we version is already a column. Preliminary and final versions of a day are two sources. A daily append is new rows keyed by source. “Which granules did this dataset use on a given date” is a question for the plan. That covers most of the need Icechunk’s versioning answers for virtual data, which is why Icechunk stays a sink rather than the store of record.
Where versioning pays: held bytes
That last point was tested on the GHRSST MUR 2024 store (935,496 references), comparing Lance, Parquet fragments with a simple manifest commit, and Icechunk with daily commits (design note).
- For references, history is mostly a column. Lance holds the mdim-refs tables unchanged, and snapshots, tags, diff, upsert and the remap view all work, with bytes matching GDAL. Nothing reading these stores needs “as of” yet. If it ever does, Lance is the better route than a thin homemade layer, which handled appends and a two-writer race and then rewrote half a year of rows on the first rescan.
- For bytes, versioning pays. This is cogcache’s model: references in a table until a chunk is needed, then the bytes held in the same row. Lance’s blob column stores each value either inline or as an external (uri, position, size) range, and reads both the same way. A half-held, half-external window over Hobart matches GDAL. Taking one held chunk by row costs 1 ms, against 29 ms from Parquet with 64-row row groups. A reprocessed day changes only the new version, and an older tag still reads the bytes that were held. A plain reference cannot do that.
- Icechunk with daily commits costs 22 MB of history for a year when manifests are split by day, and 1.84 GB with a 365-day split. Split by the unit of append.
Knowledge and catalog
The long-running AAD data stack already holds most of what a catalog needs; the virtual store gives each piece a clearer role.
| Layer | Pieces | Keeps |
|---|---|---|
| acquisition | bowerbird and its handlers, blueant (109 sources) | how to reach each provider; the dataset id |
| listing | aad-filelist, raadfiles | what files exist, and where |
| sets | raadfiles finders, date and deduplication rules | which files form a dataset; preliminary versus final |
| freshness | CMR production times, file fingerprints | what changed, and what to rescan |
| reading, old | raadtools read*() |
per-dataset code that the virtual store retires |
| reading, new | GDAL, zaro, lazy array readers | one lazy array per dataset |
| registry | sds | the string that opens it now |
- bowerbird stays the acquisition layer for data that is not FAIR. Its configurations are knowledge no aggregator holds. The bowerbird mirror is the
accesslocation; the provider URL is thepublicone. - raadfiles file sets are mosaic inputs, and raadtools’ edge cases (zero to three variables per file, a percentage becoming a fraction) are codec boundaries: eras, recorded once as groups in one store.
- sds stays strings-only and does not own the table. It gains kinds for each published form (mdim VRT, kerchunk Parquet, mdim-refs, Icechunk, STAC, CMR) and a link from each row back to its bowerbird id. Recipe functions become builder scripts that emit stores, and each recipe becomes one row.
- A row earns its place if sds adds what a catalog cannot: a recipe, a virtual store, a verified choice among versions, or a status nobody else tracks. Where STAC already serves a dataset, the row points at the collection.
- The registry is published as data, as STAC or stac-geoparquet built in CI, so that Python, GDAL and DuckDB users can contribute.
Where good tools exist they are adopted rather than rebuilt: STAC, CMR, earthdatalogin for authentication, and DMR++ as a ready-made physical layer for NASA collections.
From HPC to the public
The path for the work is to develop where the data and compute already are, publish to storage under our control, and let a cheaper and more public layer serve it.
develop Nectar (raad mirror), Pawsey Setonix + Acacia, NCI (BRAN in situ)
bowerbird sync, scans, VRTs, tables; Slurm + Singularity images
publish Pawsey object storage: the single source of truth
present source.coop: byte-exact sync of designated prefixes,
URLs rewritten only by the publisher's remap rules
code GitHub repositories and releases; Actions for small appends and checks
operate cloud functions and instances for scheduled builds and appends (to come)
The rules are simple. Pawsey is the truth; nothing appears on source.coop that is not on Pawsey first; references are never edited by hand downstream; remap rules stay with the publisher. The open leg is operations: moving from virtualisation work on research infrastructure to scheduled production that feeds stores, or indexes of stores, on source.coop. ghrsst-cogger is a live example of that shape. A block-info scan is only metadata reads, so per-granule scans and daily appends suit small serverless jobs, with larger instances for backfills.
Tooling: GDAL first, and a Rust experiment
GDAL and R
In order of leverage:
- GDAL block info. A GDAL change that makes chunk references for every multidimensional driver through one code path, with a fast path for HDF5 chunk iteration. A first pull request was paused after review and is being reworked under that name. It is written by hand.
- An R binding of
GetRawBlockInfoin gdalraster. The function is already in stock GDAL and the binding is small and independent of (1). Until it exists the scanner is Python. - GDAL as a first-class reader of the stores: the numcodecs shuffle and zlib filters on the Zarr v3 path, TIFF predictors, the LZW and floating-point predictor codecs that REMA uses, and reading Parquet key-value metadata in place of sidecar JSON.
- R readers: zaro (which already reads the mdim-refs store), support for the fsspec kerchunk layout, a batched range fetcher, and GDAL as the common face.
- A released table producer in R. The only R implementation is not yet packaged; raad.virtualize is the natural home for the store and blocklist for the kerchunk writer.
wxylib: an experiment in Rust
wxylib is a Rust project built quickly from the best available crates, with GDAL as its test oracle and design documents gating each stage. The point is to push it hard now and see whether it deserves a long-term commitment.
What it has shown for this scheme:
- the table is right: its own design principle, “a scanner that emits chunk references rather than a reader”, is this stance reached independently;
- it reads the mdim-refs Parquet export, maps codecs to Zarr v3, honours remap, writes Icechunk and reads it back, matching GDAL byte for byte on OISST; BRAN (97,920 references, many chunks per file) is structurally right, with bytes not yet verified;
- fixtures plus an oracle give the table the conformance test it never had;
- the table can be a cache with lineage, which is new and useful.
What it has not shown yet: a table producer of its own, a read path that goes directly through the flat table rather than Icechunk, a kerchunk writer, and an R binding that R users can install without a Rust toolchain. Much of it was written quickly with generated code, so correctness is guarded by the oracle more than by understanding.
It earns a long-term commitment if four things hold: the flat-table reader works with the cache inside it, with BRAN verified byte for byte against GDAL; it consumes a table producer it does not own; an R user can read a table through it without a Rust toolchain (or R is decided not to be its goal); and it introduces no new file format. If those hold, wxylib is the engine of the scheme. If not, it is a strong consumer of it, which is still worth a lot.
Reader coverage: the GHRSST COG run
GHRSST MUR analysed_sst COGs from source.coop (2024, 935,496 references) went through raad.virtualize and were published to a test bucket, then verified from the bucket: the VRT and Icechunk forms match GDAL on every window. The source encoding (Int16, ZSTD, TIFF horizontal predictor) splits the readers, and shows why several views of one table matter:
| Form | Python (imagecodecs) | GDAL and R |
|---|---|---|
| mdim VRT | yes | yes, but a full-width row took 166 s |
| Icechunk | yes, the same row in 6 s | not yet (predictor codec) |
| kerchunk Parquet | yes | not yet (predictor codec) |
| flat refs store | own decoder | through wxylib, once it decodes the predictor |
Widening that coverage for R users means a TIFF predictor filter in GDAL’s Zarr reader, predictor decoding in wxylib’s table reader, faster COG tile fetching through the VRT, and recording attribute types in the store so that packed integers are not silently decoded to doubles.
Where it stands, October 2026
Working now
- The mdim-refs/0.2 table, with OISST daily (16,387 sources in a 14 MB SQLite file) read byte for byte the same by three readers.
- All four endpoints built from that one table for OISST (mdim VRT, kerchunk Parquet, flat Parquet and Icechunk), and largely verified.
- A daily-updated OISST mdim VRT on Pawsey object storage.
- GHRSST MUR 2024 through the pipeline to object storage, with the VRT and Icechunk forms verified against GDAL.
- Lance tested as a versioned store for held bytes, with reads matching GDAL.
In progress
- GDAL block info, being reworked after a first review.
- BRAN2023 (4D, many variables, ~16.8M chunks per variable, many chunks per file) built from NCI in situ, structurally right and awaiting byte verification.
- The flat Parquet layout, and a direct reader for it with the cache inside.
- Known rough edges: a scale-factor rounding difference in the R-built BRAN index, a fill-value disagreement between VirtualiZarr and the source attributes, and GDAL 3.13’s block info erroring on unwritten chunks.
Next
- A gdalraster binding of
GetRawBlockInfo, and a released R table producer. - MUR SST three ways (DMR++, HDF5 scan, COG mirror) into one table.
- sds rows for each published store, with the registry published as STAC.
- Scheduled production on source.coop.
The end state
- A dataset’s knowledge lives in a bowerbird source, its bytes in the provider’s files, its plan in a VRT, its references in one table, and its address in one sds row.
- A daily append is a scan of the new files and an upsert of their rows, run where the data is, published to object storage and mirrored to source.coop by rule.
- R, Python, GDAL, DuckDB, Rust and the browser read the same table, and a conformance fixture proves they read the same values.
- A range cache, described by the same rows, makes repeated access to remote archives cheap and records where every derived chunk came from.
The bytes stay where they are. Everything else is a table, and a choice of how to render it.