Repository: pydata/xarray

Issue:
NotImplementedError: Vectorized indexing for <class 'xarray.core.indexing.LazilyIndexedArray'> is not implemented.
I think this is a regression in the current 0.10.0rc1:

#### Code Sample

```python
import xarray as xr
ds = xr.open_dataset('cesm_data.nc', decode_cf=False)
ds.temp.isel(time=ds.time < 274383)  # throws an error
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-18-a5c4179cd02d> in <module>()
----> 1 ds.temp.isel(time=ds.time < 274383)

~/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/dataarray.py in isel(self, drop, **indexers)
    717         DataArray.sel
    718         """
--> 719         ds = self._to_temp_dataset().isel(drop=drop, **indexers)
    720         return self._from_temp_dataset(ds)
    721 

~/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/dataset.py in isel(self, drop, **indexers)
   1278         for name, var in iteritems(self._variables):
   1279             var_indexers = {k: v for k, v in indexers_list if k in var.dims}
-> 1280             new_var = var.isel(**var_indexers)
   1281             if not (drop and name in var_indexers):
   1282                 variables[name] = new_var

~/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/variable.py in isel(self, **indexers)
    771             if dim in indexers:
    772                 key[i] = indexers[dim]
--> 773         return self[tuple(key)]
    774 
    775     def squeeze(self, dim=None):

~/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/variable.py in __getitem__(self, key)
    595         """
    596         dims, index_tuple, new_order = self._broadcast_indexes(key)
--> 597         data = self._indexable_data[index_tuple]
    598         if new_order:
    599             data = np.moveaxis(data, range(len(new_order)), new_order)

~/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/indexing.py in __getitem__(self, key)
    414 
    415     def __getitem__(self, key):
--> 416         return type(self)(_wrap_numpy_scalars(self.array[key]))
    417 
    418     def __setitem__(self, key, value):

~/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/indexing.py in __getitem__(self, key)
    394 
    395     def __getitem__(self, key):
--> 396         return type(self)(_wrap_numpy_scalars(self.array[key]))
    397 
    398     def __setitem__(self, key, value):

~/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/indexing.py in __getitem__(self, key)
    361 
    362     def __getitem__(self, key):
--> 363         return type(self)(self.array, self._updated_key(key))
    364 
    365     def __setitem__(self, key, value):

~/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/indexing.py in _updated_key(self, new_key)
    336             raise NotImplementedError(
    337                 'Vectorized indexing for {} is not implemented. Load your '
--> 338                 'data first with .load() or .compute().'.format(type(self)))
    339         new_key = iter(expanded_indexer(new_key, self.ndim))
    340         key = []

NotImplementedError: Vectorized indexing for <class 'xarray.core.indexing.LazilyIndexedArray'> is not implemented. Load your data first with .load() or .compute().

```
Here is the file: 
[cesm_data.nc.zip](https://github.com/pydata/xarray/files/1441729/cesm_data.nc.zip)

#### Expected Output

This used to work in v0.9

#### Output of ``xr.show_versions()``

<details>
INSTALLED VERSIONS
------------------
commit: None
python: 3.5.2.final.0
python-bits: 64
OS: Linux
OS-release: 4.10.0-38-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: en_US.UTF-8

xarray: 0.10.0rc1
pandas: 0.21.0
numpy: 1.13.3
scipy: 1.0.0
netCDF4: 1.3.0
h5netcdf: None
Nio: None
bottleneck: 1.2.1
cyordereddict: None
dask: 0.15.4
matplotlib: 2.1.0
cartopy: 0.15.1
seaborn: 0.8.1
setuptools: 36.6.0
pip: 9.0.1
conda: None
pytest: 3.2.3
IPython: 6.2.1
sphinx: 1.6.5
</details>


Discussion:
fujiisoup: I thought issue #1688 was already covered in `test_backend.py` but was not.
If someone could suggest an appropriate test in `test_backend.py`, it would be super nice.
fmaussion: This code reproduces the error (one can probably clean it up a little):

```
import xarray as xr
import numpy as np
a = np.arange(1e5)
da.to_netcdf('test.nc')
t = xr.open_dataarray('test.nc')
t.isel(t=t.t > 12)
```
fujiisoup: Found another bug. Not ready yet.
fujiisoup: I found and fixed the following bugs. 

1. Fix for #1579 was not perfect (indexing `dask` array by `np.int` fails).
2. `LazilyIndexingArray` always used `OuterIndexer`. 
shoyer: @fujiisoup this needs a merge from master after your last PR.
fujiisoup: @shoyer 
Thanks. Merged from the master.

But I think I should also add another test for #1694 , which seems to be solved by this PR.
fujiisoup: Found some an unfixed bug in #1676, where in some cases `MemoryCachedArray`, `CopyOnWriteArray`, `LazilyIndexedArray` has a bare `np.ndarray`, rather than `NumpyIndexingArray`.

I think this is now fixed.

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 99610640..b3cd22e3 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -15,9 +15,16 @@ What's New
 
 .. _whats-new.0.10.0:
 
+
+v0.10.0 (unreleased)
+--------------------
+
 Changes since v0.10.0 rc1 (Unreleased)
 --------------------------------------
 
+- Experimental support for the Dask collection interface (:issue:`1674`).
+  By `Matthew Rocklin <https://github.com/mrocklin>`_.
+
 Bug fixes
 ~~~~~~~~~
 
@@ -25,8 +32,16 @@ Bug fixes
   of ``.T`` attributes (:issue:`1675`).
   By `Keisuke Fujii <https://github.com/fujiisoup>`_
 
+- Fix a bug in lazily-indexing netCDF array. (:issue:`1688`)
+  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`).
+  By `Keisuke Fujii <https://github.com/fujiisoup>`_
+
+- (Internal bug) MemoryCachedArray now always wraps ``np.ndarray`` by
+  ``NumpyIndexingAdapter``. (:issue:`1694`)
+  By `Keisuke Fujii <https://github.com/fujiisoup>`_
 
 - Fix two bugs that were preventing dask arrays from being specified as
   coordinates in the DataArray constructor (:issue:`1684`).
diff --git a/xarray/backends/pynio_.py b/xarray/backends/pynio_.py
index 3e74e3bd..8b9e4974 100644
--- a/xarray/backends/pynio_.py
+++ b/xarray/backends/pynio_.py
@@ -10,6 +10,7 @@ from .. import Variable
 from ..core.utils import (FrozenOrderedDict, Frozen,
                           NdimSizeLenMixin, DunderArrayMixin)
 from ..core import indexing
+from ..core.pycompat import integer_types
 
 from .common import AbstractDataStore, DataStorePickleMixin
 
@@ -29,10 +30,12 @@ class NioArrayWrapper(NdimSizeLenMixin, DunderArrayMixin,
         return self.datastore.ds.variables[self.variable_name]
 
     def __getitem__(self, key):
-        if isinstance(key, indexing.VectorizedIndexer):
+        if isinstance(key, (indexing.VectorizedIndexer,
+                            indexing.OuterIndexer)):
             raise NotImplementedError(
-             'Vectorized indexing for {} is not implemented. Load your '
-             'data first with .load() or .compute().'.format(type(self)))
+                'Nio backend does not support vectorized / outer indexing. '
+                'Load your data first with .load() or .compute(). '
+                'Given {}'.format(key))
         key = indexing.to_tuple(key)
         with self.datastore.ensure_open(autoclose=True):
             array = self.get_array()
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index e8330ef6..1dac7233 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -576,6 +576,35 @@ class DataArray(AbstractArray, BaseDataObject):
             dataset[self.name] = self.variable
             return dataset
 
+    def __dask_graph__(self):
+        return self._to_temp_dataset().__dask_graph__()
+
+    def __dask_keys__(self):
+        return self._to_temp_dataset().__dask_keys__()
+
+    @property
+    def __dask_optimize__(self):
+        return self._to_temp_dataset().__dask_optimize__
+
+    @property
+    def __dask_scheduler__(self):
+        return self._to_temp_dataset().__dask_optimize__
+
+    def __dask_postcompute__(self):
+        func, args = self._to_temp_dataset().__dask_postcompute__()
+        return self._dask_finalize, (func, args, self.name)
+
+    def __dask_postpersist__(self):
+        func, args = self._to_temp_dataset().__dask_postpersist__()
+        return self._dask_finalize, (func, args, self.name)
+
+    @staticmethod
+    def _dask_finalize(results, func, args, name):
+        ds = func(results, *args)
+        variable = ds._variables.pop(_THIS_ARRAY)
+        coords = ds._variables
+        return DataArray(variable, coords, name=name, fastpath=True)
+
     def load(self, **kwargs):
         """Manually trigger loading of this array's data from disk or a
         remote source into memory and return this array.
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index abe32055..56c9df0a 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -493,6 +493,77 @@ class Dataset(Mapping, ImplementsDatasetReduce, BaseDataObject,
 
         return self
 
+    def __dask_graph__(self):
+        graphs = {k: v.__dask_graph__() for k, v in self.variables.items()}
+        graphs = {k: v for k, v in graphs.items() if v is not None}
+        if not graphs:
+            return None
+        else:
+            from dask import sharedict
+            return sharedict.merge(*graphs.values())
+
+    def __dask_keys__(self):
+        import dask
+        return [v.__dask_keys__() for v in self.variables.values()
+                if dask.is_dask_collection(v)]
+
+    @property
+    def __dask_optimize__(self):
+        import dask.array as da
+        return da.Array.__dask_optimize__
+
+    @property
+    def __dask_scheduler__(self):
+        import dask.array as da
+        return da.Array.__dask_scheduler__
+
+    def __dask_postcompute__(self):
+        import dask
+        info = [(True, k, v.__dask_postcompute__())
+                if dask.is_dask_collection(v) else
+                (False, k, v) for k, v in self._variables.items()]
+        return self._dask_postcompute, (info, self._coord_names, self._dims,
+                                        self._attrs, self._file_obj,
+                                        self._encoding)
+
+    def __dask_postpersist__(self):
+        import dask
+        info = [(True, k, v.__dask_postpersist__())
+                if dask.is_dask_collection(v) else
+                (False, k, v) for k, v in self._variables.items()]
+        return self._dask_postpersist, (info, self._coord_names, self._dims,
+                                        self._attrs, self._file_obj,
+                                        self._encoding)
+
+    @staticmethod
+    def _dask_postcompute(results, info, *args):
+        variables = OrderedDict()
+        results2 = list(results[::-1])
+        for is_dask, k, v in info:
+            if is_dask:
+                func, args2 = v
+                r = results2.pop()
+                result = func(r, *args2)
+            else:
+                result = v
+            variables[k] = result
+
+        final = Dataset._construct_direct(variables, *args)
+        return final
+
+    @staticmethod
+    def _dask_postpersist(dsk, info, *args):
+        variables = OrderedDict()
+        for is_dask, k, v in info:
+            if is_dask:
+                func, args2 = v
+                result = func(dsk, *args2)
+            else:
+                result = v
+            variables[k] = result
+
+        return Dataset._construct_direct(variables, *args)
+
     def compute(self, **kwargs):
         """Manually trigger loading of this dataset's data from disk or a
         remote source into memory and return a new dataset. The original is
diff --git a/xarray/core/indexing.py b/xarray/core/indexing.py
index 6b03b64e..183f4e5e 100644
--- a/xarray/core/indexing.py
+++ b/xarray/core/indexing.py
@@ -322,7 +322,7 @@ class LazilyIndexedArray(utils.NDArrayMixin, NDArrayIndexable):
         """
         # We need to avoid doubly wrapping.
         if isinstance(array, type(self)):
-            self.array = array.array
+            self.array = as_indexable(array.array)
             self.key = array.key
             if key is not None:
                 self.key = self._updated_key(key)
@@ -330,8 +330,8 @@ class LazilyIndexedArray(utils.NDArrayMixin, NDArrayIndexable):
         else:
             if key is None:
                 key = (slice(None),) * array.ndim
-                key = OuterIndexer(key)
-            self.array = array
+                key = BasicIndexer(key)
+            self.array = as_indexable(array)
             self.key = key
 
     def _updated_key(self, new_key):
@@ -347,6 +347,8 @@ class LazilyIndexedArray(utils.NDArrayMixin, NDArrayIndexable):
                 key.append(k)
             else:
                 key.append(_index_indexer_1d(k, next(new_key), size))
+        if all(isinstance(k, integer_types + (slice, )) for k in key):
+            return BasicIndexer(key)
         return OuterIndexer(key)
 
     @property
@@ -390,7 +392,7 @@ class CopyOnWriteArray(utils.NDArrayMixin, NDArrayIndexable):
 
     def _ensure_copied(self):
         if not self._copied:
-            self.array = np.array(self.array)
+            self.array = as_indexable(np.array(self.array))
             self._copied = True
 
     def __array__(self, dtype=None):
@@ -409,8 +411,8 @@ class MemoryCachedArray(utils.NDArrayMixin, NDArrayIndexable):
         self.array = _wrap_numpy_scalars(as_indexable(array))
 
     def _ensure_cached(self):
-        if not isinstance(self.array, np.ndarray):
-            self.array = np.asarray(self.array)
+        if not isinstance(self.array, NumpyIndexingAdapter):
+            self.array = NumpyIndexingAdapter(np.asarray(self.array))
 
     def __array__(self, dtype=None):
         self._ensure_cached()
@@ -521,8 +523,8 @@ class DaskIndexingAdapter(utils.NDArrayMixin, NDArrayIndexable):
     """Wrap a dask array to support xarray-style indexing.
     """
     def __init__(self, array):
-        """ This adapter is usually called in Variable.__getitem__ with
-        array=Variable._broadcast_indexes
+        """ This adapter is created in Variable.__getitem__ in
+        Variable._broadcast_indexes.
         """
         self.array = array
 
@@ -531,7 +533,8 @@ class DaskIndexingAdapter(utils.NDArrayMixin, NDArrayIndexable):
             # workaround for uint64 indexer (GH:1406)
             # TODO remove here after next dask release (0.15.3)
             return tuple([k.astype(int) if isinstance(k, np.ndarray)
-                          else k for k in key])
+                          else int(k) if isinstance(k, np.integer) else k
+                          for k in key])
 
         if isinstance(key, BasicIndexer):
             return self.array[to_int_tuple(key)]
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index 14e27708..aa64dcab 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -290,7 +290,7 @@ class Variable(common.AbstractArray, utils.NdimSizeLenMixin):
     def _in_memory(self):
         return (isinstance(self._data, (np.ndarray, np.number, PandasIndexAdapter)) or
                 (isinstance(self._data, indexing.MemoryCachedArray) and
-                 isinstance(self._data.array, np.ndarray)))
+                 isinstance(self._data.array, indexing.NumpyIndexingAdapter)))
 
     @property
     def data(self):
