My user story is that I had up to some point a DataArray with index 'a', at some point code changed and I used as index 'b', but some code was still operating on 'a'; and I expected that code to error on accessing 'a'.

If I do np.diff(x,axis=(some non-existent axis)), numpy errors.

I might be generalizing, please correct me if I'm wrong, but afaik broadcasting happens in operations with 2 or more arrays, and the resulting dimensions are the union of the dimensions of each input array; the dimensions operated on are always at least in one of the involved arrays. Broadcasting does not happen if the arrays have the same shape, as it is for the diff operation.

So by following the same logic, when using an unary operator on a Dataset, an error should be thrown if none of the contained variables have the requested dimension. Since xarray allows scalar-sized dimensions (where the index shape is () ), those could be used to keep current behavior.

I read the linked issue. imho, the operation on a DataArray should error if the dimension is not there; on a Dataset the operation should be applied to each DataArray that contains the dimension, and if there are none, an error should be thrown. i.e. (set of dimensions allowed as arguments to operators on a set of containers) = ( union of the dimensions in each container).

This looks like a bug in my opinion...

@andersy005

This runs with not issues atm.

python with dask.config.set({"tokenize.ensure-deterministic":True}): ds = xr.tutorial.open_dataset('rasm') b = ds.isel(y=0) assert dask.base.tokenize(b) == dask.base.tokenize(b)

With: xarray 2022.3.0 pyhd8ed1ab_0 conda-forge dask 2022.5.0 pyhd8ed1ab_0 conda-forge

Often I run a function over a dataset, with each call outputing a hierarchical data structure, containing fixed dimensions in the best cases and variable length in the worst.

By "variable" length, do you mean that the length of dimensions differs between variables in the same group, or just that you don't know the length of the dimension in advance?

I mean that I might have, for instance, a map from 2 variables to data, ie (x,y)->c, that I can write as a DataArray XY with two dimensions x and y and the values being c. Then I have a function f so that f(c)->d[g(c)], i.e. it yields an array whose length depends on c. I wish I could say : apply f to XY, building a variable length array as you get the output. It could be stored as sparse matrice (X,Y,G). This is a bit out of scope for this discussion; but it is related since creating a differently named group per dimension length is often mentioned as a workaround ( which does not scale when you have a 1000x(variable length dimension) data).

Is there a specific use case which you think would require explicit dimensions to solve?

The use-case is iteratively adding values to a dataset by mapping functions over multiple variables / dimensions in arbitrary compositions. This happens in the context of data analysis, where you start with some source data and then iteratively create analysis functions, and then want to query / display / do statistics/reductions on the set of original data + analysis. Explicit hierarchical dimensions allow for merging and referring to data with no collisions in a single datatree/group.

PS: in netcdf-4 dimensions are seen by children, it matches what I previously posted; in HDF5 nodes are hardlinks to the actual data , this might be exactly the xarray-datagroup posted above.

Often I run a function over a dataset, with each call outputing a hierarchical data structure, containing fixed dimensions in the best cases and variable length in the worst. For this, it would make more sense to be able to have dimensions ( with optional labels and coordinates ) assigned to nodes (and these would be inherited by any descendants). Leaf nodes would hold data. On merge, dimensions could be bubbled up as long as length (and labels) matched. Operations with dimensions would then go down to corresponding dimension level before applying the operator, i.e. container['A/B'].mean('time') would be different from container['A'].mean('time')['B'].

Datagroup and Datatree are subcases of this general structure, which could be enforced via flags/checks. Option 1 is where the extremities of the tree are a node with two sets of child nodes, dimension labels and n-dimensional arrays. Option 2 is where the extremities of the tree are a node with a child node for a n-dimensional array A, and a sibling node for each dimension of A, containing the corresponding labels.

I'm sure I'm missing some big issue with the mental model I have, for instance I haven't thought of transformations at all and about coordinates. But for clarity I tried to write it down below.

The most general structure for a dataset I can think of is a directed graph. Each node A is a n-dimensional (sparse) array, where each dimension D points optionally to a one-dimensional node B with the same length.

To get a hierarchical structure, we:

• add edges of a different color, each with a label
• restrict their graph to a tree T
• add labels to each dimension D

We can resolve D's target by (A) checking for a sibling in T with the same name, and then going up one level and goto (A).

Multindexes ( multi-dimensional (sparse) labels ) generalize this model, but require tuple labels in T's edges i.e. : h/j/a[x,y,z] has a sybling h/j/(x,y)[x,y] , with z's labels being one level above, i.e. h/z[z] ( the notation a[b] means map of index b to value a ).

Any pointers regarding where to start / modules involved to implement this? I would like to have a try.

@keewis thanks, it is a duplicate. it must work in my specific case because I was testing with just one file.

I think I should have not expected a simple slice to make the choices that are necessary for evenly sampling a possibly irregular index. Assuming most users won't, I'm closing this. Below is what I'm using now. python def make_intervals(a,b,c): return a+np.arange(0,1+int(np.floor((b-a)/c)))*c print(d.sel(dim0=make_intervals(0,1,.1),method='nearest',tolerance=.01)) <xarray.DataArray (dim0: 11)> array([ 0, 10, 20, 30, 40, 49, 59, 69, 79, 89, 99]) Coordinates: * dim0 (dim0) float64 0.0 0.101 0.202 0.303 ... 0.697 0.798 0.899 1.0

