Repository: pydata/xarray

Issue:
Orthogonal indexing and MemoryCachedArray
While working on https://github.com/pydata/xarray/pull/1260 I came upon this which looks like a bug in caching:

```python
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)
```


Discussion:
fujiisoup: @fmaussion 
How did you notice issue #1429?
I found test_backends.py already checked the orthogonal indexing.

I guess there are some untested use cases.
shoyer: Also, I should say: thank you again for diving into the internals here. This sort of clean-up work is extremely important, even though it isn't as visible as new features.
fmaussion: >  @fmaussion How did you notice issue #1429?
> I found test_backends.py already checked the orthogonal indexing.
> I guess there are some untested use cases.

I don't remember exactly, but it happened when I wrote tests for rasterio's indexing: https://github.com/pydata/xarray/blob/master/xarray/tests/test_backends.py#L1855-L1884 (I just tested all possible combinations on index operations I could think of). The caching mechanism has some non-intuitive consequences, for example that some indexing operations not allowed on a rasterio file stored on disk become available once the data is cached:  https://github.com/pydata/xarray/blob/master/xarray/tests/test_backends.py#L1916-L1927
fujiisoup: @fmaussion 
Thanks for the information. Is there any need to add indexing tests for `rasterio` backend, after this issue #1429 is resolved? 

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 7b4eca02..99610640 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -25,10 +25,19 @@ Bug fixes
   of ``.T`` attributes (:issue:`1675`).
   By `Keisuke Fujii <https://github.com/fujiisoup>`_
 
+- (Internal bug) MemoryCachedArray now supports the orthogonal indexing.
+  Also made some internal cleanups around array wrappers (:issue:`1429`).
+
 - Fix two bugs that were preventing dask arrays from being specified as
   coordinates in the DataArray constructor (:issue:`1684`).
   By `Joe Hamman <https://github.com/jhamman>`_
 
+Testing
+~~~~~~~
+
+- Remove netCDF dependency from rasterio backend tests.
+  By `Matti Eskelinen <https://github.com/maaleske>`_
+
 v0.10.0 rc1 (30 October 2017)
 -----------------------------
 
