Repository: pydata/xarray

Issue:
Regression: "ValueError: cannot unstack dimensions that do not have a MultiIndex" when unstacking a MultiIndex
<!-- Please include a self-contained copy-pastable example that generates the issue if possible.

Please be concise with code posted. See guidelines below on how to provide a good bug report:

- Craft Minimal Bug Reports: http://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports
- Minimal Complete Verifiable Examples: https://stackoverflow.com/help/mcve

Bug reports that follow these guidelines are easier to diagnose, and so are often handled much more quickly.
-->

I'm not sure if this is a bug or I'm not using `xarray` correctly, but I used to be able to do this without crashing. The new behavior seems to have been introduced some time between 0.16.2 and 0.18.2.

**What happened**:

```python traceback
Traceback (most recent call last):
  File "scripts/repro.py", line 12, in <module>
    ds = ds.unstack(['c'])
  File "/home/darsh/src/notebooks/build/venv/lib/python3.8/site-packages/xarray/core/dataset.py", line 4024, in unstack
    raise ValueError(
ValueError: cannot unstack dimensions that do not have a MultiIndex: ['c']
```

**What you expected to happen**:

The code runs without the `ValueError` exception.

**Minimal Complete Verifiable Example**:

```python
from xarray import DataArray, Dataset


a = DataArray([0], dims=['a'])
b = a.stack(b=('a',)).reset_index('b')
c = b.stack({'c': ['b']})

ds = Dataset({'d': DataArray(c.data, dims=['c'])}, coords=c.coords)
print('\nBefore:')
print(ds)

ds = ds.unstack(['c'])
print('\nAfter:')
print(ds)
```

**Anything else we need to know?**:

Here's the full output from the example on 0.18.2:
```

Before:
<xarray.Dataset>
Dimensions:  (c: 1)
Coordinates:
  * c        (c) MultiIndex
  - b        (c) int64 0
    a        (c) int64 0
Data variables:
    d        (c) int64 0
Traceback (most recent call last):
  File "scripts/repro.py", line 12, in <module>
    ds = ds.unstack(['c'])
  File "/home/darsh/src/notebooks/build/venv/lib/python3.8/site-packages/xarray/core/dataset.py", line 4024, in unstack
    raise ValueError(
ValueError: cannot unstack dimensions that do not have a MultiIndex: ['c']
```
What confuses me is that the `c` dimension is shown as a `MultiIndex`, but it still complains that it doesn't have a `MultiIndex`. Directly unstacking `ds.d` rather than the dataset itself also fails with the same exception.

Oddly, it seems to work if I assign the coordinates after constructing the dataset:
```diff
diff --git a/scripts/repro.py b/scripts/repro.py
index ed2ae7c..d5bd6a3 100644
--- a/scripts/repro.py
+++ b/scripts/repro.py
@@ -5,7 +5,7 @@ a = DataArray([0], dims=['a'])
 b = a.stack(b=('a',)).reset_index('b')
 c = b.stack({'c': ['b']})
 
-ds = Dataset({'d': DataArray(c.data, dims=['c'])}, coords=c.coords)
+ds = Dataset({'d': DataArray(c.data, dims=['c'])}).assign_coords(c.coords)
 print('\nBefore:')
 print(ds)
 

```

With that workaround, or by downgrading to 0.16.2, the example doesn't crash:

```

Before:
<xarray.Dataset>
Dimensions:  (c: 1)
Coordinates:
  * c        (c) MultiIndex
  - b        (c) int64 0
    a        (c) int64 0
Data variables:
    d        (c) int64 0

After:
<xarray.Dataset>
Dimensions:  (b: 1)
Coordinates:
    a        (b) int64 0
  * b        (b) int64 0
Data variables:
    d        (b) int64 0
```


**Environment**:

<details><summary>Output of <tt>xr.show_versions()</tt></summary>

<!-- Paste the output here xr.show_versions() here -->

INSTALLED VERSIONS
------------------
commit: None
python: 3.8.0 (default, Feb 25 2021, 22:10:10) 
[GCC 8.4.0]
python-bits: 64
OS: Linux
OS-release: 5.4.0-73-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: None
libnetcdf: None

