Can't isel on VirtualiZarr of multi-dimensional NetCDF files

[alexgleith]

I have what I would hope is a simple use case for VirtualiZarr, which is to create a lazy-loaded xarray of existing NetCDF files.

Code I'm using currently is this:

from virtualizarr import open_virtual_mfdataset
from obstore.store import HTTPStore
from virtualizarr.parsers import HDFParser

files = [
    "https://thredds.nci.org.au/thredds/fileServer/gb6/BRAN/BRAN2020/month/ocean_temp_mth_1993_01.nc",
    "https://thredds.nci.org.au/thredds/fileServer/gb6/BRAN/BRAN2020/month/ocean_temp_mth_1993_02.nc",
]

parser = HDFParser()
store = HTTPStore(url="https://thredds.nci.org.au")

ds = open_virtual_mfdataset(files, object_store=store, parser=parser)
ds

This works ok, although datetimes are in float64.

When it comes to selecting a single time and depth though, it fails.

one = ds.temp.isel(Time=0, st_ocean=0)
one

With this error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 1
----> 1 one = ds.temp.isel(Time=0, st_ocean=0)
      2 one

File /opt/homebrew/lib/python3.11/site-packages/xarray/core/dataarray.py:1555, in DataArray.isel(self, indexers, drop, missing_dims, **indexers_kwargs)
   1550     return self._from_temp_dataset(ds)
   1552 # Much faster algorithm for when all indexers are ints, slices, one-dimensional
   1553 # lists, or zero or one-dimensional np.ndarray's
-> 1555 variable = self._variable.isel(indexers, missing_dims=missing_dims)
   1556 indexes, index_variables = isel_indexes(self.xindexes, indexers)
   1558 coords = {}

File /opt/homebrew/lib/python3.11/site-packages/xarray/core/variable.py:1023, in Variable.isel(self, indexers, missing_dims, **indexers_kwargs)
   1020 indexers = drop_dims_from_indexers(indexers, self.dims, missing_dims)
   1022 key = tuple(indexers.get(dim, slice(None)) for dim in self.dims)
-> 1023 return self[key]

File /opt/homebrew/lib/python3.11/site-packages/xarray/core/variable.py:769, in Variable.__getitem__(self, key)
    766 dims, indexer, new_order = self._broadcast_indexes(key)
    767 indexable = as_indexable(self._data)
--> 769 data = indexing.apply_indexer(indexable, indexer)
    771 if new_order:
    772     data = duck_array_ops.moveaxis(data, range(len(new_order)), new_order)

File /opt/homebrew/lib/python3.11/site-packages/xarray/core/indexing.py:1029, in apply_indexer(indexable, indexer)
   1027     return indexable.oindex[indexer]
   1028 else:
-> 1029     return indexable[indexer]

File /opt/homebrew/lib/python3.11/site-packages/xarray/core/indexing.py:1555, in NumpyIndexingAdapter.__getitem__(self, indexer)
   1551 # We want 0d slices rather than scalars. This is achieved by
   1552 # appending an ellipsis (see
   1553 # https://numpy.org/doc/stable/reference/arrays.indexing.html#detailed-notes).
   1554 key = indexer.tuple + (Ellipsis,)
-> 1555 return array[key]

File /opt/homebrew/lib/python3.11/site-packages/virtualizarr/manifests/array.py:261, in ManifestArray.__getitem__(self, key)
    256     raise TypeError(
    257         f"indexer must be of type int, slice, ellipsis, None, or np.ndarray; or a tuple of such types. Got {key}"
    258     )
    260 # check value is valid
--> 261 indexer = _possibly_expand_trailing_ellipsis(indexer, self.ndim)
    262 if len(indexer) != self.ndim:
    263     raise ValueError(
    264         f"Invalid indexer for array. Indexer length must be less than or equal to the number of dimensions in the array, "
    265         f"but indexer={indexer} has length {len(indexer)} and array has {self.ndim} dimensions."
    266         f"\nIf concatenating using xarray, ensure all non-coordinate data variables to be concatenated include the concatenation dimension, "
    267         f"or consider passing `data_vars='minimal'` and `coords='minimal'` to the xarray combining function."
    268     )

File /opt/homebrew/lib/python3.11/site-packages/virtualizarr/manifests/array.py:366, in _possibly_expand_trailing_ellipsis(indexer, ndim)
    364 if final_dim_indexer == ...:
    365     if len(indexer) > ndim:
--> 366         raise ValueError(
    367             f"Invalid indexer for array. Indexer length must be less than or equal to the number of dimensions in the array, "
    368             f"but indexer={indexer} has length {len(indexer)} and array has {ndim} dimensions."
    369             f"\nIf concatenating using xarray, ensure all non-coordinate data variables to be concatenated include the concatenation dimension, "
    370             f"or consider passing `data_vars='minimal'` and `coords='minimal'` to the xarray combining function."
    371         )
    373     extra_slices_needed = ndim - (len(indexer) - 1)
    374     *indexer_as_list, ellipsis = indexer

ValueError: Invalid indexer for array. Indexer length must be less than or equal to the number of dimensions in the array, but indexer=(0, 0, slice(None, None, None), slice(None, None, None), Ellipsis) has length 5 and array has 4 dimensions.
If concatenating using xarray, ensure all non-coordinate data variables to be concatenated include the concatenation dimension, or consider passing `data_vars='minimal'` and `coords='minimal'` to the xarray combining function.