diff --git a/xarray/backends/h5netcdf_.py b/xarray/backends/h5netcdf_.py
index 5bd4fa6c..edd5b34d 100644
--- a/xarray/backends/h5netcdf_.py
+++ b/xarray/backends/h5netcdf_.py
@@ -15,6 +15,10 @@ from .netCDF4_ import (_nc4_group, _nc4_values_and_dtype,
 
 class H5NetCDFArrayWrapper(BaseNetCDF4Array):
     def __getitem__(self, key):
+        if isinstance(key, indexing.VectorizedIndexer):
+            raise NotImplementedError(
+                 'Vectorized indexing for {} is not implemented. Load your '
+                 'data first with .load() or .compute().'.format(type(self)))
         key = indexing.to_tuple(key)
         with self.datastore.ensure_open(autoclose=True):
             return self.get_array()[key]
diff --git a/xarray/backends/netCDF4_.py b/xarray/backends/netCDF4_.py
index 52785f30..87bea8fe 100644
--- a/xarray/backends/netCDF4_.py
+++ b/xarray/backends/netCDF4_.py
@@ -27,7 +27,8 @@ _endian_lookup = {'=': 'native',
                   '|': 'native'}
 
 
-class BaseNetCDF4Array(NdimSizeLenMixin, DunderArrayMixin):
+class BaseNetCDF4Array(NdimSizeLenMixin, DunderArrayMixin,
+                       indexing.NDArrayIndexable):
     def __init__(self, variable_name, datastore):
         self.datastore = datastore
         self.variable_name = variable_name
@@ -50,6 +51,11 @@ class BaseNetCDF4Array(NdimSizeLenMixin, DunderArrayMixin):
 
 class NetCDF4ArrayWrapper(BaseNetCDF4Array):
     def __getitem__(self, key):
+        if isinstance(key, indexing.VectorizedIndexer):
+            raise NotImplementedError(
+             'Vectorized indexing for {} is not implemented. Load your '
+             'data first with .load() or .compute().'.format(type(self)))
+
         key = indexing.to_tuple(key)
 
         if self.datastore.is_remote:  # pragma: no cover
diff --git a/xarray/backends/pydap_.py b/xarray/backends/pydap_.py
index 826edd2d..2ab7873c 100644
--- a/xarray/backends/pydap_.py
+++ b/xarray/backends/pydap_.py
@@ -11,7 +11,7 @@ from ..core.pycompat import integer_types
 from .common import AbstractDataStore, robust_getitem
 
 
-class PydapArrayWrapper(NDArrayMixin):
+class PydapArrayWrapper(NDArrayMixin, indexing.NDArrayIndexable):
     def __init__(self, array):
         self.array = array
 
@@ -27,6 +27,10 @@ class PydapArrayWrapper(NDArrayMixin):
             return np.dtype(t.typecode + str(t.size))
 
     def __getitem__(self, key):
+        if isinstance(key, indexing.VectorizedIndexer):
+            raise NotImplementedError(
+             'Vectorized indexing for {} is not implemented. Load your '
+             'data first with .load() or .compute().'.format(type(self)))
         key = indexing.to_tuple(key)
         if not isinstance(key, tuple):
             key = (key,)
diff --git a/xarray/backends/pynio_.py b/xarray/backends/pynio_.py
index 98abb768..3e74e3bd 100644
--- a/xarray/backends/pynio_.py
+++ b/xarray/backends/pynio_.py
@@ -14,7 +14,8 @@ from ..core import indexing
 from .common import AbstractDataStore, DataStorePickleMixin
 
 
-class NioArrayWrapper(NdimSizeLenMixin, DunderArrayMixin):
+class NioArrayWrapper(NdimSizeLenMixin, DunderArrayMixin,
+                      indexing.NDArrayIndexable):
 
     def __init__(self, variable_name, datastore):
         self.datastore = datastore
@@ -28,6 +29,10 @@ class NioArrayWrapper(NdimSizeLenMixin, DunderArrayMixin):
         return self.datastore.ds.variables[self.variable_name]
 
     def __getitem__(self, key):
+        if isinstance(key, indexing.VectorizedIndexer):
+            raise NotImplementedError(
+             'Vectorized indexing for {} is not implemented. Load your '
+             'data first with .load() or .compute().'.format(type(self)))
         key = indexing.to_tuple(key)
         with self.datastore.ensure_open(autoclose=True):
             array = self.get_array()
diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py
index 938cea6b..b39299f3 100644
--- a/xarray/backends/rasterio_.py
+++ b/xarray/backends/rasterio_.py
@@ -1,6 +1,5 @@
 import os
 from collections import OrderedDict
-from distutils.version import LooseVersion
 import numpy as np
 
 from .. import DataArray
@@ -18,7 +17,8 @@ _ERROR_MSG = ('The kind of indexing operation you are trying to do is not '
               'first.')
 
 
-class RasterioArrayWrapper(NdimSizeLenMixin, DunderArrayMixin):
+class RasterioArrayWrapper(NdimSizeLenMixin, DunderArrayMixin,
+                           indexing.NDArrayIndexable):
     """A wrapper around rasterio dataset objects"""
     def __init__(self, rasterio_ds):
         self.rasterio_ds = rasterio_ds
@@ -38,8 +38,11 @@ class RasterioArrayWrapper(NdimSizeLenMixin, DunderArrayMixin):
         return self._shape
 
     def __getitem__(self, key):
+        if isinstance(key, indexing.VectorizedIndexer):
+            raise NotImplementedError(
+             'Vectorized indexing for {} is not implemented. Load your '
+             'data first with .load() or .compute().'.format(type(self)))
         key = indexing.to_tuple(key)
-
         # bands cannot be windowed but they can be listed
         band_key = key[0]
         n_bands = self.shape[0]
diff --git a/xarray/backends/scipy_.py b/xarray/backends/scipy_.py
index 798b41aa..346caf76 100644
--- a/xarray/backends/scipy_.py
+++ b/xarray/backends/scipy_.py
@@ -11,7 +11,7 @@ from .. import Variable
 from ..core.pycompat import iteritems, OrderedDict, basestring
 from ..core.utils import (Frozen, FrozenOrderedDict, NdimSizeLenMixin,
                           DunderArrayMixin)
-from ..core.indexing import NumpyIndexingAdapter
+from ..core.indexing import NumpyIndexingAdapter, NDArrayIndexable
 
 from .common import WritableCFDataStore, DataStorePickleMixin
 from .netcdf3 import (is_valid_nc3_name, encode_nc3_attr_value,
@@ -31,7 +31,7 @@ def _decode_attrs(d):
                        for (k, v) in iteritems(d))
 
 
-class ScipyArrayWrapper(NdimSizeLenMixin, DunderArrayMixin):
+class ScipyArrayWrapper(NdimSizeLenMixin, DunderArrayMixin, NDArrayIndexable):
 
     def __init__(self, variable_name, datastore):
         self.datastore = datastore
diff --git a/xarray/conventions.py b/xarray/conventions.py
index e78374b8..94c963da 100644
--- a/xarray/conventions.py
+++ b/xarray/conventions.py
@@ -333,7 +333,7 @@ def encode_cf_timedelta(timedeltas, units=None):
     return (num, units)
 
 
-class MaskedAndScaledArray(utils.NDArrayMixin):
+class MaskedAndScaledArray(utils.NDArrayMixin, indexing.NDArrayIndexable):
     """Wrapper around array-like objects to create a new indexable object where
     values, when accessed, are automatically scaled and masked according to
     CF conventions for packed and missing data values.
@@ -370,7 +370,7 @@ class MaskedAndScaledArray(utils.NDArrayMixin):
             After applying scale_factor, add this number to entries in the
             original array.
         """
-        self.array = array
+        self.array = indexing.as_indexable(array)
         self.fill_value = fill_value
         self.scale_factor = scale_factor
         self.add_offset = add_offset
@@ -391,13 +391,13 @@ class MaskedAndScaledArray(utils.NDArrayMixin):
                  self.scale_factor, self.add_offset, self._dtype))
 
 
