Understood. Thank you for your prompt replies.

I'll read up on ask again if I have any questions.

I guess I was trying to accommodate in the past users that were not using our wrappers to to_netcdf

nice!

I know it is "stale" but aligning to these "surprise dimensions" creates "late stage" bugs that are hard to pinpoint.

I'm not sure if it is possible to mark these dimensions as "unnamed" and as such, they should be "merged" into new "unnamed" dimensions that the user isn't tracking at this point in time.

Our workaround have included calling these dimensions something related to the datarray d1_i, or simply making small small "arrays" a countable number of scalar variables (d1_min, d1_max) instead of a single array containing two values d1_limits.

```python import xarray as xr d1 = xr.DataArray(data=[1, 2]) assert 'dim_0' in d1.dims d2 = xr.DataArray(data=[1, 2, 3]) assert 'dim_0' in d2.dims

xr.Dataset({'d1': d1, 'd2': d2}) ```

cc: @claydugo

I kinda think this reminds me of

https://github.com/pydata/xarray/discussions/7359

As a followup question, is the LazilyIndexedArray part of the 'public api'. That is when you do decide to refactor, https://docs.xarray.dev/en/stable/generated/xarray.core.indexing.LazilyIndexedArray.html

Will you try to warn us users that choose to

from xarray.core.indexing import LazilyIndexedArray

I'm going to say, the LazilyIndexedArray is pretty cool.

Ultimately, I'm not sure how you want to manage resources. This zarr store could be considered a resource and thus, may have an owner. Or maybe zarr should close itself upon garbage cleanup.

I'm not sure. I decided not to use zarr (not now) so i lost interest sorry.

Thank you for your explination.

Do you think it is safe to "strip" encoding after "loading" the data? or is it still used after the initial call to open_dataset?

Kind bump

Any chance of a release, this is quite breaking for large datasets that can only be out of memory.

👍🏾

Yes, without chunks of anything

Very smart test!

The data is loaded from an NetCDF store through open_dataset

No explicit test was added to ensure that the data wasn't loaded. I just experienced this bug enough (we would accidentally load 100GB files in our code base) that I knew exactly how to fix it.

If you want i can add a test to ensure that future optimizations to nbytes do not trigger a data load.

I was hoping the 1 line fix would be a shoe in.

It seems that checking hasattr on the _data variable achieves both purposes.

I think that at the very lease, the current implementation works as well as the old one for arrays that are defined by the sparse package.

Looking into the history a little more. I seem to be proposing to revert: https://github.com/pydata/xarray/commit/60f8c3d3488d377b0b21009422c6121e1c8f1f70

I think this is important since many users have arrays that are larger than memory. For me, I found this bug when trying to access the number of bytes in a 16GB dataset that I'm trying to load on my wimpy laptop. Not fun to start swapping. I feel like others might be hitting this too.

xref: https://github.com/pydata/xarray/pull/6797 https://github.com/pydata/xarray/issues/4842

I personally do not even think the hasattr is really that useful. You might as well use size and itemsize

I'm really not sure. It seems to happen with a large swath of versions from my recent search.

Also running from the python REPL, i don't see the warning. which makes me feel like numpy/cython/netcdf4 are trying to suppress the harmless warning.

https://github.com/cython/cython/blob/0.29.x/Cython/Utility/ImportExport.c#L365

one or the other.

I think it is a numpy thing bash mamba create --name np numpy netcdf4 --channel conda-forge --override-channels conda activate np python -c "import numpy; import warnings; warnings.filterwarnings('error'); import netCDF4" python Traceback (most recent call last): File "<string>", line 1, in <module> File "/home/mark/mambaforge/envs/np/lib/python3.11/site-packages/netCDF4/__init__.py", line 3, in <module> from ._netCDF4 import * File "src/netCDF4/_netCDF4.pyx", line 1, in init netCDF4._netCDF4 RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility. Expected 16 from C header, got 96 from PyObject

mamba create --name xr numpy pandas xarray netcdf4 --channel conda-forge --override-channels conda activate xr python -c "import xarray; import warnings; warnings.filterwarnings('error'); import netCDF4"

Traceback (most recent call last): File "<string>", line 1, in <module> File "/home/mark/mambaforge/envs/xr/lib/python3.11/site-packages/netCDF4/__init__.py", line 3, in <module> from ._netCDF4 import * File "src/netCDF4/_netCDF4.pyx", line 1, in init netCDF4._netCDF4 RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility. Expected 16 from C header, got 96 from PyObject