xarray: 0.18.2
pandas: 1.2.4
numpy: 1.20.3
scipy: 1.6.3
netCDF4: None
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: None
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: None
dask: 2021.05.0
distributed: None
matplotlib: 3.4.2
cartopy: None
seaborn: None
numbagg: None
pint: 0.17
setuptools: 39.0.1
pip: 21.1.1
conda: None
pytest: 6.2.4
IPython: 7.23.1
sphinx: None
None


</details>


Discussion:

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 9ba9a7e9..fb96ed62 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -21,11 +21,14 @@ v0.18.3 (unreleased)
 
 New Features
 ~~~~~~~~~~~~
-
 - Allow assigning values to a subset of a dataset using positional or label-based
   indexing (:issue:`3015`, :pull:`5362`). By `Matthias Göbel <https://github.com/matzegoebel>`_.
 - Attempting to reduce a weighted object over missing dimensions now raises an error (:pull:`5362`).
   By `Mattia Almansi <https://github.com/malmans2>`_.
+- :py:func:`xarray.cov` and :py:func:`xarray.corr` now lazily check for missing
+  values if inputs are dask arrays (:issue:`4804`, :pull:`5284`).
+  By `Andrew Williams <https://github.com/AndrewWilliams3142>`_.
+
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
@@ -41,6 +44,9 @@ Bug fixes
   :py:class:`CFTimeIndex` and upcoming pandas version 1.3.0 (:issue:`5356`,
   :pull:`5359`).
   By `Spencer Clark <https://github.com/spencerkclark>`_.
+- Fix 1-level multi-index incorrectly converted to single index (:issue:`5384`,
+  :pull:`5385`).
+  By `Benoit Bovy <https://github.com/benbovy>`_.
 
 
 Documentation
diff --git a/xarray/core/computation.py b/xarray/core/computation.py
index 6010d502..07b7bca2 100644
--- a/xarray/core/computation.py
+++ b/xarray/core/computation.py
@@ -1352,12 +1352,23 @@ def _cov_corr(da_a, da_b, dim=None, ddof=0, method=None):
 
     # 2. Ignore the nans
     valid_values = da_a.notnull() & da_b.notnull()
+    valid_count = valid_values.sum(dim) - ddof
 
-    if not valid_values.all():
-        da_a = da_a.where(valid_values)
-        da_b = da_b.where(valid_values)
+    def _get_valid_values(da, other):
+        """
+        Function to lazily mask da_a and da_b
+        following a similar approach to
+        https://github.com/pydata/xarray/pull/4559
+        """
+        missing_vals = np.logical_or(da.isnull(), other.isnull())
+        if missing_vals.any():
+            da = da.where(~missing_vals)
+            return da
+        else:
+            return da
 
-    valid_count = valid_values.sum(dim) - ddof
+    da_a = da_a.map_blocks(_get_valid_values, args=[da_b])
+    da_b = da_b.map_blocks(_get_valid_values, args=[da_a])
 
     # 3. Detrend along the given dim
     demeaned_da_a = da_a - da_a.mean(dim=dim)
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index 22cb03fd..2ab60b89 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -1102,7 +1102,7 @@ class DataArray(AbstractArray, DataWithCoords, DataArrayArithmetic):
             What to do if dimensions that should be selected from are not present in the
             DataArray:
             - "raise": raise an exception
-            - "warning": raise a warning, and ignore the missing dimensions
+            - "warn": raise a warning, and ignore the missing dimensions
             - "ignore": ignore the missing dimensions
         **indexers_kwargs : {dim: indexer, ...}, optional
             The keyword arguments form of ``indexers``.
@@ -2222,7 +2222,7 @@ class DataArray(AbstractArray, DataWithCoords, DataArrayArithmetic):
             What to do if dimensions that should be selected from are not present in the
             DataArray:
             - "raise": raise an exception
-            - "warning": raise a warning, and ignore the missing dimensions
+            - "warn": raise a warning, and ignore the missing dimensions
             - "ignore": ignore the missing dimensions
 
         Returns
@@ -4452,7 +4452,7 @@ class DataArray(AbstractArray, DataWithCoords, DataArrayArithmetic):
             What to do if dimensions that should be selected from are not present in the
             Dataset:
             - "raise": raise an exception
-            - "warning": raise a warning, and ignore the missing dimensions
+            - "warn": raise a warning, and ignore the missing dimensions
             - "ignore": ignore the missing dimensions
         **queries_kwargs : {dim: query, ...}, optional
             The keyword arguments form of ``queries``.