-class DecodedCFDatetimeArray(utils.NDArrayMixin):
+class DecodedCFDatetimeArray(utils.NDArrayMixin, indexing.NDArrayIndexable):
     """Wrapper around array-like objects to create a new indexable object where
     values, when accessed, are automatically converted into datetime objects
     using decode_cf_datetime.
     """
     def __init__(self, array, units, calendar=None):
-        self.array = array
+        self.array = indexing.as_indexable(array)
         self.units = units
         self.calendar = calendar
 
@@ -430,13 +430,13 @@ class DecodedCFDatetimeArray(utils.NDArrayMixin):
                                   calendar=self.calendar)
 
 
-class DecodedCFTimedeltaArray(utils.NDArrayMixin):
+class DecodedCFTimedeltaArray(utils.NDArrayMixin, indexing.NDArrayIndexable):
     """Wrapper around array-like objects to create a new indexable object where
     values, when accessed, are automatically converted into timedelta objects
     using decode_cf_timedelta.
     """
     def __init__(self, array, units):
-        self.array = array
+        self.array = indexing.as_indexable(array)
         self.units = units
 
     @property
@@ -447,7 +447,7 @@ class DecodedCFTimedeltaArray(utils.NDArrayMixin):
         return decode_cf_timedelta(self.array[key], units=self.units)
 
 
-class StackedBytesArray(utils.NDArrayMixin):
+class StackedBytesArray(utils.NDArrayMixin, indexing.NDArrayIndexable):
     """Wrapper around array-like objects to create a new indexable object where
     values, when accessed, are automatically stacked along the last dimension.
 
@@ -465,7 +465,7 @@ class StackedBytesArray(utils.NDArrayMixin):
         if array.dtype != 'S1':
             raise ValueError(
                 "can only use StackedBytesArray if argument has dtype='S1'")
-        self.array = array
+        self.array = indexing.as_indexable(array)
 
     @property
     def dtype(self):
@@ -493,7 +493,7 @@ class StackedBytesArray(utils.NDArrayMixin):
         return char_to_bytes(self.array[key])
 
 
-class BytesToStringArray(utils.NDArrayMixin):
+class BytesToStringArray(utils.NDArrayMixin, indexing.NDArrayIndexable):
     """Wrapper that decodes bytes to unicode when values are read.
 
     >>> BytesToStringArray(np.array([b'abc']))[:]