It is also reproducible on binder:

It seems that the binder uses conda-forge, which is why i'm commenting here.

It is really strange in the sense that xarray doesn't compile anything.

https://github.com/conda-forge/xarray-feedstock/blob/main/recipe/meta.yaml#L16

So it must be something that gets lazy loaded that triggers things.

As I've been recently going down this performance rabbit hole, I think the discussion around https://github.com/pydata/xarray/issues/7045 is relevant and provides some additional historical context as to "why" this performance penalty might be happening.

While the above "fix" addresses the issues with renaming coordinates, I think there are plenty of usecases where we would still end up with strange, or unexpected results. For example.

1. Load a dataset with many non-indexing coordinates.
2. Dropping variables (that happen to be coordinates).
3. Then adding back a variable with the same name.
4. Upon save, encoding would dictate that it is a coordinate of a particular variable and will promote it to a coordinate instead of data.

We could apply the "fix" to the drop_vars method as well, but I think it may be hard (though not impossible) to hit all the cases.

I think a more "generic", albeit breaking" fix would be to remove the "coordinates" entirely from encoding after the dataset has been loaded. That said, this only "works" if dataset['variable_name'].encoding['coordinates'] is considered a private variable. That is, users are not supposed to be adding to it at will.

And if you want to have a clean encoding dictionary, you may want to do the following:

python names = set(names) for _, variable in obj._variables.items(): if 'coordinates' in variable.encoding: coords_in_encoding = set(variable.encoding.get('coordinates').split(' ')) remaining_coords = coords_in_encoding - names if len(remaining_coords) == 0: del variable.encoding['coordinates'] else: variable.encoding['coordinates'] = ' '.join(remaining_coords)

I think these are good alternatives.

