Repository: pydata/xarray

Issue:
Xarray with cfgrib backend errors with .where() when drop=True
<!-- 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.
-->

**What happened**:
When loading a HRRR GRIBv2 file in this manner:
`ds = xr.open_dataset(gribfile_path,engine='cfgrib',backend_kwargs={'filter_by_keys':{'typeOfLevel':'hybrid'}})`

I have trouble using the `.where()` method when `drop=True`. If I set `drop=False`, it works fine. I am attempting to subset via latitude and longitude like this:

`ds_sub = ds.where((ds.latitude>=40.0)&(ds.latitude<=50.0)&(ds.longitude>=252.0)&(ds.longitude<=280.0),drop=True)`

However I receive the following errors:
```
Traceback (most recent call last):
  File "read_interp_format_grib.py", line 77, in <module>
    test = ds.where(mask_lat & mask_lon,drop=True)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/common.py", line 1268, in where
    return ops.where_method(self, cond, other)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/ops.py", line 193, in where_method
    return apply_ufunc(
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/computation.py", line 1092, in apply_ufunc
    return apply_dataset_vfunc(
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/computation.py", line 410, in apply_dataset_vfunc
    result_vars = apply_dict_of_variables_vfunc(
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/computation.py", line 356, in apply_dict_of_variables_vfunc
    result_vars[name] = func(*variable_args)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/computation.py", line 606, in apply_variable_ufunc
    input_data = [
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/computation.py", line 607, in <listcomp>
    broadcast_compat_data(arg, broadcast_dims, core_dims)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/computation.py", line 522, in broadcast_compat_data
    data = variable.data
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/variable.py", line 359, in data
    return self.values
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/variable.py", line 510, in values
    return _as_array_or_item(self._data)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/variable.py", line 272, in _as_array_or_item
    data = np.asarray(data)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/numpy/core/_asarray.py", line 83, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/indexing.py", line 685, in __array__
    self._ensure_cached()
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/indexing.py", line 682, in _ensure_cached
    self.array = NumpyIndexingAdapter(np.asarray(self.array))
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/numpy/core/_asarray.py", line 83, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/indexing.py", line 655, in __array__
    return np.asarray(self.array, dtype=dtype)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/numpy/core/_asarray.py", line 83, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/indexing.py", line 560, in __array__
    return np.asarray(array[self.key], dtype=None)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/backends/cfgrib_.py", line 23, in __getitem__
    return indexing.explicit_indexing_adapter(
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/indexing.py", line 848, in explicit_indexing_adapter
    result = NumpyIndexingAdapter(np.asarray(result))[numpy_indices]
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/indexing.py", line 1294, in __getitem__
    return array[key]
IndexError: too many indices for array: array is 2-dimensional, but 3 were indexed
```

**What you expected to happen**:
I expect the dataset to be reduced in the x and y (latitude/longitude) dimensions where `cond=False` from `.where()` when `drop=True`.

**Minimal Complete Verifiable Example**:

```python
# Put your MCVE code here
```

**Anything else we need to know?**:
I was able to confirm cfgrib is where the issue lies by doing the following:
```
ds = xr.open_dataset(gribfile_path,engine='cfgrib',backend_kwargs={'filter_by_keys':{'typeOfLevel':'hybrid'}})
ds.to_netcdf('test.nc')
dsnc = xr.open_dataset('test.nc')
ds_sub = dsnc.where((dsnc.latitude>=40.0)&(dsnc.latitude<=50.0)&(dsnc.longitude>=252.0)&(dsnc.longitude<=280.0),drop=True)
```

That correctly gives me:
`Dimensions:     (hybrid: 50, x: 792, y: 414)`

Originally x=1799 and y=1059.

**Environment**:

<details><summary>Output of <tt>xr.show_versions()</tt></summary>

INSTALLED VERSIONS
------------------
commit: None
python: 3.8.6 | packaged by conda-forge | (default, Oct  7 2020, 19:08:05) 
[GCC 7.5.0]
python-bits: 64
OS: Linux
OS-release: 4.9.0-14-amd64
machine: x86_64
processor: 
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: en_US.UTF-8
libhdf5: 1.10.6
libnetcdf: 4.7.4

