Regression from 0.10.8 to 0.10.9, "shape mismatch" IndexError in reindex_like() after open_rasterio().squeeze()

[Huite]

Code Sample

import xarray as xr
import numpy as np
import rasterio
from affine import Affine


def write_tif(path, epsg):
    nrow = 5
    ncol = 8
    values = 10.0 * np.random.rand(nrow, ncol)

    profile = dict()
    profile["crs"] = rasterio.crs.CRS.from_epsg(epsg)
    profile["transform"] = Affine.translation(0.0, 0.0)
    profile["driver"] = "GTiff"
    profile["height"] = nrow
    profile["width"] = ncol
    profile["count"] = 1
    profile["dtype"] = np.float64

    with rasterio.Env():
        with rasterio.open(path, "w", **profile) as dst:
            dst.write(values, 1)


write_tif("basic.tif", epsg=28992)
da = xr.open_rasterio("basic.tif").squeeze("band", drop=True)
likevalues = np.empty((10, 16))
likecoords = {"y": np.arange(0.25, 5.0, 0.5), "x": np.arange(0.25, 8.0, 0.5)}
likedims = ("y", "x")
like = xr.DataArray(likevalues, likecoords, likedims)

newda = da.reindex_like(like)

Problem description

The code above executes without issues in xarray version 0.10.8. However, it results in a "shape-mismatch" error in 0.10.9:

Traceback (most recent call last):

  File "<ipython-input-57-f099072e9750>", line 8, in <module>
    newda = da.reindex_like(like)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\dataarray.py", line 919, in reindex_like
    **indexers)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\dataarray.py", line 969, in reindex
    indexers=indexers, method=method, tolerance=tolerance, copy=copy)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\dataset.py", line 1924, in reindex
    tolerance, copy=copy)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\alignment.py", line 377, in reindex_variables
    new_var = var._getitem_with_mask(key)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\variable.py", line 655, in _getitem_with_mask
    data = duck_array_ops.where(mask, fill_value, data)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\duck_array_ops.py", line 183, in where
    return _where(condition, *as_shared_dtype([x, y]))

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\duck_array_ops.py", line 115, in as_shared_dtype
    arrays = [asarray(x) for x in scalars_or_arrays]

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\duck_array_ops.py", line 115, in <listcomp>
    arrays = [asarray(x) for x in scalars_or_arrays]

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\duck_array_ops.py", line 110, in asarray
    return data if isinstance(data, dask_array_type) else np.asarray(data)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\numpy\core\numeric.py", line 492, in asarray
    return array(a, dtype, copy=False, order=order)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\indexing.py", line 624, in __array__
    self._ensure_cached()

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\indexing.py", line 621, in _ensure_cached
    self.array = NumpyIndexingAdapter(np.asarray(self.array))

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\numpy\core\numeric.py", line 492, in asarray
    return array(a, dtype, copy=False, order=order)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\indexing.py", line 602, in __array__
    return np.asarray(self.array, dtype=dtype)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\numpy\core\numeric.py", line 492, in asarray
    return array(a, dtype, copy=False, order=order)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\indexing.py", line 508, in __array__
    return np.asarray(array[self.key], dtype=None)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\backends\rasterio_.py", line 119, in __getitem__
    key, self.shape, indexing.IndexingSupport.OUTER, self._getitem)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\core\indexing.py", line 776, in explicit_indexing_adapter
    result = raw_indexing_method(raw_key.tuple)

  File "d:\bootsma\miniconda3\envs\main\lib\site-packages\xarray\backends\rasterio_.py", line 115, in _getitem
    return out[np_inds]

IndexError: shape mismatch: indexing arrays could not be broadcast together with shapes (10,) (16,) 

I only get it to occur when loading a raster using open_rasterio(). Defining da directly does not lead to issues, using either numpy or dask.array.

values = 10.0 * np.random.rand(1, 5, 8)
coords = {"band": [1.0], "y": np.arange(0.5, 5.0, 1.0), "x": np.arange(0.5, 8.0, 1.0)}
dims = ("band", "y", "x")
da = xr.DataArray(values, coords, dims)
da = da.squeeze("band", drop=True)

import dask.array
values = 10.0 * dask.array.random.random((1, 5, 8), (1, 5, 8))
coords = {"band": [1.0], "y": np.arange(0.5, 5.0, 1.0), "x": np.arange(0.5, 8.0, 1.0)}
dims = ("band", "y", "x")
da = xr.DataArray(values, coords, dims)
da = da.squeeze("band", drop=True)

The "shape mismatch" makes me think that maybe the squeeze() isn't executed properly or in time. This works:

write_tif("basic.tif", epsg=28992)
da = xr.open_rasterio("basic.tif")
likevalues = np.empty((1, 10, 16))
likecoords = {"band":[1], "y": np.arange(0.25, 5.0, 0.5), "x": np.arange(0.25, 8.0, 0.5)}
likedims = ("band", "y", "x")
like = xr.DataArray(likevalues, likecoords, likedims)

newda = da.reindex_like(like)

So it looks to me like some lazy evaluation is going awry when open_rasterio is involved, in this case.
Anything which forces a compute() before reindex_like in this case seems to make the problem go away, like calling da in the console, print(da) after opening the .tif, etc.

write_tif("basic.tif", epsg=28992)
da = xr.open_rasterio("basic.tif").compute().squeeze("band", drop=True)
likevalues = np.empty((10, 16))
likecoords = {"y": np.arange(0.25, 5.0, 0.5), "x": np.arange(0.25, 8.0, 0.5)}
likedims = ("y", "x")
like = xr.DataArray(likevalues, likecoords, likedims)

newda = da.reindex_like(like)

Interestingly, I can't easily inspect it, because inspecting makes the problem go away!

Expected Output

A succesfully reindexed DataArray.

Output of xr.show_versions()

INSTALLED VERSIONS
------------------
commit: None
python: 3.6.5.final.0
python-bits: 64
OS: Windows
OS-release: 7
machine: AMD64
processor: Intel64 Family 6 Model 60 Stepping 3, GenuineIntel
byteorder: little
LC_ALL: None
LANG: en
LOCALE: None.None

xarray: 0.10.9
pandas: 0.23.3
numpy: 1.14.5
scipy: 1.1.0
netCDF4: 1.3.1
h5netcdf: 0.6.1
h5py: 2.8.0
Nio: None
zarr: None
cftime: 1.0.0
PseudonetCDF: None
rasterio: 1.0.0
iris: None
bottleneck: 1.2.1
cyordereddict: None
dask: 0.18.1
distributed: 1.22.0
matplotlib: 2.2.2
cartopy: 0.16.0
seaborn: 0.9.0
setuptools: 40.0.0
pip: 18.0
conda: None
pytest: 3.6.3
IPython: 6.4.0
sphinx: 1.7.5

[shoyer]

Thanks for the clear code to reproduce. I can verify this on the latest development version of xarray.

Anything which forces a compute() before reindex_like in this case seems to make the problem go away, like calling da in the console, print(da) after opening the .tif, etc.

This is because open_rasterio() defaults to cache=True. If you turn off caching, this bug reproduces itself reliably:
da = xr.open_rasterio("basic.tif", cache=False).squeeze("band", drop=True)