From my experiments (and I'm still trying to create a minimum reproducible code that shows the real problem behind the slowdowns) reindexing can be quite an expensive. We used to have many coordinates (to ensure that critical metdata stays with data_variables) and those coordinates were causing slowdowns on reindexing operations.

Thus the two calls update and expand_dims might cause two reindex merges to occur.

However, for this particular issue, I think that documenting the strategies proposed in the docstring is good enough. I have a feeling if one can get to the bottom of 7224, the performance concerns here will be mitigated too.

We can leave the performance discussion to: https://github.com/pydata/xarray/issues/7224

Hmm...I was kind of hoping we could avoid something like adding a _stacklevel_increment argument.

Right. thank you for finding that example. I was going to try to construct one.

xref: https://github.com/pandas-dev/pandas/pull/49393

Ok, I don't think I have the right tools to really get to the bottom of this. The spyder profiler just seems to slowdown code too much. Any other tools to recommend?

Well now the benchmarks look like they make more sense: [ 50.00%] ··· ==================== ========== ========== ========== ========== ========== -- count -------------------- ------------------------------------------------------ strategy 0 1 10 100 1000 ==================== ========== ========== ========== ========== ========== dict_of_DataArrays 1.56±0ms 3.60±0ms 5.83±0ms 16.6±0ms 67.3±0ms dict_of_Variables 1.65±0ms 3.11±0ms 4.03±0ms 6.22±0ms 18.9±0ms dict_of_Tuples 2.42±0ms 3.11±0ms 982±0μs 5.17±0ms 17.2±0ms ==================== ========== ========== ========== ========== ==========

[ 50.00%] ··· ==================== ========== ========== ========== ========== ========== -- count -------------------- ------------------------------------------------------ strategy 0 1 10 100 1000 ==================== ========== ========== ========== ========== ========== dict_of_DataArrays 1.65±0ms 3.83±0ms 4.03±0ms 6.14±0ms 16.6±0ms dict_of_Variables 3.04±0ms 3.24±0ms 3.38±0ms 4.04±0ms 9.91±0ms dict_of_Tuples 2.90±0ms 3.03±0ms 3.32±0ms 3.22±0ms 3.22±0ms ==================== ========== ========== ========== ========== ==========

as you though, the numbers improve quite a bit.

I kinda want to understand why a no-op takes 1 ms! ^_^

What about just specifying "dims"?

With the right window size it looks like: [ 50.00%] ··· ==================== ========== ========== ========== ========== ========== -- count -------------------- ------------------------------------------------------ strategy 0 1 10 100 1000 ==================== ========== ========== ========== ========== ========== dict_of_DataArrays 1.32±0ms 5.87±0ms 7.58±0ms 18.7±0ms 98.6±0ms dict_of_Variables 2.70±0ms 2.91±0ms 3.01±0ms 3.91±0ms 7.04±0ms dict_of_Tuples 2.84±0ms 3.02±0ms 3.22±0ms 3.42±0ms 3.02±0ms ==================== ========== ========== ========== ========== ==========

On the CI, it reports similar findings: ``` [ 67.73%] ··· ...dVariable.time_dict_of_dataarrays_to_dataset ok [ 67.73%] ··· =================== ============= existing_elements
------------------- ------------- 0 269±0.9μs
10 2.21±0.01ms 100 16.5±0.07ms 1000 153±0.9ms
=================== =============

[ 67.88%] ··· ...etAddVariable.time_dict_of_tuples_to_dataset ok [ 67.88%] ··· =================== =========== existing_elements
------------------- ----------- 0 269±0.6μs 10 289±0.4μs 100 293±1μs
1000 346±0.4μs =================== ===========

[ 68.02%] ··· ...ddVariable.time_dict_of_variables_to_dataset ok [ 68.02%] ··· =================== ============= existing_elements
------------------- ------------- 0 270±1μs
10 329±0.6μs
100 636±1μs
1000 3.70±0.01ms =================== =============

[ 68.17%] ··· ...e.DatasetAddVariable.time_merge_two_datasets ok [ 68.17%] ··· =================== ============= existing_elements
------------------- ------------- 0 104±0.5μs
10 235±0.6μs
100 1.05±0ms
1000 9.02±0.02ms =================== =============

[ 68.31%] ··· ...e.DatasetAddVariable.time_variable_insertion ok [ 68.31%] ··· =================== ============= existing_elements
------------------- ------------- 0 119±1μs
10 225±0.7μs
100 1.04±0ms
1000 9.03±0.03ms =================== ============= ```

``` $ asv run -E existing --quick --bench merge · Discovering benchmarks · Running 5 total benchmarks (1 commits * 1 environments * 5 benchmarks) [ 0.00%] ·· Benchmarking existing-py_home_mark_mambaforge_envs_mcam_dev_bin_python [ 10.00%] ··· merge.DatasetAddVariable.time_dict_of_dataarrays_to_dataset ok [ 10.00%] ··· =================== ========== existing_elements ------------------- ---------- 0 762±0μs 10 7.18±0ms 100 12.6±0ms 1000 89.1±0ms =================== ==========

[ 20.00%] ··· merge.DatasetAddVariable.time_dict_of_tuples_to_dataset ok [ 20.00%] ··· =================== ========== existing_elements ------------------- ---------- 0 889±0μs 10 2.01±0ms 100 1.34±0ms 1000 605±0μs =================== ==========

[ 30.00%] ··· merge.DatasetAddVariable.time_dict_of_variables_to_dataset ok [ 30.00%] ··· =================== ========== existing_elements ------------------- ---------- 0 2.48±0ms 10 2.06±0ms 100 2.13±0ms 1000 2.38±0ms =================== ==========

[ 40.00%] ··· merge.DatasetAddVariable.time_merge_two_datasets ok [ 40.00%] ··· =================== ========== existing_elements ------------------- ---------- 0 814±0μs 10 945±0μs 100 2.42±0ms 1000 5.23±0ms =================== ==========

[ 50.00%] ··· merge.DatasetAddVariable.time_variable_insertion ok [ 50.00%] ··· =================== ========== existing_elements ------------------- ---------- 0 1.10±0ms 10 954±0μs 100 1.88±0ms 1000 5.29±0ms =================== ========== ```

Exciting improvements on usability for the next version!

but that would be a lot of work especially for such a critical piece of code in Xarray.

Agreed. I'll take the small wins where I can :D.

Great! I think this will be a good addition with: https://github.com/pydata/xarray/pull/7223#discussion_r1007023769

ok naming is always hard. I tried to pick a good name.

I know it is not comparable, but I was really curious what "dictionary insertion" costs, in order to be able to understand if my comparisons were fair:

I think xarray gives me 3 order of magnitude of "thinking" benefit, so I'll take it! python --version Python 3.9.13

Ok. I'll want to rethink them.

I know it looks quadratic time, but i really would like to test n=1000 and i have an idea

I'm somewhat ocnfused, I can run the benchmark locally

``` [ 1.80%] ··· dataset_creation.Creation.time_dataset_creation 4.37±0s

```

:/ not fun, the benchmark is failing. not sure why.

I think the rapid return, helps by about 40% is still pretty good.

hmm ok. it seems i can't blatently avoid the copy like that.

Would be interesting to see whether this was covered by our existing asv benchmarks.

I wasn't able to find something that really benchmarked "large" datasets.

Would be a good benchmark to add if we don't have one already.

Added one.

The reason this is a separate merge request, is that I agree that this is more contentious as a change.

However, I will argue that Aligner should really not be a class.

Using ripgrep you find that the only instances of Aligner exist internally: ` xarray/core/dataset.py 2775: aligner: alignment.Aligner, 2783: """Callback called fromAligner`` to create a new reindexed Dataset."""

xarray/core/alignment.py 107:class Aligner(Generic[DataAlignable]): 114: aligner = Aligner(objects, *kwargs) <------- Example 767: aligner = Aligner( <----------- Used and consumed for the method align 881: aligner = Aligner( <----------- Used and consumed for the method reindex 909: # This check is not performed in Aligner.

xarray/core/dataarray.py 1752: aligner: alignment.Aligner, 1760: """Callback called from Aligner to create a new reindexed DataArray.""" ```

out of interest, how did you find this?

Spyder profiler

Thank you all for taking the time to study, and worry about these improvements.

Now i have to figure out how my software went from 2 sec loading time to 12 ;) Totally unrelated to this. But one day I'll have benchmarking in place to monitor it :D.