@@ -509,7 +509,7 @@ class BytesToStringArray(utils.NDArrayMixin):
         encoding : str
             String encoding to use.
         """
-        self.array = array
+        self.array = indexing.as_indexable(array)
         self.encoding = encoding
 
     @property
@@ -532,7 +532,7 @@ class BytesToStringArray(utils.NDArrayMixin):
         return decode_bytes_array(self.array[key], self.encoding)
 
 
-class NativeEndiannessArray(utils.NDArrayMixin):
+class NativeEndiannessArray(utils.NDArrayMixin, indexing.NDArrayIndexable):
     """Decode arrays on the fly from non-native to native endianness
 
     This is useful for decoding arrays from netCDF3 files (which are all
@@ -551,7 +551,7 @@ class NativeEndiannessArray(utils.NDArrayMixin):
     dtype('int16')
     """
     def __init__(self, array):
-        self.array = array
+        self.array = indexing.as_indexable(array)
 
     @property
     def dtype(self):
@@ -561,7 +561,7 @@ class NativeEndiannessArray(utils.NDArrayMixin):
         return np.asarray(self.array[key], dtype=self.dtype)
 
 
-class BoolTypeArray(utils.NDArrayMixin):
+class BoolTypeArray(utils.NDArrayMixin, indexing.NDArrayIndexable):
     """Decode arrays on the fly from integer to boolean datatype
 
     This is useful for decoding boolean arrays from integer typed netCDF
@@ -579,7 +579,7 @@ class BoolTypeArray(utils.NDArrayMixin):
     dtype('bool')
     """
     def __init__(self, array):
-        self.array = array
+        self.array = indexing.as_indexable(array)
 
     @property
     def dtype(self):
@@ -589,7 +589,7 @@ class BoolTypeArray(utils.NDArrayMixin):
         return np.asarray(self.array[key], dtype=self.dtype)
 
 
-class UnsignedIntTypeArray(utils.NDArrayMixin):
+class UnsignedIntTypeArray(utils.NDArrayMixin, indexing.NDArrayIndexable):
     """Decode arrays on the fly from signed integer to unsigned
     integer. Typically used when _Unsigned is set at as a netCDF
     attribute on a signed integer variable.
@@ -606,7 +606,7 @@ class UnsignedIntTypeArray(utils.NDArrayMixin):
     array([  0,   1, 127, 128, 255], dtype=uint8)
     """
     def __init__(self, array):
-        self.array = array
+        self.array = indexing.as_indexable(array)
         self.unsigned_dtype = np.dtype('u%s' % array.dtype.itemsize)
 
     @property
diff --git a/xarray/core/indexing.py b/xarray/core/indexing.py
index 956e5e2d..6b03b64e 100644
--- a/xarray/core/indexing.py
+++ b/xarray/core/indexing.py
@@ -303,7 +303,11 @@ class VectorizedIndexer(IndexerTuple):
     """ Tuple for vectorized indexing """
 
 
-class LazilyIndexedArray(utils.NDArrayMixin):
+class NDArrayIndexable(object):
+    """ Mixin to mark support for IndexerTuple subclasses in indexing."""
+
+
+class LazilyIndexedArray(utils.NDArrayMixin, NDArrayIndexable):
     """Wrap an array that handles orthogonal indexing to make indexing lazy
     """
     def __init__(self, array, key=None):
@@ -356,7 +360,7 @@ class LazilyIndexedArray(utils.NDArrayMixin):
         return tuple(shape)
 
     def __array__(self, dtype=None):
-        array = xarray_indexable(self.array)
+        array = as_indexable(self.array)
         return np.asarray(array[self.key], dtype=None)
 
     def __getitem__(self, key):
@@ -379,9 +383,9 @@ def _wrap_numpy_scalars(array):
         return array
 
 
-class CopyOnWriteArray(utils.NDArrayMixin):
+class CopyOnWriteArray(utils.NDArrayMixin, NDArrayIndexable):
     def __init__(self, array):
-        self.array = array
+        self.array = as_indexable(array)
         self._copied = False
 
     def _ensure_copied(self):
@@ -400,9 +404,9 @@ class CopyOnWriteArray(utils.NDArrayMixin):
         self.array[key] = value
 
 
-class MemoryCachedArray(utils.NDArrayMixin):
+class MemoryCachedArray(utils.NDArrayMixin, NDArrayIndexable):
     def __init__(self, array):
-        self.array = _wrap_numpy_scalars(array)
+        self.array = _wrap_numpy_scalars(as_indexable(array))
 
     def _ensure_cached(self):
         if not isinstance(self.array, np.ndarray):
@@ -419,14 +423,21 @@ class MemoryCachedArray(utils.NDArrayMixin):
         self.array[key] = value
 
 
-def xarray_indexable(array):
+def as_indexable(array):
+    """
+    This function always returns a NDArrayIndexable subclass,
+    so that the vectorized indexing is always possible with the returned
+    object.
+    """
+    if isinstance(array, NDArrayIndexable):
+        return array
     if isinstance(array, np.ndarray):
         return NumpyIndexingAdapter(array)
     if isinstance(array, pd.Index):
         return PandasIndexAdapter(array)
     if isinstance(array, dask_array_type):
         return DaskIndexingAdapter(array)
-    return array
+    raise TypeError('Invalid array type: {}'.format(type(array)))
 
 
 def _outer_to_numpy_indexer(key, shape):
@@ -467,10 +478,14 @@ def _outer_to_numpy_indexer(key, shape):
     return tuple(new_key)
 
 
-class NumpyIndexingAdapter(utils.NDArrayMixin):
+class NumpyIndexingAdapter(utils.NDArrayMixin, NDArrayIndexable):
     """Wrap a NumPy array to use broadcasted indexing
     """
     def __init__(self, array):
+        # In NumpyIndexingAdapter we only allow to store bare np.ndarray
+        if not isinstance(array, np.ndarray):
+            raise TypeError('NumpyIndexingAdapter only wraps np.ndarray. '
+                            'Trying to wrap {}'.format(type(array)))
         self.array = array
 
     def _ensure_ndarray(self, value):
@@ -502,7 +517,7 @@ class NumpyIndexingAdapter(utils.NDArrayMixin):
         array[key] = value
 
 
-class DaskIndexingAdapter(utils.NDArrayMixin):
+class DaskIndexingAdapter(utils.NDArrayMixin, NDArrayIndexable):
     """Wrap a dask array to support xarray-style indexing.
     """
     def __init__(self, array):
@@ -522,6 +537,8 @@ class DaskIndexingAdapter(utils.NDArrayMixin):
             return self.array[to_int_tuple(key)]
         elif isinstance(key, VectorizedIndexer):
             return self.array.vindex[to_int_tuple(tuple(key))]
+        elif key is Ellipsis:
+            return self.array
         else:
             assert isinstance(key, OuterIndexer)
             key = to_int_tuple(tuple(key))
@@ -543,7 +560,7 @@ class DaskIndexingAdapter(utils.NDArrayMixin):
                         'method or accessing its .values attribute.')
 
 
-class PandasIndexAdapter(utils.NDArrayMixin):
+class PandasIndexAdapter(utils.NDArrayMixin, NDArrayIndexable):
     """Wrap a pandas.Index to be better about preserving dtypes and to handle
     indexing by length 1 tuples like numpy
     """
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index 66b09f0d..14e27708 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -19,7 +19,7 @@ from . import ops
 from . import utils
 from .pycompat import (basestring, OrderedDict, zip, integer_types,
                        dask_array_type)
-from .indexing import (PandasIndexAdapter, xarray_indexable, BasicIndexer,
+from .indexing import (PandasIndexAdapter, as_indexable, BasicIndexer,
                        OuterIndexer, VectorizedIndexer)
 from .utils import OrderedSet
 
@@ -31,6 +31,8 @@ except ImportError:
     pass
 
 
+NON_NUMPY_SUPPORTED_ARRAY_TYPES = (
+    indexing.NDArrayIndexable, pd.Index) + dask_array_type
 BASIC_INDEXING_TYPES = integer_types + (slice,)
 
 
@@ -159,12 +161,7 @@ def as_compatible_data(data, fastpath=False):
     if isinstance(data, Variable):
         return data.data
 
-    # add a custom fast-path for dask.array to avoid expensive checks for the
-    # dtype attribute
-    if isinstance(data, dask_array_type):
-        return data
-
-    if isinstance(data, pd.Index):
+    if isinstance(data, NON_NUMPY_SUPPORTED_ARRAY_TYPES):
         return _maybe_wrap_data(data)
 
     if isinstance(data, tuple):
@@ -177,14 +174,6 @@ def as_compatible_data(data, fastpath=False):
     if isinstance(data, timedelta):
         data = np.timedelta64(getattr(data, 'value', data), 'ns')
 
-    if (not hasattr(data, 'dtype') or not isinstance(data.dtype, np.dtype) or
-            not hasattr(data, 'shape') or
-            isinstance(data, (np.string_, np.unicode_,
-                              np.datetime64, np.timedelta64))):
-        # data must be ndarray-like
-        # don't allow non-numpy dtypes (e.g., categories)
-        data = np.asarray(data)
-
     # we don't want nested self-described arrays
     data = getattr(data, 'values', data)
 
@@ -197,6 +186,9 @@ def as_compatible_data(data, fastpath=False):
         else:
             data = np.asarray(data)
 
+    # validate whether the data is valid data types
+    data = np.asarray(data)
+
     if isinstance(data, np.ndarray):
         if data.dtype.kind == 'O':
             data = _possibly_convert_objects(data)
@@ -317,7 +309,7 @@ class Variable(common.AbstractArray, utils.NdimSizeLenMixin):
 
     @property
     def _indexable_data(self):
-        return xarray_indexable(self._data)
+        return as_indexable(self._data)
 
     def load(self, **kwargs):
         """Manually trigger loading of this variable's data from disk or a


Test patch:
diff --git a/xarray/tests/__init__.py b/xarray/tests/__init__.py
index 13cdd99a..d4d1b092 100644
--- a/xarray/tests/__init__.py
+++ b/xarray/tests/__init__.py
@@ -14,6 +14,7 @@ import pytest
 
 from xarray.core import utils
 from xarray.core.pycompat import PY3
+from xarray.core.indexing import NDArrayIndexable
 from xarray.testing import assert_equal, assert_identical, assert_allclose
 from xarray.plot.utils import import_seaborn
 
@@ -198,7 +199,7 @@ class UnexpectedDataAccess(Exception):
     pass
 
 
-class InaccessibleArray(utils.NDArrayMixin):
+class InaccessibleArray(utils.NDArrayMixin, NDArrayIndexable):
 
     def __init__(self, array):
         self.array = array
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index eb3176ed..f8a2ff9c 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -31,7 +31,7 @@ from . import (TestCase, requires_scipy, requires_netCDF4, requires_pydap,
                assert_identical, raises_regex)
 from .test_dataset import create_test_data
 
-from xarray.tests import mock
+from xarray.tests import mock, assert_identical
 
 try:
     import netCDF4 as nc4
@@ -393,11 +393,34 @@ class DatasetIOTestCases(object):
                         'dim3': np.arange(5)}
             expected = in_memory.isel(**indexers)
             actual = on_disk.isel(**indexers)