@@ -355,6 +355,41 @@ class Variable(common.AbstractArray, utils.NdimSizeLenMixin):
         new = self.copy(deep=False)
         return new.load(**kwargs)
 
+    def __dask_graph__(self):
+        if isinstance(self._data, dask_array_type):
+            return self._data.__dask_graph__()
+        else:
+            return None
+
+    def __dask_keys__(self):
+        return self._data.__dask_keys__()
+
+    @property
+    def __dask_optimize__(self):
+        return self._data.__dask_optimize__
+
+    @property
+    def __dask_scheduler__(self):
+        return self._data.__dask_scheduler__
+
+    def __dask_postcompute__(self):
+        array_func, array_args = self._data.__dask_postcompute__()
+        return self._dask_finalize, (array_func, array_args, self._dims,
+                                     self._attrs, self._encoding)
+
+    def __dask_postpersist__(self):
+        array_func, array_args = self._data.__dask_postpersist__()
+        return self._dask_finalize, (array_func, array_args, self._dims,
+                                     self._attrs, self._encoding)
+
+    @staticmethod
+    def _dask_finalize(results, array_func, array_args, dims, attrs, encoding):
+        if isinstance(results, dict):  # persist case
+            name = array_args[0]
+            results = {k: v for k, v in results.items() if k[0] == name}  # cull
+        data = array_func(results, *array_args)
+        return Variable(dims, data, attrs=attrs, encoding=encoding)
+
     @property
     def values(self):
         """The variable's data as a numpy.ndarray"""