Separate issue, but do these need to be imported into xarray/init.py

At this point removing testing and tutorial would be strange and break things. Stefan in the discussion linked above speaks about the reasoning behind importing submodules in the top level namespace.

In developing https://github.com/pydata/xarray/pull/7172, there are also some places where class types are used to check for features: https://github.com/pydata/xarray/blob/main/xarray/core/pycompat.py#L35

Dask and sparse and big contributors due to their need to resolve the class name in question.

Ultimately. I think it is important to maybe constrain the problem.

Are we ok with 100 ms over numpy + pandas? 20 ms?

On my machines, the 0.5 s that xarray is close to seems long... but everytime I look at it, it seems to "just be a python problem".

I guess the code: ```python import xarray as xr dataset = xr.Dataset()

my_variable = np.asarray(dataset.get('my_variable', np.asarray(1.0))) ``` coerces things as an array.

Talking things out made me find this one. Though it doesn't read very well.

Feel free to close.

Are the functions you are considering using this functions that never had keyword arguments before? When I wrote a similar decorator before i had an explicit list of arguments that were allowed to be converted.

These decorators are kinda fun to write and are quite taylored to a certain release philosophy.

It might be warranted to just write your own ;)

Interesting.

I think in my readme i suggest vedoring the code.

Happy to give you a license for it so you don't need to credit me in addition to your own license.

@max-sixty unfortunately, I think the way hdf5 is designed, it doesn't try to be too smart about what would be the best fine tuning for your particular system. In some ways, this is the correct approach.

The current constructor pathway: https://github.com/pydata/xarray/blob/main/xarray/backends/h5netcdf_.py#L164

Doesn't provide a user with a catchall-kwargs. I think this would be an acceptable solution.

I should say that the the performance of the direct driver is terrible without aligned data: https://github.com/Unidata/netcdf-c/pull/2206#issuecomment-1054855769

I have to elaborate that this may be even more important for users that READ the data back alot. Reading with the standard Xarray operands hits other limits, but one limit that it definitely hits is that of the HDF5 driver used.

Thanks for merging so quickly

Thank you @kmuehlbauer for the explicit PR link.

I do plan on adding alignment features to h5py then to bring it toward h5netcdf. So I think something like this will be useful in the future.

Feature request link: https://github.com/h5py/h5py/issues/2034

ImportError is a superset of ModuleNotFoundError. https://github.com/python/cpython/blob/f4c03484da59049eb62a9bf7777b963e2267d187/Lib/test/exception_hierarchy.txt#L19

So it depends what question you care about asking:

1. Does the python exist? You should test for ModuleNotFoundError
2. Is the package usable? You should probably test for ImportError

I think question 2 is friendlier to xarray users.

This is likely true. Thanks for looking back into this.

I hope the following can help users that struggle with the speed of xarray:

I've found that when doing numerical computation, I often use the xarray to grab all the metadata relevant to my computation. Scale, chromaticity, experimental information.