-            self.assertDatasetIdentical(expected, actual)
+            assert_identical(expected, actual)
             # do it twice, to make sure we're switched from orthogonal -> numpy
             # when we cached the values
             actual = on_disk.isel(**indexers)
-            self.assertDatasetIdentical(expected, actual)
+            assert_identical(expected, actual)
+
+    def _test_vectorized_indexing(self, vindex_support=True):
+        # Make sure vectorized_indexing works or at least raises
+        # NotImplementedError
+        in_memory = create_test_data()
+        with self.roundtrip(in_memory) as on_disk:
+            indexers = {'dim1': DataArray([0, 2, 0], dims='a'),
+                        'dim2': DataArray([0, 2, 3], dims='a')}
+            expected = in_memory.isel(**indexers)
+            if vindex_support:
+                actual = on_disk.isel(**indexers)
+                assert_identical(expected, actual)
+                # do it twice, to make sure we're switched from
+                # orthogonal -> numpy when we cached the values
+                actual = on_disk.isel(**indexers)
+                assert_identical(expected, actual)
+            else:
+                with raises_regex(NotImplementedError, 'Vectorized indexing '):
+                    actual = on_disk.isel(**indexers)
+
+    def test_vectorized_indexing(self):
+        # This test should be overwritten if vindex is supported
+        self._test_vectorized_indexing(vindex_support=False)
 
 
 class CFEncodedDataTest(DatasetIOTestCases):
