Orthogonal indexing and MemoryCachedArray

[fmaussion]

While working on #1260 I came upon this which looks like a bug in caching:

import numpy as np
import xarray as xr
from xarray.core import indexing

nx, ny = 8, 10
data = np.arange(nx*ny).reshape(ny, nx)
cached = indexing.MemoryCachedArray(data)

data = xr.DataArray(data=data, dims=('y', 'x'))
cached = xr.DataArray(data=cached, dims=('y', 'x'))

a = data.isel(x=[2, 4], y=[3, 5])
b = cached.isel(x=[2, 4], y=[3, 5])

The last line raises:

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-13-45cd1493cf6b> in <module>()
     11 
     12 a = data.isel(x=[2, 4], y=[3, 5])
---> 13 b = cached.isel(x=[2, 4], y=[3, 5])

/home/mowglie/Documents/git/xarray/xarray/core/dataarray.py in isel(self, drop, **indexers)
    668         DataArray.sel
    669         """
--> 670         ds = self._to_temp_dataset().isel(drop=drop, **indexers)
    671         return self._from_temp_dataset(ds)
    672 

/home/mowglie/Documents/git/xarray/xarray/core/dataset.py in isel(self, drop, **indexers)
   1141         for name, var in iteritems(self._variables):
   1142             var_indexers = dict((k, v) for k, v in indexers if k in var.dims)
-> 1143             new_var = var.isel(**var_indexers)
   1144             if not (drop and name in var_indexers):
   1145                 variables[name] = new_var

/home/mowglie/Documents/git/xarray/xarray/core/variable.py in isel(self, **indexers)
    547             if dim in indexers:
    548                 key[i] = indexers[dim]
--> 549         return self[tuple(key)]
    550 
    551     def squeeze(self, dim=None):

/home/mowglie/Documents/git/xarray/xarray/core/variable.py in __getitem__(self, key)
    380         # orthogonal indexing should ensure the dimensionality is consistent
    381         if hasattr(values, 'ndim'):
--> 382             assert values.ndim == len(dims), (values.ndim, len(dims))
    383         else:
    384             assert len(dims) == 0, len(dims)

AssertionError: (1, 2)

[shoyer]

I think I have a fix for this on my laptop -- I'll push it to GitHub tonight or tomorrow

[shoyer]

I thought this would be fixed #1433, but this is actually a different issue.

We need to wrap the indexing result here for orthogonal indexing, but that probably shouldn't happen in the MemoryCachedArray layer. For SciPy, this happens inside ScipyArrayWrapper.__getitem__. For NumPy, this currently happens inside indexing.orthogonally_indexable.

Right now, xarray assumes that unknown array type that makes it into Variable._data supports orthogonal indexing. This is probably a mistake: we really should be more explicit and make fewer assumptions, perhaps by explicitly listing allowed array types or adding a registration system.

[fujiisoup]

I noticed that PR #1639 does not fully support indexing for MemoryCachedArray.

I don't think I understand how MemoryCachedArray is being used in xarray.
What object would be stored in this array wrapper?

[shoyer]

I don't think I understand how MemoryCachedArray is being used in xarray.
What object would be stored in this array wrapper?

When not using dask and cache=True in open_dataset, we use it to wrap arrays loaded from netCDF files.

Currently, it looks like the object loaded into it is usually a CopyOnWriteArray, which we use to ensure that modifying an xarray object never modifies the file on disk:

In [24]: ds2 = xr.open_dataset('foo.nc', engine='netcdf4')

In [25]: ds2.x.variable._data
Out[25]: MemoryCachedArray(array=CopyOnWriteArray(array=LazilyIndexedArray(array=<xarray.backends.netCDF4_.NetCDF4ArrayWrapper object at 0x11b3f25c0>, key=OuterIndexer(slice(None, None, None),))))

[fujiisoup]

Thanks, @shoyer.

These array wrappers should not store a bare np.array but current code is error-prone.
I will send a fix.