@@ -453,9 +488,8 @@ class Variable(common.AbstractArray, utils.NdimSizeLenMixin):
                 if len(k.dims) > 1:
                     return self._broadcast_indexes_vectorized(key)
                 dims.append(k.dims[0])
-            if not isinstance(k, integer_types):
+            elif not isinstance(k, integer_types):
                 dims.append(d)
-
         if len(set(dims)) == len(dims):
             return self._broadcast_indexes_outer(key)
 


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index f8a2ff9c..1701b8d4 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -22,7 +22,8 @@ from xarray import (Dataset, DataArray, open_dataset, open_dataarray,
 from xarray.backends.common import robust_getitem
 from xarray.backends.netCDF4_ import _extract_nc4_variable_encoding
 from xarray.core import indexing
-from xarray.core.pycompat import iteritems, PY2, ExitStack, basestring
+from xarray.core.pycompat import (iteritems, PY2, ExitStack, basestring,
+                                  dask_array_type)
 
 from . import (TestCase, requires_scipy, requires_netCDF4, requires_pydap,
                requires_scipy_or_netCDF4, requires_dask, requires_h5netcdf,
@@ -422,6 +423,51 @@ class DatasetIOTestCases(object):
         # This test should be overwritten if vindex is supported
         self._test_vectorized_indexing(vindex_support=False)
 
+    def test_isel_dataarray(self):
+        # Make sure isel works lazily. GH:issue:1688
+        in_memory = create_test_data()
+        with self.roundtrip(in_memory) as on_disk:
+            expected = in_memory.isel(dim2=in_memory['dim2'] < 3)
+            actual = on_disk.isel(dim2=on_disk['dim2'] < 3)
+            self.assertDatasetIdentical(expected, actual)
+
+    def validate_array_type(self, ds):
+        # Make sure that only NumpyIndexingAdapter stores a bare np.ndarray.
+        def find_and_validate_array(obj):
+            # recursively called function. obj: array or array wrapper.
+            if hasattr(obj, 'array'):
+                if isinstance(obj.array, indexing.NDArrayIndexable):
+                    find_and_validate_array(obj.array)
+                else:
+                    if isinstance(obj.array, np.ndarray):
+                        assert isinstance(obj, indexing.NumpyIndexingAdapter)
+                    elif isinstance(obj.array, dask_array_type):
+                        assert isinstance(obj, indexing.DaskIndexingAdapter)
+                    elif isinstance(obj.array, pd.Index):
+                        assert isinstance(obj, indexing.PandasIndexAdapter)
+                    else:
+                        raise TypeError('{} is wrapped by {}'.format(
+                                            type(obj.array), type(obj)))
+
+        for k, v in ds.variables.items():
+            find_and_validate_array(v._data)
+
+    def test_array_type_after_indexing(self):
+        in_memory = create_test_data()
+        with self.roundtrip(in_memory) as on_disk:
+            self.validate_array_type(on_disk)
+            indexers = {'dim1': [1, 2, 0], 'dim2': [3, 2, 0, 3],
+                        'dim3': np.arange(5)}
+            expected = in_memory.isel(**indexers)
+            actual = on_disk.isel(**indexers)
+            assert_identical(expected, actual)
+            self.validate_array_type(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)
+            self.validate_array_type(actual)
+
 
 class CFEncodedDataTest(DatasetIOTestCases):
 
@@ -1192,6 +1238,10 @@ class H5NetCDFDataTest(BaseNetCDF4Test, TestCase):
         # doesn't work for h5py (without using dask as an intermediate layer)
         pass
 
+    def test_array_type_after_indexing(self):
+        # pynio also does not support list-like indexing
+        pass
+
     def test_complex(self):
         expected = Dataset({'x': ('y', np.ones(5) + 1j * np.ones(5))})
         with self.roundtrip(expected) as actual:
@@ -1701,6 +1751,15 @@ class TestPyNio(CFEncodedDataTest, NetCDF3Only, TestCase):
         pass
 
     def test_orthogonal_indexing(self):
+        # pynio also does not support list-like indexing
+        with raises_regex(NotImplementedError, 'Nio backend does not '):
+            super(TestPyNio, self).test_orthogonal_indexing()
+
+    def test_isel_dataarray(self):
+        with raises_regex(NotImplementedError, 'Nio backend does not '):
+            super(TestPyNio, self).test_isel_dataarray()
+
+    def test_array_type_after_indexing(self):
         # pynio also does not support list-like indexing
         pass
 
diff --git a/xarray/tests/test_dask.py b/xarray/tests/test_dask.py
index 6122aaf4..4e9b0250 100644
--- a/xarray/tests/test_dask.py
+++ b/xarray/tests/test_dask.py
@@ -206,6 +206,34 @@ class TestVariable(DaskTestCase):
         self.assertLazyAndAllClose(np.maximum(u, 0), xu.maximum(v, 0))
         self.assertLazyAndAllClose(np.maximum(u, 0), xu.maximum(0, v))
 
+    @pytest.mark.skipif(LooseVersion(dask.__version__) <= '0.15.4',
+                        reason='Need dask 0.16 for new interface')
+    def test_compute(self):
+        u = self.eager_var
+        v = self.lazy_var
+
+        assert dask.is_dask_collection(v)
+        (v2,) = dask.compute(v + 1)
+        assert not dask.is_dask_collection(v2)
+
+        assert ((u + 1).data == v2.data).all()
+
+    @pytest.mark.skipif(LooseVersion(dask.__version__) <= '0.15.4',
+                        reason='Need dask 0.16 for new interface')
+    def test_persist(self):
+        u = self.eager_var
+        v = self.lazy_var + 1
+
+        (v2,) = dask.persist(v)
+        assert v is not v2
+        assert len(v2.__dask_graph__()) < len(v.__dask_graph__())
+        assert v2.__dask_keys__() == v.__dask_keys__()
+        assert dask.is_dask_collection(v)
+        assert dask.is_dask_collection(v2)
+
+        self.assertLazyAndAllClose(u + 1, v)
+        self.assertLazyAndAllClose(u + 1, v2)
+
 
 class TestDataArrayAndDataset(DaskTestCase):
     def assertLazyAndIdentical(self, expected, actual):
@@ -251,6 +279,34 @@ class TestDataArrayAndDataset(DaskTestCase):
         actual = xr.concat([v[:2], v[2:]], 'x')
         self.assertLazyAndAllClose(u, actual)
 
+    @pytest.mark.skipif(LooseVersion(dask.__version__) <= '0.15.4',
+                        reason='Need dask 0.16 for new interface')
+    def test_compute(self):
+        u = self.eager_array
+        v = self.lazy_array
+
+        assert dask.is_dask_collection(v)
+        (v2,) = dask.compute(v + 1)
+        assert not dask.is_dask_collection(v2)
+
+        assert ((u + 1).data == v2.data).all()
+
+    @pytest.mark.skipif(LooseVersion(dask.__version__) <= '0.15.4',
+                        reason='Need dask 0.16 for new interface')
+    def test_persist(self):
+        u = self.eager_array
+        v = self.lazy_array + 1
+
+        (v2,) = dask.persist(v)
+        assert v is not v2
+        assert len(v2.__dask_graph__()) < len(v.__dask_graph__())
+        assert v2.__dask_keys__() == v.__dask_keys__()
+        assert dask.is_dask_collection(v)
+        assert dask.is_dask_collection(v2)
+
+        self.assertLazyAndAllClose(u + 1, v)
+        self.assertLazyAndAllClose(u + 1, v2)
+
     def test_concat_loads_variables(self):
         # Test that concat() computes not-in-memory variables at most once
         # and loads them in the output, while leaving the input unaltered.
@@ -402,28 +458,6 @@ class TestDataArrayAndDataset(DaskTestCase):
             ds.load()
         self.assertEqual(count[0], 1)
 
-    def test_persist_Dataset(self):
-        ds = Dataset({'foo': ('x', range(5)),
-                      'bar': ('x', range(5))}).chunk()
-        ds = ds + 1
-        n = len(ds.foo.data.dask)
-
-        ds2 = ds.persist()
-
-        assert len(ds2.foo.data.dask) == 1
-        assert len(ds.foo.data.dask) == n  # doesn't mutate in place
-
-    def test_persist_DataArray(self):
-        x = da.arange(10, chunks=(5,))
-        y = DataArray(x)
-        z = y + 1
-        n = len(z.data.dask)
-
-        zz = z.persist()
-
-        assert len(z.data.dask) == n
-        assert len(zz.data.dask) == zz.data.npartitions
-
     def test_stack(self):
         data = da.random.normal(size=(2, 3, 4), chunks=(1, 3, 4))
         arr = DataArray(data, dims=('w', 'x', 'y'))
@@ -737,3 +771,62 @@ def build_dask_array(name):
     return dask.array.Array(
         dask={(name, 0): (kernel, name)}, name=name,
         chunks=((1,),), dtype=np.int64)
+
+
+# test both the perist method and the dask.persist function
+# the dask.persist function requires a new version of dask
+@pytest.mark.parametrize('persist', [
+    lambda x: x.persist(),
+    pytest.mark.skipif(LooseVersion(dask.__version__) <= '0.15.4',
+                       lambda x: dask.persist(x)[0],
+                       reason='Need Dask 0.16')
+])
+def test_persist_Dataset(persist):
+    ds = Dataset({'foo': ('x', range(5)),
+                  'bar': ('x', range(5))}).chunk()
+    ds = ds + 1
+    n = len(ds.foo.data.dask)
+
+    ds2 = persist(ds)
+
+    assert len(ds2.foo.data.dask) == 1
+    assert len(ds.foo.data.dask) == n  # doesn't mutate in place
+
+
+@pytest.mark.parametrize('persist', [
+    lambda x: x.persist(),
+    pytest.mark.skipif(LooseVersion(dask.__version__) <= '0.15.4',
+                       lambda x: dask.persist(x)[0],
+                       reason='Need Dask 0.16')
+])
+def test_persist_DataArray(persist):
+    x = da.arange(10, chunks=(5,))
+    y = DataArray(x)
+    z = y + 1
+    n = len(z.data.dask)
+
+    zz = persist(z)
+
+    assert len(z.data.dask) == n
+    assert len(zz.data.dask) == zz.data.npartitions
+
+
+@pytest.mark.skipif(LooseVersion(dask.__version__) <= '0.15.4',
+                    reason='Need dask 0.16 for new interface')
+def test_dataarray_with_dask_coords():
+    import toolz
+    x = xr.Variable('x', da.arange(8, chunks=(4,)))
+    y = xr.Variable('y', da.arange(8, chunks=(4,)) * 2)
+    data = da.random.random((8, 8), chunks=(4, 4)) + 1
+    array = xr.DataArray(data, dims=['x', 'y'])
+    array.coords['xx'] = x
+    array.coords['yy'] = y
+
+    assert dict(array.__dask_graph__()) == toolz.merge(data.__dask_graph__(),
+                                                       x.__dask_graph__(),
+                                                       y.__dask_graph__())
+
+    (array2,) = dask.compute(array)
+    assert not dask.is_dask_collection(array2)
+
+    assert all(isinstance(v._variable.data, np.ndarray) for v in array2.coords.values())
diff --git a/xarray/tests/test_distributed.py b/xarray/tests/test_distributed.py
index 1868486b..9999ed9a 100644
--- a/xarray/tests/test_distributed.py
+++ b/xarray/tests/test_distributed.py
@@ -4,7 +4,9 @@ from xarray.core.pycompat import suppress
 
 distributed = pytest.importorskip('distributed')
 da = pytest.importorskip('dask.array')
-from distributed.utils_test import cluster, loop
+import dask
+from distributed.utils_test import cluster, loop, gen_cluster
+from distributed.client import futures_of, wait
 
 from xarray.tests.test_backends import create_tmp_file, ON_WINDOWS
 from xarray.tests.test_dataset import create_test_data
@@ -32,3 +34,33 @@ def test_dask_distributed_integration_test(loop, engine):
                     assert isinstance(restored.var1.data, da.Array)
                     computed = restored.compute()
                     assert_allclose(original, computed)
+
+
+@pytest.mark.skipif(distributed.__version__ <= '1.19.3',
+                    reason='Need recent distributed version to clean up get')
+@gen_cluster(client=True, timeout=None)
+def test_async(c, s, a, b):
+    x = create_test_data()
+    assert not dask.is_dask_collection(x)
+    y = x.chunk({'dim2': 4}) + 10
+    assert dask.is_dask_collection(y)
+    assert dask.is_dask_collection(y.var1)
+    assert dask.is_dask_collection(y.var2)
+
+    z = y.persist()
+    assert str(z)
+
+    assert dask.is_dask_collection(z)
+    assert dask.is_dask_collection(z.var1)
+    assert dask.is_dask_collection(z.var2)
+    assert len(y.__dask_graph__()) > len(z.__dask_graph__())
+
+    assert not futures_of(y)
+    assert futures_of(z)
+
+    future = c.compute(z)
+    w = yield future
+    assert not dask.is_dask_collection(w)
+    assert_allclose(x + 10, w)
+
+    assert s.task_state
diff --git a/xarray/tests/test_indexing.py b/xarray/tests/test_indexing.py
index 4b78efb2..553c3012 100644
--- a/xarray/tests/test_indexing.py
+++ b/xarray/tests/test_indexing.py
@@ -11,6 +11,7 @@ import pandas as pd
 from xarray import Dataset, DataArray, Variable
 from xarray.core import indexing
 from xarray.core import nputils
+from xarray.core.pycompat import integer_types
 from . import TestCase, ReturnItem, raises_regex
 
 
@@ -169,6 +170,16 @@ class TestLazyArray(TestCase):
                         self.assertArrayEqual(expected, actual)
                         assert isinstance(actual._data,
                                           indexing.LazilyIndexedArray)
+
+                        # make sure actual.key is appropriate type
+                        if all(isinstance(k, integer_types + (slice, ))
+                               for k in v_lazy._data.key):
+                            assert isinstance(v_lazy._data.key,
+                                              indexing.BasicIndexer)
+                        else:
+                            assert isinstance(v_lazy._data.key,
+                                              indexing.OuterIndexer)
+
         # test sequentially applied indexers
         indexers = [(3, 2), (I[:], 0), (I[:2], -1), (I[:4], [0]), ([4, 5], 0),
                     ([0, 1, 2], [0, 1]), ([0, 3, 5], I[:2])]
@@ -211,7 +222,7 @@ class TestMemoryCachedArray(TestCase):
         original = indexing.LazilyIndexedArray(np.arange(10))
         wrapped = indexing.MemoryCachedArray(original)
         self.assertArrayEqual(wrapped, np.arange(10))
-        self.assertIsInstance(wrapped.array, np.ndarray)
+        self.assertIsInstance(wrapped.array, indexing.NumpyIndexingAdapter)
 
     def test_sub_array(self):
         original = indexing.LazilyIndexedArray(np.arange(10))
@@ -219,7 +230,7 @@ class TestMemoryCachedArray(TestCase):
         child = wrapped[:5]
         self.assertIsInstance(child, indexing.MemoryCachedArray)
         self.assertArrayEqual(child, np.arange(5))
-        self.assertIsInstance(child.array, np.ndarray)
+        self.assertIsInstance(child.array, indexing.NumpyIndexingAdapter)
         self.assertIsInstance(wrapped.array, indexing.LazilyIndexedArray)
 
     def test_setitem(self):
diff --git a/xarray/tests/test_variable.py b/xarray/tests/test_variable.py
index 7dd121d3..e07f01af 100644
--- a/xarray/tests/test_variable.py
+++ b/xarray/tests/test_variable.py
@@ -16,8 +16,10 @@ 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,
-                                  NumpyIndexingAdapter, CopyOnWriteArray,
-                                  MemoryCachedArray, DaskIndexingAdapter)
+                                  BasicIndexer, OuterIndexer,
+                                  VectorizedIndexer, 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
@@ -578,6 +580,9 @@ class VariableSubclassTestCases(object):
         v_new = v[np.array([0], dtype="uint64")]
         self.assertArrayEqual(v_new, v_data[[0], :])
 
+        v_new = v[np.uint64(0)]
+        self.assertArrayEqual(v_new, v_data[0, :])
+
     def test_getitem_0d_array(self):
         # make sure 0d-np.array can be used as an indexer
         v = self.cls(['x'], [0, 1, 2])
@@ -941,6 +946,39 @@ class TestVariable(TestCase, VariableSubclassTestCases):
         _, ind, _ = v._broadcast_indexes((vind, 3))
         assert type(ind) == indexing.VectorizedIndexer
 
+    def test_indexer_type(self):
+        # GH:issue:1688. Wrong indexer type induces NotImplementedError
+        data = np.random.random((10, 11))
+        v = Variable(['x', 'y'], data)
+
+        def assert_indexer_type(key, object_type):
+            dims, index_tuple, new_order = v._broadcast_indexes(key)
+            assert isinstance(index_tuple, object_type)
+
+        # should return BasicIndexer
+        assert_indexer_type((0, 1), BasicIndexer)
+        assert_indexer_type((0, slice(None, None)), BasicIndexer)
+        assert_indexer_type((Variable([], 3), slice(None, None)), BasicIndexer)
+        assert_indexer_type((Variable([], 3), (Variable([], 6))), BasicIndexer)
+
+        # should return OuterIndexer
+        assert_indexer_type(([0, 1], 1), OuterIndexer)
+        assert_indexer_type(([0, 1], [1, 2]), OuterIndexer)
+        assert_indexer_type((Variable(('x'), [0, 1]), 1), OuterIndexer)
+        assert_indexer_type((Variable(('x'), [0, 1]), slice(None, None)),
+                            OuterIndexer)
+        assert_indexer_type((Variable(('x'), [0, 1]), Variable(('y'), [0, 1])),
+                            OuterIndexer)
+
+        # should return VectorizedIndexer
+        assert_indexer_type((Variable(('y'), [0, 1]), [0, 1]),
+                            VectorizedIndexer)
+        assert_indexer_type((Variable(('z'), [0, 1]), Variable(('z'), [0, 1])),
+                            VectorizedIndexer)
+        assert_indexer_type((Variable(('a', 'b'), [[0, 1], [1, 2]]),
+                             Variable(('a', 'b'), [[0, 1], [1, 2]])),
+                            VectorizedIndexer)
+
     def test_items(self):
         data = np.random.random((10, 11))
         v = Variable(['x', 'y'], data)