@@ -623,6 +646,9 @@ class CFEncodedDataTest(DatasetIOTestCases):
             with self.open(tmp_file) as actual:
                 self.assertDatasetIdentical(data, actual)
 
+    def test_vectorized_indexing(self):
+        self._test_vectorized_indexing(vindex_support=False)
+
 
 _counter = itertools.count()
 
@@ -947,6 +973,9 @@ class NetCDF4ViaDaskDataTest(NetCDF4DataTest):
         # caching behavior differs for dask
         pass
 
+    def test_vectorized_indexing(self):
+        self._test_vectorized_indexing(vindex_support=True)
+
 
 class NetCDF4ViaDaskDataTestAutocloseTrue(NetCDF4ViaDaskDataTest):
     autoclose = True
@@ -1523,7 +1552,6 @@ class DaskTest(TestCase, DatasetIOTestCases):
         with raises_regex(TypeError, 'supports writing Dataset'):
             save_mfdataset([da], ['dataarray'])
 
-
     @requires_pathlib
     def test_save_mfdataset_pathlib_roundtrip(self):
         original = Dataset({'foo': ('x', np.random.randn(10))})
@@ -1608,6 +1636,9 @@ class DaskTest(TestCase, DatasetIOTestCases):
         self.assertTrue(computed._in_memory)
         self.assertDataArrayAllClose(actual, computed, decode_bytes=False)
 
