Slicing error when concatenating virtual datasets of a netCDF4 root group

[danielfromearth]

Description

I'm working to generate a combined virtual dataset from grouped netCDF4 files in the NASA Earthdata Archive. To do so, I am opening each group individually for each dataset (and this works). However, when (for each group) I try to concatenate the datasets along an existing dimension, it fails on the 'root' group, even though it works on other subgroups.

The error being raised is NotImplementedError: Doesn't support slicing with (None, slice(None, None, None)).

Related issues

I see there are other issues related to groups, but most seem to be about developing more native handling of groups in line with Datatree, so not sure if they are helpful: #84, #336, #402.

To Reproduce - Example 1

import earthaccess
import xarray as xr
from virtualizarr import open_virtual_dataset  # seems to requires virtualizarr>=1.3.0 to open the root group for a single netCDF.

# 1. Authenticate for NASA Earthdata (requires a free Earthdata Login account).
earthaccess.login()

# 2. Get AWS credentials for the data archive.
fs = earthaccess.get_s3_filesystem(daac="ASDC")
reader_opts = {"storage_options": fs.storage_options}

# 3. For each of the following two netCDF links: the 'root' and 'product' groups are each opened as a virtual dataset.
data_s3links = [
    's3://asdc-prod-protected/TEMPO/TEMPO_NO2_L3_V03/2025.01.11/TEMPO_NO2_L3_V03_20250111T125225Z_S002.nc',
    's3://asdc-prod-protected/TEMPO/TEMPO_NO2_L3_V03/2025.01.11/TEMPO_NO2_L3_V03_20250111T133233Z_S003.nc'
]

vdataset_groups = {"": [], "product": []}
for url in data_s3links:
    vdataset_groups[""].append(
        open_virtual_dataset(url, group="", indexes={}, reader_options=reader_opts)
    )
    vdataset_groups["product"].append(
        open_virtual_dataset(url, group="product", indexes={}, reader_options=reader_opts)
    )

# 4. Concatenate together the two virtual datasets along the 'time' dimension.
xr_concat_kwargs = {
    "coords": "minimal",
    "compat": "override",
    "combine_attrs": "override",
}

# *** THIS DOES NOT RAISE THE ERROR ***
vds = xr.concat(vdataset_groups["product"], dim="time", **xr_concat_kwargs)

# *** THIS RAISES THE ERROR: `NotImplementedError: Doesn't support slicing with (None, slice(None, None, None))` ***
vds = xr.concat(vdataset_groups[""], dim="time", **xr_concat_kwargs)

Error traceback

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[3], line 38
     35 vds = xr.concat(vdataset_groups["product"], dim="time", **xr_concat_kwargs)
     37 # THIS RAISES THE ERROR: `NotImplementedError: Doesn't support slicing with (None, slice(None, None, None))`