Doing a regular load like this works fine (with the caveat that it's reading from opendap URLs):

ds = xr.open_mfdataset(
    files_od,
    decode_timedelta=True,
    engine="netcdf4"
) 
one = ds.temp.isel(Time=0, st_ocean=0)
one.plot.imshow()

VirtualiZarr version is virtualizarr-1.3.3.dev81+ga5d04d7 (latest dev as of 1 July 2025)

[TomNicholas]

Hi Alex, thanks for trying out VirtualiZarr. However what you're trying to do here (index into and load a "virtual" dataset directly) betrays an important misunderstanding of what VirtualiZarr does and how it works.

Please read the "Important" box that starts with "Virtual datasets are not normal Xarray datasets" in this section of the docs.

[TomNicholas]

(The error you actually see is because of #51)

[alexgleith]

The only use case for a virtual dataset is combining references to files before writing out those references to disk.

So, wait, VitualiZarr is not aimed at using the data it's virtualising, but preparing a representation before it's being written to disk? Or does this mean that it's meant to be used to create a Kerchunk reference, which can then be actually used?

I guess I am missing the whole "why" of the tool! I thought it was to create a virtual Zarr, which could be used in xarray.

[TomNicholas]

So, wait, VitualiZarr is not aimed at using the data it's virtualising, but preparing a representation before it's being written to disk? Or does this mean that it's meant to be used to create a Kerchunk reference, which can then be actually used?

This is correct.

I thought it was to create a virtual Zarr, which could be used in xarray.

As the zarr spec does not have virtual references, currently a "virtual zarr" is either a set of kerchunk references, an icechunk store, or an (in-memory) ManifestStore. So there's no contradiction here, just a confusingly-named range of format choices.

I guess I am missing the whole "why" of the tool!

Creating a virtual Zarr is useful because you're providing cloud-optimized access to potentially non-cloud-optimized formats. Outside of object storage this is less useful.

VirtualiZarr provides a general framework for doing so for any file format that can possibly be accessed as "virtual" zarr. By using a parser, VirtualiZarr can create a virtual zarr of one file. Recently we added the ability to load this as a "normal" xarray dataset without serializing to disk first (the ManifestStore). Note that the utility of that is diminished for a file format which is already cloud-optimized on a per-file basis (e.g. COG).

VirtualiZarr also provides functionality for combining the virtual references for many files into a single big datacube efficiently, using xarray concat/merge etc. as the combining API. You do this using open_virtual_dataset/open_virtual_mfdataset + xr.concat/merge. Once you serialize the combined virtual references to kerchunk/icechunk now you have one big virtual zarr store pointing at potentially millions of individual archival files, i.e. you have cloud-optimized an entire collection of related files. Note at that point if you wrote to kerchunk the result should be exactly the same as calling kerchunk's MultiZarrToZarr, it's just less painful, more generalizable, and you could have alternatively used the same steps to write to Icechunk.

[mdsumner]

Ok so you just have to write it first, then load it, then it works. That could be stated more plainly and perhaps framed as a ux inconvenience. I take it we're only really meant to use methods on .virtualize before serializing and not much else.

Where in the onion are the byte references when loaded from kerchunk? Prior to serializing they are in ds.var.data.manifest.

[maxrjones]

thanks for raising the issue!

what I would hope is a simple use case for VirtualiZarr, which is to create a lazy-loaded xarray of existing NetCDF files.

this statement is absolutely true - it is a simple use case. For that reason, I think we should seriously consider this addition since it would add little maintenance burden and support a real use-case. In addition to the original example in this issue, it would be the simplest path for an xarray backend for async-tiff via virtual-tiff (though many performance benefits will be lost without an async virtualizarr API).

Please see #654 for a rough proof of concept.

[alexgleith]

Hey @maxrjones, here's my naive take on a COG-backed xarray, which didn't work, unsurprisingly!

#648

I'm honestly not sure how odc.stac.load actually frames up the .tif access and passes it on to Rasterio, but it works, and it works very fast. So as a happy, naive and enthusiastic user, I like it a lot. I'm interested to know what the virtual-tiff does that is different to this method?

[TomNicholas]

Ok so you just have to write it first, then load it, then it works.

Well now that we have ManifestStore, you actually don't need to write it first, if you just don't bother creating a virtual dataset. This already works without writing anything:

manifest_store = HDFParser()('file.nc')
ds = xr.open_zarr(manifest_store)

This is now a normal xarray dataset, and we never had to save it to kerchunk/icechunk first.

That could be stated more plainly and perhaps framed as a ux inconvenience.

That's fair, we definitely need to document this better, as the choice between calling the parser directly and passing the parser to open_virtual_dataset is confusing.

you just have to write it first, then load it

We have discussed at length the idea of making loading into memory very simple (e.g. we could make vds.load() work), but the downside of making it easy is that it creates a massive footgun: you could have a virtual dataset that refers to a PetaByte of data, and unknowingly attempt to fetch all that data. To some extent it's better that we deviate from the behaviour of normal xarray because you're not supposed to use this exactly like how you normally use xarray.