Is this change still of potential value?

In my original PR, I wrote:

In putting this together, I noted that open_zarr to re-read the data triggers this bug, while open_dataset(..., engine='zarr') does not. I'm not sure if my proposed fix is a band-aid, or if something in open_zarr is the real culprit.

I just want to be sure we think this is the correct fix, and there won't be any unintended consequences. I'm really not familiar with the interaction of chunks between zarr, dask, and xarray. I don't feel 100% sure this is the correct solution, though it fixes the immediate solution.

@jhamman Ready for review, though not urgent.

Mostly, I'm curious if this the right way to go about solving the linked issue, or if there is something deeper in xarray to update.

Sure thing. I'll revive it.

Is this solution desired, or is there a deeper fix?

It may run even deeper -- there seem to be several checks on dimension sizes that would need special casing. Even simply doing a variable[dim] lookup fails!

There seems to be some support, but now you have me worried. I have a used xarray mainly for labelling, but not for much computation -- I'm dropping into dask because I need map_overlap.

FWIW, calling dask.compute(arr) works with unknown chunk sizes, but now I see arr.compute() does not. This fooled me into thinking I could use unknown chunk sizes. Now I see that writing to zarr does not work, either. This might torpedo my current design.

I see the compute_chunk_sizes method, but that seems to trigger computation. I'm running on a dask cluster -- is there anything I can do to salvage the pattern arr_with_nan_shape.to_dataset().to_zarr(compute=False) (with our without xarray)?

I do not see the failures in my client that are seen in the checks, so there must be some mismatch in the environments.

In putting this together, I noted that open_zarr to re-read the data triggers this bug, while open_dataset(..., engine='zarr') does not. I'm not sure if my proposed fix is a band-aid, or if something in open_zarr is the real culprit.

My mental model of what's happening may not be correct. I did want sel(), isel(), and squeeze() to all operate the same way (and maybe someday even work on non-dim coordinates!). Replacing squeeze() with isel() in my initial example gives the same failure, which I would want it to work:

``` import numpy as np import xarray as xr

arr1 = xr.DataArray(np.zeros((1,5)), dims=['y', 'x'], coords={'e': ('y', [10])}) arr2 = arr1.isel(y=0).expand_dims('y') xr.testing.assert_identical(arr1, arr2) ```

``` AssertionError: Left and right DataArray objects are not identical

Differing coordinates: L e (y) int64 10 R e int64 10 ```

The non-dim coordinate e has forgotten that it was associated with y. I'd prefer that this association remained.

Where it gets really interesting is in the following example where the non-dim coordinate moves from one dim to another. I understand the logic here (since the isel() were done in a way that correlates 'y' and 'z'). In my proposal, this would not happen without explicit user intervention -- which may actually be desired here (it's sorta surprising):

``` import numpy as np import xarray as xr

arr = xr.DataArray(np.zeros((2, 2, 5)), dims=['z', 'y', 'x'], coords={'e': ('y', [10, 20])}) print(arr.coords) print()

arr0 = arr.isel(z=0,y=0) arr1 = arr.isel(z=1,y=1)

arr_concat = xr.concat([arr0, arr1], 'z') print(arr_concat.coords) ```

``` Coordinates: e (y) int64 10 20

Coordinates: e (z) int64 10 20 ```

One problem with this -- at least for now -- is that xarray currently doesn't allow coordinates on DataArray objects to have dimensions that don't appear on the DataArray itself.

Ah, then that would be the desire here.

It might also be surprising that this would make squeeze('y') inconsistent with isel(y=0)

The suggestion here is that both of these would behave the same. The MCVE was just for the squeeze case, but I expect that isel and sel would both allow non-dim coords to maintain the reference to their original dim (even if it becomes a non-dim coord itself).

I'm not a huge fan of adding arguments for a case that rarely comes up (I presume). 

One difference in your example is that the 'e' coord is never based on 'y', so I would not want it expanded -- so I'd still like that test to pass.

The case I'm interested in is where the non-dimension coords are based on existing dimension coords that gets squeezed.

So in this example: ``` import numpy as np import xarray as xr

arr1 = xr.DataArray(np.zeros((1,5)), dims=['y', 'x'], coords={'y': [42], 'e': ('y', [10])}) arr1.squeeze() ```

The squeezed array looks like:

<xarray.DataArray (x: 5)> array([0., 0., 0., 0., 0.]) Coordinates: y int64 42 e int64 10 Dimensions without coordinates: x

What I think would be more useful:

<xarray.DataArray (x: 5)> array([0., 0., 0., 0., 0.]) Coordinates: y int64 42 e (y) int64 10 <---- Note the (y) Dimensions without coordinates: x

There has been some discussion on the dask chunking issue here: https://github.com/dask/dask/issues/3650 https://github.com/dask/dask/issues/5544

Regarding the position of the inserted variable, it is not related to the chunking. It seems possible to do this. Would this be an acceptable change? If so, the first problem is the signature, as multiple dimensions may be passed in. :/

Hi! This work is interesting to me, as I was implementing in dask an image processing algo which needs an intermediate 1-d linear interpolation step. This bottlenecks the calculation through a single node. Your work here on distributed interpolation is intriguing, and I'm wondering if it would be useful in my work and if it could possibly become part of dask itself.

Here is the particular function, which you'll note has a dask.delayed wrapper around np.interp.

What about the case of no missing values, when other wouldn't be needed? Could the same dtype be returned then? This is my case, since I'm re-purposing where to do sel for non-dimension coordinates.

Could you give a concrete example of what this would look like?

It seems rather unlikely to me to have an example of where with drop=True where the condition is exactly aligned with the grid, such that there are no missing values.

I guess it could happen if you're trying to index out exactly one element along a dimension?

That's exactly right. I am just selecting one slice of a data array, using data.where(data.coords['stain'] == 'DAPI').

In the long term, the cleaner solution for this will be some form for support for more flexibly / multi-dimensional indexing.

Agreed. Once I actually get things running, I'll be ready to try and contribute fixes for all my TODOs that reference xarray github issues. :)

What about the case of no missing values, when other wouldn't be needed? Could the same dtype be returned then? This is my case, since I'm re-purposing where to do sel for non-dimension coordinates.

I'm capable of just recasting for my use case, if this is becoming an idea that would be difficult to maintain/document.

If drop=True, would it be problematic to return the same dtype or allow other?

My use case is a simple slicing of a dataset -- no missing values. The use of where is due to one of selections being on a non-dimension coordinate (#2028).

I can workaround using astype, but will say I was mildly surprised by this feature. I now understand why it's there. Our code is old and the data is intermediate and never deeply inspected -- I only noticed this when we started using a memory-intensive algorithm and surprised how much space was taken by our supposed uint16 data. :)

If there is a non-dimension coordinate, the error is also tickled.

``` import xarray as xr import numpy as np ds=xr.Dataset({'foo': (['bar'], np.zeros((100,)))})

Problem also affects non-dimension coords

ds.coords['baz'] = ('bar', ['mah'] * 100)

ds.to_zarr('test.zarr', mode='w') ds2 = xr.open_zarr('test.zarr')

ds3 = ds2.chunk({'bar': 2})

ds3.foo.encoding = {} ds3.coords['baz'].encoding = {} # Need this, too.

ds3.to_zarr('test3.zarr', mode='w') ```

Just flying by and dropping a note because I just ran into this with Imaris Open files being created by a microscope camera. I wanted to use one of my favorite packages (xarray) to dig into the data, and noted the dimension reuse. Not a big blocker, but this functionality of the data format might be growing in usage.

More details on Imaris. http://open.bitplane.com/Default.aspx?tabid=268