---> 38 vds = xr.concat(vdataset_groups[""], dim="time", **xr_concat_kwargs)

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/structure/concat.py:277](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/structure/concat.py#line=276), in concat(objs, dim, data_vars, coords, compat, positions, fill_value, join, combine_attrs, create_index_for_new_dim)
    264     return _dataarray_concat(
    265         objs,
    266         dim=dim,
   (...)    274         create_index_for_new_dim=create_index_for_new_dim,
    275     )
    276 elif isinstance(first_obj, Dataset):
--> 277     return _dataset_concat(
    278         objs,
    279         dim=dim,
    280         data_vars=data_vars,
    281         coords=coords,
    282         compat=compat,
    283         positions=positions,
    284         fill_value=fill_value,
    285         join=join,
    286         combine_attrs=combine_attrs,
    287         create_index_for_new_dim=create_index_for_new_dim,
    288     )
    289 else:
    290     raise TypeError(
    291         "can only concatenate xarray Dataset and DataArray "
    292         f"objects, got {type(first_obj)}"
    293     )

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/structure/concat.py:676](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/structure/concat.py#line=675), in _dataset_concat(datasets, dim, data_vars, coords, compat, positions, fill_value, join, combine_attrs, create_index_for_new_dim)
    674         result_vars[k] = v
    675 else:
--> 676     combined_var = concat_vars(
    677         vars, dim_name, positions, combine_attrs=combine_attrs
    678     )
    679     # reindex if variable is not present in all datasets
    680     if len(variable_index) < concat_index_size:

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/core/variable.py:3018](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/core/variable.py#line=3017), in concat(variables, dim, positions, shortcut, combine_attrs)
   2970 def concat(
   2971     variables,
   2972     dim="concat_dim",
   (...)   2975     combine_attrs="override",
   2976 ):
   2977     """Concatenate variables along a new or existing dimension.
   2978 
   2979     Parameters
   (...)   3016         along the given dimension.
   3017     """
-> 3018     variables = list(variables)
   3019     if all(isinstance(v, IndexVariable) for v in variables):
   3020         return IndexVariable.concat(variables, dim, positions, shortcut, combine_attrs)

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/structure/concat.py:598](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/structure/concat.py#line=597), in _dataset_concat.<locals>.ensure_common_dims(vars, concat_dim_lengths)
    596 if var.dims != common_dims:
    597     common_shape = tuple(dims_sizes.get(d, dim_len) for d in common_dims)
--> 598     var = var.set_dims(common_dims, common_shape)
    599 yield var

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/util/deprecation_helpers.py:143](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/util/deprecation_helpers.py#line=142), in deprecate_dims.<locals>.wrapper(*args, **kwargs)
    135     emit_user_level_warning(
    136         f"The `{old_name}` argument has been renamed to `dim`, and will be removed "
    137         "in the future. This renaming is taking place throughout xarray over the "
   (...)    140         PendingDeprecationWarning,
    141     )
    142     kwargs["dim"] = kwargs.pop(old_name)
--> 143 return func(*args, **kwargs)

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/core/variable.py:1395](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/xarray/core/variable.py#line=1394), in Variable.set_dims(self, dim, shape)
   1388 elif shape is None or all(
   1389     s == 1 for s, e in zip(shape, dim, strict=True) if e not in self_dims
   1390 ):
   1391     # "Trivial" broadcasting, i.e. simply inserting a new dimension
   1392     # This is typically easier for duck arrays to implement
   1393     # than the full "broadcast_to" semantics
   1394     indexer = (None,) * (len(expanded_dims) - self.ndim) + (...,)
-> 1395     expanded_data = self.data[indexer]
   1396 else:  # elif shape is not None:
   1397     dims_map = dict(zip(dim, shape, strict=True))

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/virtualizarr/manifests/array.py:226](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/virtualizarr/manifests/array.py#line=225), in ManifestArray.__getitem__(self, key)
    224     return self
    225 else:
--> 226     raise NotImplementedError(f"Doesn't support slicing with {indexer}")

NotImplementedError: Doesn't support slicing with (None, slice(None, None, None))

[TomNicholas]

Hi @danielfromearth , thanks for raising this, and thank you for the reproducer.

requires a free Earthdata Login account

However, I do not have an Earthdata Login, and would rather not make one. Could you please help make a minimal reproducer either by:

1. Simplifying your example to the point that it uses fake data, not actual files that need to be downloaded from NASA,
2. Showing us a printout of the contents of the problematic datasets so that we can try to assemble such a minimal reproducer.

I suspect that the minimal reproducer for a group-related issue like this might only need one group with one small variable in it...

[danielfromearth]

Some good news

@TomNicholas, okay, in the process of creating a smaller example, I think I at least discovered a bit more of what is leading to the error. Here's what I'm seeing and what I've learned thus far...

What does not work

The below Example 2, which creates minimal netCDF files, generates the same NotImplementedError error in my environment (virtualizarr=1.3.2, xarray=2025.6.1, netCDF4=1.7.2). (Note that the error message is shorter (NotImplementedError: Doesn't support slicing with (None,)), though I suspect that's simply because this is an example with fewer variables (or dimensions), as opposed to the more extensive data from my original post's Example 1.)

What does work

In order to create this smaller example, I found the error was not being generated — until I added a data variable (in this case, "weight") into the root group. Note that the "weight" variable does not include the concatenation-dimension. So, that points to a possible solution, which is to add data_vars="minimal" to the xarray concatenation arguments. And that works! So, by adding that keyword argument to the code at the end of Example 2, the code no longer generates the error.

Example 2

import netCDF4  # not needed to be imported, but needs to be available in the environment for `DataTree.to_netcdf()`
import xarray as xr
from virtualizarr import open_virtual_dataset  # seems to requires virtualizarr>=1.3.0 to open the root group for a single netCDF.

# 1. Create a first grouped Datatree from two example datasets.
ds_root = xr.Dataset.from_dict({
    "time": {"dims": ("time"), "data": [0]},
    "longitude": {"dims": ("longitude"), "data": [0, 90]},
    "weight": {"dims": ("longitude"), "data": [50, 51]}
})
ds_product = xr.Dataset.from_dict({"a": {"dims": ("time", "longitude"), "data": [[10, 20]]}})
datatree1 = xr.DataTree.from_dict({"/": ds_root, "/product": ds_product})

# 2. Create a second grouped Datatree from two example datasets.
ds_root = xr.Dataset.from_dict({
    "time": {"dims": ("time"), "data": [1]},
    "longitude": {"dims": ("longitude"), "data": [0, 90]},
    "weight": {"dims": ("longitude"), "data": [50, 51]}
})
ds_product = xr.Dataset.from_dict({"a": {"dims": ("time", "longitude"), "data": [[11, 22]]}})
datatree2 = xr.DataTree.from_dict({"/": ds_root, "/product": ds_product})

# 3. Write Datatrees to netCDF.
data_paths = ["dataset1.nc", "dataset2.nc"]
datatree1.to_netcdf(data_paths[0])
datatree2.to_netcdf(data_paths[1])

# 4. For each group, open virtual datasets.
vdataset_groups = {"": [], "product": []}
for path in data_paths:
    vdataset_groups[""].append(open_virtual_dataset(path, group="", indexes={}))
    vdataset_groups["product"].append(open_virtual_dataset(path, group="product", indexes={}))

# 5. Concatenate together the two virtual datasets along the 'time' dimension.
xr_concat_kwargs = {
    "coords": "minimal",
    "compat": "override",
}

# *** THIS DOES NOT RAISE THE ERROR ***
vds_product = xr.concat(vdataset_groups["product"], dim="time", **xr_concat_kwargs)

# *** THIS RAISES THE ERROR: `NotImplementedError: Doesn't support slicing with (None, slice(None, None, None))` ***
vds_root = xr.concat(vdataset_groups[""], dim="time", **xr_concat_kwargs)

What does this mean?

I think this all makes sense to me, i.e., that there is no clear way to handle a variable for which the concatenation-dimension does not appear but is present across all data sets. And, for VirtualiZarr, since the arrays are represented as ManifestArrays, the variables also can't be compared without loading these variables (right?) — therefore, we must assume that the data variables that aren't being concatenated are either equal or unimportant beforehand... we can't lean on xarray to handle this with either data_vars="different" or with compat options other than "override"

Moreover, perhaps the error message or handling can be adjusted to more clearly indicate the possibility of these data variable handling options being a place to check.

[TomNicholas]

@danielfromearth I think you've just written up a good explanation of why the docs tell you to use coords='minimal', compat='override' (though they should presumably also tell you to use data_vars='minimal' too). 😄

Moreover, perhaps the error message or handling can be adjusted to more clearly indicate the possibility of these data variable handling options being a place to check.

I would welcome a PR to improve the error message that VirtualiZarr raises here.

[TomNicholas]

Also note our longer-term plan to make this easier for users is to change the defaults in xarray to be less of a bear trap: pydata/xarray#10062

[danielfromearth]

I would welcome a PR to improve the error message that VirtualiZarr raises here.

Can't guarantee anything but will try to get one started!

Also note our longer-term plan to make this easier for users is to change the defaults in xarray to be less of a bear trap: pydata/xarray#10062

This seems like a great direction to me :)

[danielfromearth]

Well, I started trying to work on a PR on the develop branch. I'm trying to catch up with the latest changes, but having trouble with the new API — running into an error when using the example usage code in the 'latest' docs (as of this writing) to open my locally saved netCDF:

from virtualizarr import open_virtual_dataset
from virtualizarr.parsers import HDFParser
from obstore.store import LocalStore

vds = open_virtual_dataset('dataset1.nc', object_store=LocalStore, parser=HDFParser())

Raises ValueError: Expected an object store instance, got <class 'obstore.store.LocalStore'>

Full traceback

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[12], line 5
      2 from virtualizarr.parsers import HDFParser
      3 from obstore.store import LocalStore
----> 5 vds = open_virtual_dataset('dataset1.nc', object_store=LocalStore, parser=HDFParser())

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/virtualizarr/xarray.py:49](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/virtualizarr/xarray.py#line=48), in open_virtual_dataset(file_url, object_store, parser, drop_variables, loadable_variables, decode_times, cftime_variables, indexes)
     37 def open_virtual_dataset(
     38     file_url: str,
     39     object_store: ObjectStore,
   (...)     45     indexes: Mapping[str, xr.Index] | None = None,
     46 ) -> xr.Dataset:
     47     filepath = validate_and_normalize_path_to_uri(file_url, fs_root=Path.cwd().as_uri())
---> 49     manifest_store = parser(
     50         file_url=filepath,
     51         object_store=object_store,
     52     )
     54     ds = manifest_store.to_virtual_dataset(
     55         loadable_variables=loadable_variables,
     56         decode_times=decode_times,
     57         indexes=indexes,
     58     )
     59     return ds.drop_vars(list(drop_variables or ()))

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/virtualizarr/parsers/hdf/hdf.py:143](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/virtualizarr/parsers/hdf/hdf.py#line=142), in HDFParser.__call__(self, file_url, object_store)
    138 def __call__(
    139     self,
    140     file_url: str,
    141     object_store: ObjectStore,
    142 ) -> ManifestStore:
--> 143     reader = ObstoreReader(store=object_store, path=file_url)
    144     manifest_group = _construct_manifest_group(
    145         filepath=file_url,
    146         reader=reader,
    147         group=self.group,
    148         drop_variables=self.drop_variables,
    149     )
    150     registry = ObjectStoreRegistry({get_store_prefix(file_url): object_store})

File [/srv/conda/envs/earthaccess/lib/python3.13/site-packages/virtualizarr/utils.py:38](https://openscapes.2i2c.cloud/srv/conda/envs/earthaccess/lib/python3.13/site-packages/virtualizarr/utils.py#line=37), in ObstoreReader.__init__(self, store, path)
     35 else:
     36     filepath = parsed.path
---> 38 self._reader = obs.open_reader(store, filepath)

ValueError: Expected an object store instance, got <class 'obstore.store.LocalStore'>

Perhaps I'm missing something quite simple?

[TomNicholas]

You need to instantiate the LocalStore before passing it in.

Thanks for being out guinea pig 😅

[danielfromearth]

Hah, of course. Happy to be the guinea pig here 😄