+    def test_vectorized_indexing(self):
+        self._test_vectorized_indexing(vindex_support=True)
+
 
 class DaskTestAutocloseTrue(DaskTest):
     autoclose = True
@@ -1701,7 +1732,35 @@ class TestPyNioAutocloseTrue(TestPyNio):
 @requires_rasterio
 class TestRasterio(TestCase):
 
-    def test_serialization_utm(self):
+    @requires_scipy_or_netCDF4
+    def test_serialization(self):
+        import rasterio
+        from rasterio.transform import from_origin
+
+        # Create a geotiff file in utm proj
+        with create_tmp_file(suffix='.tif') as tmp_file:
+            # data
+            nx, ny, nz = 4, 3, 3
+            data = np.arange(nx*ny*nz,
+                             dtype=rasterio.float32).reshape(nz, ny, nx)
+            transform = from_origin(5000, 80000, 1000, 2000.)
+            with rasterio.open(
+                    tmp_file, 'w',
+                    driver='GTiff', height=ny, width=nx, count=nz,
+                    crs={'units': 'm', 'no_defs': True, 'ellps': 'WGS84',
+                         'proj': 'utm', 'zone': 18},
+                    transform=transform,
+                    dtype=rasterio.float32) as s:
+                s.write(data)
+
+            # Write it to a netcdf and read again (roundtrip)
+            with xr.open_rasterio(tmp_file) as rioda:
+                with create_tmp_file(suffix='.nc') as tmp_nc_file:
+                    rioda.to_netcdf(tmp_nc_file)
+                    with xr.open_dataarray(tmp_nc_file) as ncds:
+                        assert_identical(rioda, ncds)
+
+    def test_utm(self):
         import rasterio
         from rasterio.transform import from_origin
 
@@ -1740,13 +1799,7 @@ class TestRasterio(TestCase):
                 assert 'transform' in rioda.attrs
                 assert isinstance(rioda.attrs['transform'], tuple)
 
-                # Write it to a netcdf and read again (roundtrip)
-                with create_tmp_file(suffix='.nc') as tmp_nc_file:
-                    rioda.to_netcdf(tmp_nc_file)
-                    with xr.open_dataarray(tmp_nc_file) as ncds:
-                        assert_identical(rioda, ncds)
-
-    def test_serialization_platecarree(self):
+    def test_platecarree(self):
 
         import rasterio
         from rasterio.transform import from_origin
@@ -1784,12 +1837,6 @@ class TestRasterio(TestCase):
                 assert 'transform' in rioda.attrs
                 assert isinstance(rioda.attrs['transform'], tuple)
 
-                # Write it to a netcdf and read again (roundtrip)
-                with create_tmp_file(suffix='.nc') as tmp_nc_file:
-                    rioda.to_netcdf(tmp_nc_file)
-                    with xr.open_dataarray(tmp_nc_file) as ncds:
-                        assert_identical(rioda, ncds)
-
     def test_indexing(self):
 
         import rasterio
@@ -1851,7 +1898,6 @@ class TestRasterio(TestCase):
                     actual.isel(x=[4, 2]).values
                 with raises_regex(IndexError, err_msg):
                     actual.isel(x=slice(5, 2, -1)).values
-
                 # Integer indexing
                 ex = expected.isel(band=1)
                 ac = actual.isel(band=1)
diff --git a/xarray/tests/test_variable.py b/xarray/tests/test_variable.py
index 35c64040..7dd121d3 100644
--- a/xarray/tests/test_variable.py
+++ b/xarray/tests/test_variable.py
@@ -15,9 +15,12 @@ import pandas as pd
 from xarray import Variable, IndexVariable, Coordinate, Dataset
 from xarray.core import indexing
 from xarray.core.variable import as_variable, as_compatible_data
-from xarray.core.indexing import PandasIndexAdapter, LazilyIndexedArray
+from xarray.core.indexing import (PandasIndexAdapter, LazilyIndexedArray,
+                                  NumpyIndexingAdapter, CopyOnWriteArray,
+                                  MemoryCachedArray, DaskIndexingAdapter)
 from xarray.core.pycompat import PY3, OrderedDict
 from xarray.core.common import full_like, zeros_like, ones_like
+from xarray.core.utils import NDArrayMixin
 
 from . import TestCase, source_ndarray, requires_dask, raises_regex
 
@@ -1676,8 +1679,97 @@ class TestAsCompatibleData(TestCase):
         self.assertVariableIdentical(ones_like(orig, dtype=int),
                                      full_like(orig, 1, dtype=int))
 