@dcherian I tried to reproduce, with this minimal example I couldn't, so I'm closing the issue. ```python

import xarray as xr import numpy as np

n0 = 10 n1 = 3 x1 = xr.DataArray(np.empty((n0,n1),dtype=np.float64),dims=('dim0','dim1')).chunk({'dim0':2}) x2 = xr.DataArray(np.empty(n0,dtype=bool),dims=('dim0',)).chunk({'dim0':2})

n2 = 10

def f(x1,x2): return np.empty(n2,dtype=x1.dtype),np.empty(n2,dtype=np.min_scalar_type(n2))

m,w = xr.apply_ufunc( f, x1,x2, input_core_dims=[('dim0','dim1'),('dim0',)], output_core_dims=[('dim2',),('dim2',)], vectorize=True, dask='parallelized', dask_gufunc_kwargs={ 'output_sizes':{'dim2':n2}, 'allow_rechunk':True,

'meta':(np.empty((M,),dtype=p.dtype),np.empty((M,),dtype=np.min_scalar_type(M)))

'output_dtypes':[p.dtype,np.min_scalar_type(M)],

},
output_dtypes=[x1.dtype,np.min_scalar_type(n2)] # now works

output_dtypes=(x.dtype,np.min_scalar_type(ny)) # works

) m.compute(),w.compute() (<xarray.DataArray (dim2: 10)> array([1.e-323, 2.e-323, 3.e-323, 4.e-323, 5.e-323, 6.e-323, 7.e-323, 8.e-323, 9.e-323, 1.e-322]) Dimensions without coordinates: dim2, <xarray.DataArray (dim2: 10)> array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint8) Dimensions without coordinates: dim2) ```

I did a ctrl-f for zarr in this issue, found nothing, so here's my two cents: it should be possible to write a Datagroup with either zarr or netcdf. I wonder if @emilbiju (posted https://github.com/pydata/xarray/issues/4118 ) has any of that code laying around, could be a starting point. In general, a tree structure to which I can broadcast operations in the same dimensions to different datasets that do not necessary share dimension lengths would solve my use case . This corresponds to bullet point number 3 in https://github.com/pydata/xarray/issues/1092#issuecomment-290159834 . My use case is a set of experiments, that have: the same parameter variables, with different values; the same dimensions with different lengths for their data. The parameters and data would benefit of having a hierarchical naming structure. Currently I build a master dataset containing experiment datasets, with a coordinate for each parameter. Then I map functions over it.

I arrived here due to a different use case / problem, which ultimately I solved, but I think there's value in documenting it here. My use case is the following workflow: 1 . take raw data, build a dataset, append it to a zarr store Z 2 . analyze the data on Z, then maybe goto 1. Step 2's performance is much better when data on Z is chunked properly along the appending dimension 'frame' (chunks of size 50), however step 1 only adds 1 element along it. I end up with Z having chunks (1,1,1,1,1...) on 'frame'. On xarray 0.16.0, this seems solvable via the encoding parameter, if we take care to only use it on the store creation. Before that version, I was using something like the monkey patch posted by @chrisbarber . Code: ```python import shutil import xarray as xr import numpy as np import tempfile zarr_path = tempfile.mkdtemp()

def append_test(ds,chunks): shutil.rmtree(zarr_path)

for i in range(21):
    d = ds.isel(frame=slice(i,i+1))
    d = d.chunk(chunks)
    d.to_zarr(zarr_path,consolidated=True,**(dict(mode='a',append_dim='frame') if i>0 else {}))
dsa = xr.open_zarr(str(zarr_path),consolidated=True)
print(dsa.chunks,dsa.dims)

sometime before 0.16.0

import contextlib @contextlib.contextmanager def change_determine_zarr_chunks(chunks): orig_determine_zarr_chunks = xr.backends.zarr._determine_zarr_chunks try: def new_determine_zarr_chunks( enc_chunks, var_chunks, ndim, name): da = ds[name] zchunks = tuple(chunks[dim] if (dim in chunks and chunks[dim] is not None) else da.shape[i] for i,dim in enumerate(da.dims)) return zchunks xr.backends.zarr._determine_zarr_chunks = new_determine_zarr_chunks yield finally: xr.backends.zarr._determine_zarr_chunks = orig_determine_zarr_chunks chunks = {'frame':10,'other':50} ds = xr.Dataset({'data':xr.DataArray(data=np.random.rand(100,100),dims=('frame','other'))})

append_test(ds,chunks) with change_determine_zarr_chunks(chunks): append_test(ds,chunks)

with 0.16.0

def append_test_encoding(ds,chunks): shutil.rmtree(zarr_path)

encoding = {}
for k,v in ds.variables.items():
    encoding[k]={'chunks':tuple(chunks[dk] if dk in chunks else v.shape[i] for i,dk in enumerate(v.dims))}

for i in range(21):
    d = ds.isel(frame=slice(i,i+1))
    d = d.chunk(chunks)
    d.to_zarr(zarr_path,consolidated=True,**(dict(mode='a',append_dim='frame') if i>0 else dict(encoding = encoding)))
dsa = xr.open_zarr(str(zarr_path),consolidated=True)
print(dsa.chunks,dsa.dims)

append_test_encoding(ds,chunks) ```

Frozen(SortedKeysDict({'frame': (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), 'other': (50, 50)})) Frozen(SortedKeysDict({'frame': 21, 'other': 100})) Frozen(SortedKeysDict({'frame': (10, 10, 1), 'other': (50, 50)})) Frozen(SortedKeysDict({'frame': 21, 'other': 100})) Frozen(SortedKeysDict({'frame': (10, 10, 1), 'other': (50, 50)})) Frozen(SortedKeysDict({'frame': 21, 'other': 100}))