Eventually, i create a function that acts as a barrier: - Xarray input (high level experimental data) - Computation parameters output (low level implementation detail relevant information).

The low level implementation can operate on the fast numpy arrays. I've found this to be the struggle with creating high level APIs that do things like sanitize inputs (xarray routines like _validate_indexers and _broadcast_indexes) and low level APIs that are simply interested in moving and computing data.

For the example that @nbren12 brought up originally, it might be better to create xarray routines (if they don't exist already) that can create fast iterators for the underlying numpy arrays given a set of dimensions that the user cares about.

I think you should be able to define your own custom encoder if you want it to be a datetime.

But inevitably, you will have to define your own save and load functions.

Python, by definition of being such a loose language, allows you to do things that the original developers never really imagined.

this can sometimes lead to silent corruption.like the one you've experienced.

It isn't really part of any library. I don't really have plans of making it into a public library. I think the discussion is really around the xarray API, and what functions to implement at first.

Then somebody can take the code and integrate it into the decided upon API.

yeah, i'm not too sure. I think the idea is that this breaks compatibility with netcdf times, so the resulting file is thus not standard.

For my application, us timing is enough.

Small prototype, but maybe it can help boost the development.

i ran into this problem trying to round trip time to the nanosecond (even though i don't need it, sub micro second would be nice)

but unfrotunately, you run into the fact that cftime doesn't support nanoseconds https://github.com/Unidata/cftime/blob/master/cftime/_cftime.pyx

Seems like they discussed a nanosecond issue a while back too https://github.com/Unidata/cftime/issues/77

Their ultimate point was that there was little point in having precision down to the nano second given that python datetime objects only have microseconds. I guess they are right.

I think I got a basic prototype working.

That said, I think a real challenge lies in supporting the numerous backends and lazy arrays.

For example, I was only able to add data in peculiar fashions using the netcdf4 library which may trigger complex computations many times.

Is this a use case that we must optimize for now?

I'm not too sure about this anymore.

with the way the test is written now, it is unclear to me if the store should be closed afterward.

I'm also unsure of how to deal with the case where the user passed it a ZipStore instead of a string. Will have to keep thinking.

Sorry for noise. It seems that 1D arrays are still supported.

I still had a 2D array lingering in my codebase.

With the netcdf4 back end, I'm not able to save a 1D attr dataset.

I can save my dataset with the h5netcdf backend

Sorry for posting on such an old thread.

Are attrs supposed to support 1D arrays?

xfail just gets forgotten, so i'll leave it for now.

And actually, zarr provides a data argument in create_dataset that actually encounters the same bug python import zarr import numpy as np name = 'hello' data = np.array('world', dtype='<U5') store = zarr.ZipStore('test_store.zip', mode='w') root = zarr.open(store , mode='w') zarr_array = root.create_dataset(name, data=data, shape=data.shape, dtype=data.dtype) zarr_array[...]

I guess i can open upstream in zarr, but I think for catching the 0 sized array case, it is probably best to use the data argument.

I think the reason is that for zero sized arrays, you technically aren't allowed to write data to them.

https://github.com/pydata/xarray/blob/6378a711d50ba7f1ba9b2a451d4d1f5e1fb37353/xarray/backends/zarr.py#L449

This means that when you create the 0 sized array, you can't actually change the value.

Here is a reproducer without xarray python import zarr import numpy as np name = 'hello' data = np.array('world', dtype='<U5') store = zarr.ZipStore('test_store.zip', mode='w') root = zarr.open(store , mode='w') zarr_array = root.create_dataset(name, shape=data.shape, dtype=data.dtype) root[name][...] = data zarr_array[...]

Though the code path follows what xarray does in the backend.

Alright, it probably makes more sense to reopen this when the issue gets fixed.

wouldn't this be a useful to test to have? I think the ability to save things in a zip store is quite useful.

@jakirkham not sure if you have any thoughts on why the code above is bugging out.

Honestly, i've found that __ coordinates are treated as less strict coordinates, so I've been using those.

Keeping things like: * The version of different software that was used.

Seems more like an attribute, but really, it is all data, so :/.

Regarding the originally issue, it seems that you are right in the sense that a 0 dimension string might be buggy in zarr itself.

I guess we (when we have time) will have to dig down to find a MVC example that reproduces the issue without xarray to submit to Zarr.

Hmm, interesting!

I've avoided attrs since they often get "lost" in computation, and don't get dragged along as rigorously as coordinates.

I do have some real coordinates that are stored as strings.

Thanks for the quickfeedback.

Here is the reproducing code without using context managers (which auto clsoe things you know)

```python import xarray as xr import zarr x = xr.Dataset() x['hello'] = 'world' x with zarr.ZipStore('test_store.zip', mode='w') as store: x.to_zarr(store)

read_store = zarr.ZipStore('test_store.zip', mode='r')
x_read = xr.open_zarr(read_store).compute()

The error will happen before this line is executed

read_store.close()

```

I will have to try the debugging things you mentionned some later time :/

Not sure if the builds in https://github.com/pydata/xarray/pull/3888 help reproduce things or not?

hmm i didn't realize this. I"m running from conda-forge + linux.

Let me try on your CIs.

See the zipstore example in my first comment

I thought I provided it, but in either case, here is my traceback

My guess is that that xarray might be trying to write to the store character by character???

Otherwise, not too sure.

Sure, I just wanted to make the note that this operation should be more or less constant time, as opposed to dependent on the size of the array. Somebody had mentionned it should increase with the size of the array.

Hmm, slicing should basically be a no-op.

The fact that xarray makes it about 100x slower is a real killer. It seems from this conversation that it might be hard to workaround

mind blown!!!! thanks for that pointer I haven't touched my serialization code in a while, kinda scared to go back to it now, but I will keep that library in mind.

I saw Zarr a while back, looks cool. I hope to see it grow.

no need to be sorry. These two functions were easy enough for me to do myself in my own codebase.

There are few issues that I've found doing this though. Mainly, I can't find a good way to serialize numpy arrays in a round-trippable fashion. It is difficult to get back lists of arrays, or arrays of unit8. I don't know if you have a good way to solvle this problem.

@shoyer, @fmaussion thank you for your answers. I'm OK with this issue being closed. I'm no expert on netcdf4, so I don't know if I could express the issue in a concise manner there.

Thanks!

Maybe the issue that I am facing is that I want to deal with the storage of my metadata and data seperately.

I used to have my own library that was replicating much of xarray's functionality, but your code is much nicer than anything I would be able to write in a finite time. :smile:

Following the information here: http://xarray.pydata.org/en/stable/data-structures.html#coordinates-methods

Currently, my serialization pipeline is: ```python import xarray as xr import numpy as np

Setup an array with coordinates

n = np.zeros(3) coords={'x': np.arange(3)} m = xr.DataArray(n, dims=['x'], coords=coords)

coords_dataset_dict = m.coords.to_dataset().to_dict() coords_dict = coords_dataset_dict['coords']

Read/Write dictionary to JSON file

This works, but I'm essentially creating an emtpy dataset for it

coords_set = xr.Dataset.from_dict(coords_dataset_dict) coords2 = coords_set.coords # so many coords :D m2 = xr.DataArray(np.zeros(shape=m.shape), dims=m.dims, coords=coords2)

I used to just pass the dataset to "coords"

m3 = xr.DataArray(np.zeros(shape=m.shape), dims=m.dims, coords=coords_set) ```

Is there a better way to serialize coordinates only?

It seems like this warning isn't benign though. I will take your suggestion though (coords=dataset.coords)

I feel like I'm not the only one who probably did this. Should you raise an other warning explicitly?

I came across this when serializing/deserializing my coordinates to a json file.

Finishing up this line of though: without the assumption that the relative order of dimensions is maintained across arrays in a set, this feature is impossible to implement as a neat function call. You would have to specify exactly how to expand each of the coordinates which can get pretty long.

I wrote some code, that I think should have worked if relative ordering was a valid assumption: Here it is for reference https://github.com/hmaarrfk/xarray/pull/1

To obtain the desired effect, you have to expand the dimensions of the coordinates individually: ```python import xarray as xr import numpy as np

Setup an array with coordinates

n = np.arange(1, 13).reshape(3, 2, 2) coords={'y': np.arange(1, 4), 'x': np.arange(1, 3), 'xi': np.arange(2)}

%%

z = xr.DataArray(n[..., 0]2, dims=['y', 'x']) a = xr.DataArray(n, dims=['y', 'x', 'xi'], coords={*coords, 'z': z})

sliced = a[0] print("The original xarray") print(a.z) print("The sliced xarray") print(sliced.z)

%%

expanded = sliced.expand_dims('y', 0) expanded['z'] = expanded.z.expand_dims('y', 0) print(expanded) ```