+    def test_unsupported_type(self):
+        # Non indexable type
+        class CustomArray(NDArrayMixin):
+            def __init__(self, array):
+                self.array = array
+
+        class CustomIndexable(CustomArray, indexing.NDArrayIndexable):
+            pass
+
+        array = CustomArray(np.arange(3))
+        orig = Variable(dims=('x'), data=array, attrs={'foo': 'bar'})
+        assert isinstance(orig._data, np.ndarray)  # should not be CustomArray
+
+        array = CustomIndexable(np.arange(3))
+        orig = Variable(dims=('x'), data=array, attrs={'foo': 'bar'})
+        assert isinstance(orig._data, CustomIndexable)
+
 
 def test_raise_no_warning_for_nan_in_binary_ops():
     with pytest.warns(None) as record:
         Variable('x', [1, 2, np.NaN]) > 0
     assert len(record) == 0
+
+
+class TestBackendIndexing(TestCase):
+    """    Make sure all the array wrappers can be indexed. """
+    def setUp(self):
+        self.d = np.random.random((10, 3)).astype(np.float64)
+
+    def check_orthogonal_indexing(self, v):
+        assert np.allclose(v.isel(x=[8, 3], y=[2, 1]),
+                           self.d[[8, 3]][:, [2, 1]])
+
+    def check_vectorized_indexing(self, v):
+        ind_x = Variable('z', [0, 2])
+        ind_y = Variable('z', [2, 1])
+        assert np.allclose(v.isel(x=ind_x, y=ind_y), self.d[ind_x, ind_y])
+
+    def test_NumpyIndexingAdapter(self):
+        v = Variable(dims=('x', 'y'), data=NumpyIndexingAdapter(self.d))
+        self.check_orthogonal_indexing(v)
+        self.check_vectorized_indexing(v)
+        # could not doubly wrapping
+        with raises_regex(TypeError, 'NumpyIndexingAdapter only wraps '):
+            v = Variable(dims=('x', 'y'), data=NumpyIndexingAdapter(
+                                            NumpyIndexingAdapter(self.d)))
+
+    def test_LazilyIndexedArray(self):
+        v = Variable(dims=('x', 'y'), data=LazilyIndexedArray(self.d))
+        self.check_orthogonal_indexing(v)
+        with raises_regex(NotImplementedError, 'Vectorized indexing for '):
+            self.check_vectorized_indexing(v)
+        # doubly wrapping
+        v = Variable(dims=('x', 'y'),
+                     data=LazilyIndexedArray(LazilyIndexedArray(self.d)))
+        self.check_orthogonal_indexing(v)
+        # hierarchical wrapping
+        v = Variable(dims=('x', 'y'),
+                     data=LazilyIndexedArray(NumpyIndexingAdapter(self.d)))
+        self.check_orthogonal_indexing(v)
+
+    def test_CopyOnWriteArray(self):
+        v = Variable(dims=('x', 'y'), data=CopyOnWriteArray(self.d))
+        self.check_orthogonal_indexing(v)
+        self.check_vectorized_indexing(v)
+        # doubly wrapping
+        v = Variable(dims=('x', 'y'),
+                     data=CopyOnWriteArray(LazilyIndexedArray(self.d)))
+        self.check_orthogonal_indexing(v)
+        with raises_regex(NotImplementedError, 'Vectorized indexing for '):
+            self.check_vectorized_indexing(v)
+
+    def test_MemoryCachedArray(self):
+        v = Variable(dims=('x', 'y'), data=MemoryCachedArray(self.d))
+        self.check_orthogonal_indexing(v)
+        self.check_vectorized_indexing(v)
+        # doubly wrapping
+        v = Variable(dims=('x', 'y'),
+                     data=CopyOnWriteArray(MemoryCachedArray(self.d)))
+        self.check_orthogonal_indexing(v)
+        self.check_vectorized_indexing(v)
+
+    @requires_dask
+    def test_DaskIndexingAdapter(self):
+        import dask.array as da
+        da = da.asarray(self.d)
+        v = Variable(dims=('x', 'y'), data=DaskIndexingAdapter(da))
+        self.check_orthogonal_indexing(v)
+        self.check_vectorized_indexing(v)
+        # doubly wrapping
+        v = Variable(dims=('x', 'y'),
+                     data=CopyOnWriteArray(DaskIndexingAdapter(da)))
+        self.check_orthogonal_indexing(v)
+        self.check_vectorized_indexing(v)