xarray: 0.16.1
pandas: 1.1.3
numpy: 1.19.1
scipy: 1.5.2
netCDF4: 1.5.5.1
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.2.1
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: 0.9.8.5
iris: None
bottleneck: None
dask: 2.30.0
distributed: 2.30.0
matplotlib: 3.3.2
cartopy: 0.17.0
seaborn: None
numbagg: None
pint: 0.16.1
setuptools: 49.6.0.post20200917
pip: 20.2.3
conda: None
pytest: None
IPython: None
sphinx: None


</details>


Discussion:

Patch:
diff --git a/doc/conf.py b/doc/conf.py
index fcd301ba..d83e966f 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -115,7 +115,7 @@ napoleon_google_docstring = False
 napoleon_numpy_docstring = True
 
 napoleon_use_param = False
-napoleon_use_rtype = True
+napoleon_use_rtype = False
 napoleon_preprocess_types = True
 napoleon_type_aliases = {
     # general terms
@@ -418,6 +418,7 @@ intersphinx_mapping = {
     "matplotlib": ("https://matplotlib.org", None),
     "dask": ("https://docs.dask.org/en/latest", None),
     "cftime": ("https://unidata.github.io/cftime", None),
+    "rasterio": ("https://rasterio.readthedocs.io/en/latest", None),
     "sparse": ("https://sparse.pydata.org/en/latest/", None),
 }
 
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 00bceed4..151af2de 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -40,6 +40,8 @@ Bug fixes
   By `Richard Kleijn <https://github.com/rhkleijn>`_ .
 - Remove dictionary unpacking when using ``.loc`` to avoid collision with ``.sel`` parameters (:pull:`4695`).
   By `Anderson Banihirwe <https://github.com/andersy005>`_
+- Fix a crash in orthogonal indexing on geographic coordinates with ``engine='cfgrib'`` (:issue:`4733` :pull:`4737`).
+  By `Alessandro Amici <https://github.com/alexamici>`_
 
 Documentation
 ~~~~~~~~~~~~~
@@ -115,7 +117,7 @@ Bug fixes
 ~~~~~~~~~
 
 - Fix bug where reference times without padded years (e.g. ``since 1-1-1``) would lose their units when
-  being passed by :py:func:`encode_cf_datetime` (:issue:`4422`, :pull:`4506`). Such units are ambiguous
+  being passed by ``encode_cf_datetime`` (:issue:`4422`, :pull:`4506`). Such units are ambiguous
   about which digit represents the years (is it YMD or DMY?). Now, if such formatting is encountered,
   it is assumed that the first digit is the years, they are padded appropriately (to e.g. ``since 0001-1-1``)
   and a warning that this assumption is being made is issued. Previously, without ``cftime``, such times
@@ -319,8 +321,10 @@ Internal Changes
 
   All are up from 6 months (:issue:`4295`)
   `Guido Imperiale <https://github.com/crusaderky>`_.
-- Use :py:func:`dask.array.apply_gufunc` instead of :py:func:`dask.array.blockwise` in
-  :py:func:`xarray.apply_ufunc` when using ``dask='parallelized'``. (:pull:`4060`, :pull:`4391`, :pull:`4392`)
+- Use :py:func:`dask.array.apply_gufunc <dask.array.gufunc.apply_gufunc>` instead of
+  :py:func:`dask.array.blockwise` in :py:func:`xarray.apply_ufunc` when using
+  ``dask='parallelized'``. (:pull:`4060`, :pull:`4391`, :pull:`4392`)
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
 - Align ``mypy`` versions to ``0.782`` across ``requirements`` and
   ``.pre-commit-config.yml`` files. (:pull:`4390`)
   By `Maximilian Roos <https://github.com/max-sixty>`_
diff --git a/xarray/backends/cfgrib_.py b/xarray/backends/cfgrib_.py
index 57cc4d6f..d4933e37 100644
--- a/xarray/backends/cfgrib_.py
+++ b/xarray/backends/cfgrib_.py
@@ -24,7 +24,7 @@ class CfGribArrayWrapper(BackendArray):
 
     def __getitem__(self, key):
         return indexing.explicit_indexing_adapter(
-            key, self.shape, indexing.IndexingSupport.OUTER, self._getitem
+            key, self.shape, indexing.IndexingSupport.BASIC, self._getitem
         )
 
     def _getitem(self, key):
diff --git a/xarray/core/combine.py b/xarray/core/combine.py
index 38a32758..86ed1870 100644
--- a/xarray/core/combine.py
+++ b/xarray/core/combine.py
@@ -518,7 +518,6 @@ def combine_nested(
     --------
     concat
     merge
-    auto_combine
     """
     if isinstance(concat_dim, (str, DataArray)) or concat_dim is None:
         concat_dim = [concat_dim]
diff --git a/xarray/core/common.py b/xarray/core/common.py
index 0a626361..302a8e52 100644
--- a/xarray/core/common.py
+++ b/xarray/core/common.py
@@ -375,8 +375,7 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
         try:
             return self.indexes[key]
         except KeyError:
-            # need to ensure dtype=int64 in case range is empty on Python 2
-            return pd.Index(range(self.sizes[key]), name=key, dtype=np.int64)
+            return pd.Index(range(self.sizes[key]), name=key)
 
     def _calc_assign_results(
         self: C, kwargs: Mapping[Hashable, Union[T, Callable[[C], T]]]
@@ -817,8 +816,9 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
 
         Returns
         -------
-        Rolling object (core.rolling.DataArrayRolling for DataArray,
-        core.rolling.DatasetRolling for Dataset.)
+        core.rolling.DataArrayRolling or core.rolling.DatasetRolling
+            A rolling object (``DataArrayRolling`` for ``DataArray``,
+            ``DatasetRolling`` for ``Dataset``)
 
         Examples
         --------
@@ -929,8 +929,9 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
 
         Returns
         -------
-        Coarsen object (core.rolling.DataArrayCoarsen for DataArray,
-        core.rolling.DatasetCoarsen for Dataset.)
+        core.rolling.DataArrayCoarsen or core.rolling.DatasetCoarsen
+            A coarsen object (``DataArrayCoarsen`` for ``DataArray``,
+            ``DatasetCoarsen`` for ``Dataset``)
 
         Examples
         --------
diff --git a/xarray/core/computation.py b/xarray/core/computation.py
index c332ecb0..e0d9ff4b 100644
--- a/xarray/core/computation.py
+++ b/xarray/core/computation.py
@@ -1148,24 +1148,24 @@ def cov(da_a, da_b, dim=None, ddof=1):
 
     Parameters
     ----------
-    da_a: DataArray
+    da_a : DataArray
         Array to compute.
-    da_b: DataArray
+    da_b : DataArray
         Array to compute.
     dim : str, optional
         The dimension along which the covariance will be computed
-    ddof: int, optional
+    ddof : int, optional
         If ddof=1, covariance is normalized by N-1, giving an unbiased estimate,
         else normalization is by N.
 
     Returns
     -------
-    covariance: DataArray
+    covariance : DataArray
 
     See also
     --------
-    pandas.Series.cov: corresponding pandas function
-    xr.corr: respective function to calculate correlation
+    pandas.Series.cov : corresponding pandas function
+    xarray.corr: respective function to calculate correlation
 
     Examples
     --------
@@ -1229,11 +1229,11 @@ def corr(da_a, da_b, dim=None):
 
     Parameters
     ----------
-    da_a: DataArray
+    da_a : DataArray
         Array to compute.
-    da_b: DataArray
+    da_b : DataArray
         Array to compute.
-    dim: str, optional
+    dim : str, optional
         The dimension along which the correlation will be computed
 
     Returns
@@ -1242,8 +1242,8 @@ def corr(da_a, da_b, dim=None):
 
     See also
     --------
-    pandas.Series.corr: corresponding pandas function
-    xr.cov: underlying covariance function
+    pandas.Series.corr : corresponding pandas function
+    xarray.cov : underlying covariance function
 
     Examples
     --------
diff --git a/xarray/core/concat.py b/xarray/core/concat.py
index a4b9065e..1275d002 100644
--- a/xarray/core/concat.py
+++ b/xarray/core/concat.py
@@ -160,7 +160,6 @@ def concat(
     See also
     --------
     merge
-    auto_combine
 
     Examples
     --------
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index 6316afd5..489642d0 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -856,11 +856,11 @@ class DataArray(AbstractArray, DataWithCoords):
         Parameters
         ----------
         **kwargs : dict
-            Additional keyword arguments passed on to ``dask.array.compute``.
+            Additional keyword arguments passed on to ``dask.compute``.
 
         See Also
         --------
-        dask.array.compute
+        dask.compute
         """
         ds = self._to_temp_dataset().load(**kwargs)
         new = self._from_temp_dataset(ds)
@@ -881,11 +881,11 @@ class DataArray(AbstractArray, DataWithCoords):
         Parameters
         ----------
         **kwargs : dict
-            Additional keyword arguments passed on to ``dask.array.compute``.
+            Additional keyword arguments passed on to ``dask.compute``.
 
         See Also
         --------
-        dask.array.compute
+        dask.compute
         """
         new = self.copy(deep=False)
         return new.load(**kwargs)
@@ -2994,10 +2994,10 @@ class DataArray(AbstractArray, DataWithCoords):
         difference : same type as caller
             The n-th order finite difference of this object.
 
-        .. note::
-
-            `n` matches numpy's behavior and is different from pandas' first
-            argument named `periods`.
+        Notes
+        -----
+        `n` matches numpy's behavior and is different from pandas' first argument named
+        `periods`.
 
 
         Examples
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 2476cb26..3bb5cd8b 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -779,11 +779,11 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         Parameters
         ----------
         **kwargs : dict
-            Additional keyword arguments passed on to ``dask.array.compute``.
+            Additional keyword arguments passed on to ``dask.compute``.
 
         See Also
         --------
-        dask.array.compute
+        dask.compute
         """
         # access .data to coerce everything to numpy or dask arrays
         lazy_data = {
@@ -951,11 +951,11 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         Parameters
         ----------
         **kwargs : dict
-            Additional keyword arguments passed on to ``dask.array.compute``.
+            Additional keyword arguments passed on to ``dask.compute``.
 
         See Also
         --------
-        dask.array.compute
+        dask.compute
         """
         new = self.copy(deep=False)
         return new.load(**kwargs)
diff --git a/xarray/core/rolling.py b/xarray/core/rolling.py
index 38cb11b5..39d88924 100644
--- a/xarray/core/rolling.py
+++ b/xarray/core/rolling.py
@@ -41,10 +41,10 @@ class Rolling:
 
     See Also
     --------
-    Dataset.groupby
-    DataArray.groupby
-    Dataset.rolling
-    DataArray.rolling
+    xarray.Dataset.groupby
+    xarray.DataArray.groupby
+    xarray.Dataset.rolling
+    xarray.DataArray.rolling
     """
 
     __slots__ = ("obj", "window", "min_periods", "center", "dim", "keep_attrs")
@@ -209,10 +209,10 @@ class DataArrayRolling(Rolling):
 
         See Also
         --------
-        DataArray.rolling
-        DataArray.groupby
-        Dataset.rolling
-        Dataset.groupby
+        xarray.DataArray.rolling
+        xarray.DataArray.groupby
+        xarray.Dataset.rolling
+        xarray.Dataset.groupby
         """
         super().__init__(
             obj, windows, min_periods=min_periods, center=center, keep_attrs=keep_attrs
@@ -526,10 +526,10 @@ class DatasetRolling(Rolling):
 
         See Also
         --------
-        Dataset.rolling
-        DataArray.rolling
-        Dataset.groupby
-        DataArray.groupby
+        xarray.Dataset.rolling
+        xarray.DataArray.rolling
+        xarray.Dataset.groupby
+        xarray.DataArray.groupby
         """
         super().__init__(obj, windows, min_periods, center, keep_attrs)
         if any(d not in self.obj.dims for d in self.dim):


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index b13b86fc..2d20f7f6 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -3575,6 +3575,19 @@ class TestCfGrib:
             assert list(ds.data_vars) == ["t"]
             assert ds["t"].min() == 231.0
 
+    def test_read_outer(self):
+        expected = {
+            "number": 2,
+            "time": 3,
+            "isobaricInhPa": 2,
+            "latitude": 2,
+            "longitude": 3,
+        }
+        with open_example_dataset("example.grib", engine="cfgrib") as ds:
+            res = ds.isel(latitude=[0, 2], longitude=[0, 1, 2])
+            assert res.dims == expected
+            assert res["t"].min() == 231.0
+
 
 @requires_pseudonetcdf
 @pytest.mark.filterwarnings("ignore:IOAPI_ISPH is assumed to be 6370000")
diff --git a/xarray/tests/test_interp.py b/xarray/tests/test_interp.py
index 5f79f81d..20d5fb12 100644
--- a/xarray/tests/test_interp.py
+++ b/xarray/tests/test_interp.py
@@ -816,7 +816,7 @@ def test_interpolate_chunk_1d(method, data_ndim, interp_ndim, nscalar, chunked):
 
                 assert_identical(actual, expected)
 
-                # all the combinations are usualy not necessary
+                # all the combinations are usually not necessary
                 break
             break
         break