diff --git a/xarray/core/utils.py b/xarray/core/utils.py
index aa51366f..3a9f4180 100644
--- a/xarray/core/utils.py
+++ b/xarray/core/utils.py
@@ -107,6 +107,8 @@ def safe_cast_to_index(array: Any) -> pd.Index:
         index = array
     elif hasattr(array, "to_index"):
         index = array.to_index()
+    elif hasattr(array, "to_pandas_index"):
+        index = array.to_pandas_index()
     else:
         kwargs = {}
         if hasattr(array, "dtype") and array.dtype.kind == "O":


Test patch:
diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py
index cbfa61a4..b7ae1ca9 100644
--- a/xarray/tests/test_computation.py
+++ b/xarray/tests/test_computation.py
@@ -22,7 +22,7 @@ from xarray.core.computation import (
     unified_dim_sizes,
 )
 
-from . import has_dask, requires_dask
+from . import has_dask, raise_if_dask_computes, requires_dask
 
 dask = pytest.importorskip("dask")
 
@@ -1386,6 +1386,7 @@ def arrays_w_tuples():
         da.isel(time=range(2, 20)).rolling(time=3, center=True).mean(),
         xr.DataArray([[1, 2], [1, np.nan]], dims=["x", "time"]),
         xr.DataArray([[1, 2], [np.nan, np.nan]], dims=["x", "time"]),
+        xr.DataArray([[1, 2], [2, 1]], dims=["x", "time"]),
     ]
 
     array_tuples = [
@@ -1394,12 +1395,40 @@ def arrays_w_tuples():
         (arrays[1], arrays[1]),
         (arrays[2], arrays[2]),
         (arrays[2], arrays[3]),
+        (arrays[2], arrays[4]),
+        (arrays[4], arrays[2]),
         (arrays[3], arrays[3]),
+        (arrays[4], arrays[4]),
     ]
 
     return arrays, array_tuples
 
 
+@pytest.mark.parametrize("ddof", [0, 1])
+@pytest.mark.parametrize(
+    "da_a, da_b",
+    [
+        arrays_w_tuples()[1][3],
+        arrays_w_tuples()[1][4],
+        arrays_w_tuples()[1][5],
+        arrays_w_tuples()[1][6],
+        arrays_w_tuples()[1][7],
+        arrays_w_tuples()[1][8],
+    ],
+)
+@pytest.mark.parametrize("dim", [None, "x", "time"])
+def test_lazy_corrcov(da_a, da_b, dim, ddof):
+    # GH 5284
+    from dask import is_dask_collection
+
+    with raise_if_dask_computes():
+        cov = xr.cov(da_a.chunk(), da_b.chunk(), dim=dim, ddof=ddof)
+        assert is_dask_collection(cov)
+
+        corr = xr.corr(da_a.chunk(), da_b.chunk(), dim=dim)
+        assert is_dask_collection(corr)
+
+
 @pytest.mark.parametrize("ddof", [0, 1])
 @pytest.mark.parametrize(
     "da_a, da_b",
diff --git a/xarray/tests/test_utils.py b/xarray/tests/test_utils.py
index 8b5bcbed..483fa723 100644
--- a/xarray/tests/test_utils.py
+++ b/xarray/tests/test_utils.py
@@ -7,6 +7,7 @@ import pytest
 
 from xarray.coding.cftimeindex import CFTimeIndex
 from xarray.core import duck_array_ops, utils
+from xarray.core.indexes import PandasIndex
 from xarray.core.utils import either_dict_or_kwargs
 
 from . import assert_array_equal, requires_cftime, requires_dask
@@ -28,11 +29,13 @@ def test_safe_cast_to_index():
     dates = pd.date_range("2000-01-01", periods=10)
     x = np.arange(5)
     td = x * np.timedelta64(1, "D")
+    midx = pd.MultiIndex.from_tuples([(0,)], names=["a"])
     for expected, array in [
         (dates, dates.values),
         (pd.Index(x, dtype=object), x.astype(object)),
         (pd.Index(td), td),
         (pd.Index(td, dtype=object), td.astype(object)),
+        (midx, PandasIndex(midx)),
     ]:
         actual = utils.safe_cast_to_index(array)
         assert_array_equal(expected, actual)