Repository: pydata/xarray

Issue:
keepdims=True for xarray reductions
For operations where arrays are aggregated but then combined, the `keepdims=True` option for NumPy aggregations is convenient.

We should consider supporting this in xarray as well. Aggregating a DataArray/Dataset with `keepdims=True` (or maybe `keep_dims=True`) would remove all original coordinates along aggregated dimensions and return a result with a dimension of size 1 without any coordinates, e.g.,
```
>>> array = xr.DataArray([1, 2, 3], dims='x', coords={'x': ['a', 'b', 'c']})
>>> array.mean(keepdims=True)
<xarray.DataArray (x: 1)>
array([2.])
Dimensions without coordinates: x
```

In case, `array.mean(keepdims=True()` is equivalent to `array.mean().expand_dims('x')` but in general this equivalent does not hold, because the location of the original dimension is lost.

Implementation-wise, we have two options:
1. Pass on `keepdims=True` to NumPy functions like `numpy.mean()`, or
2. Implement `keepdims=True` ourselves, in `Variable.reduce()`.

I think I like option 2 a little better, because it places fewer requirements on aggregation functions. For example, functions like `bottleneck.nanmean()` don't accept a `keepdims` argument.

Discussion:
pep8speaks: Hello @ScottWales! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:

* In the file [`xarray/core/dataarray.py`](https://github.com/pydata/xarray/blob/bd2dc61a8ff1062c5daf03c77f633d7e3219468c/xarray/core/dataarray.py):

> [Line 1646:80](https://github.com/pydata/xarray/blob/bd2dc61a8ff1062c5daf03c77f633d7e3219468c/xarray/core/dataarray.py#L1646): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (80 > 79 characters)












ScottWales: Thanks for the review @max-sixty, happy to help out

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index e62c7e87..373cb8d1 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -22,6 +22,8 @@ Enhancements
 ~~~~~~~~~~~~
 
 
+- Add ``keepdims`` argument for reduce operations (:issue:`2170`)
+  By `Scott Wales <https://github.com/ScottWales>`_.
 - netCDF chunksizes are now only dropped when original_shape is different,
   not when it isn't found. (:issue:`2207`)
   By `Karel van de Plassche <https://github.com/Karel-van-de-Plassche>`_.
@@ -59,6 +61,11 @@ Enhancements
   formatted datetimes. By `Alan Brammer <https://github.com/abrammer>`_.
 - Add ``.str`` accessor to DataArrays for string related manipulations.
   By `0x0L <https://github.com/0x0L>`_.
+- Add ``errors`` keyword argument to :py:meth:`Dataset.drop` and :py:meth:`Dataset.drop_dims`
+  that allows ignoring errors if a passed label or dimension is not in the dataset
+  (:issue:`2994`).
+  By `Andrew Ross <https://github.com/andrew-c-ross>`_.
+
 
 Bug fixes
 ~~~~~~~~~
@@ -70,8 +77,6 @@ Bug fixes
   By `Mayeul d'Avezac <https://github.com/mdavezac>`_.
 - Return correct count for scalar datetime64 arrays (:issue:`2770`)
   By `Dan Nowacki <https://github.com/dnowacki-usgs>`_.
-- Fix facetgrid colormap bug when ``extend=True``. (:issue:`2932`)
-  By `Deepak Cherian <https://github.com/dcherian`_.
 - A deep copy deep-copies the coords (:issue:`1463`)
   By `Martin Pletcher <https://github.com/pletchm>`_.
 - Increased support for `missing_value` (:issue:`2871`)
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index 2746c32a..ff77a6ab 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -259,8 +259,14 @@ class DataArray(AbstractArray, DataWithCoords):
         return type(self)(variable, coords, name=name, fastpath=True)
 
     def _replace_maybe_drop_dims(self, variable, name=__default):
-        if variable.dims == self.dims:
+        if variable.dims == self.dims and variable.shape == self.shape:
             coords = self._coords.copy()
+        elif variable.dims == self.dims:
+            # Shape has changed (e.g. from reduce(..., keepdims=True)
+            new_sizes = dict(zip(self.dims, variable.shape))
+            coords = OrderedDict((k, v) for k, v in self._coords.items()
+                                 if v.shape == tuple(new_sizes[d]
+                                                     for d in v.dims))
         else:
             allowed_dims = set(variable.dims)
             coords = OrderedDict((k, v) for k, v in self._coords.items()
@@ -1462,7 +1468,7 @@ class DataArray(AbstractArray, DataWithCoords):
     def T(self) -> 'DataArray':
         return self.transpose()
 
-    def drop(self, labels, dim=None):
+    def drop(self, labels, dim=None, *, errors='raise'):
         """Drop coordinates or index labels from this DataArray.
 
         Parameters
@@ -1472,14 +1478,18 @@ class DataArray(AbstractArray, DataWithCoords):
         dim : str, optional
             Dimension along which to drop index labels. By default (if
             ``dim is None``), drops coordinates rather than index labels.
-
+        errors: {'raise', 'ignore'}, optional
+            If 'raise' (default), raises a ValueError error if
+            any of the coordinates or index labels passed are not
+            in the array. If 'ignore', any given labels that are in the
+            array are dropped and no error is raised.
         Returns
         -------
         dropped : DataArray
         """
         if utils.is_scalar(labels):
             labels = [labels]
-        ds = self._to_temp_dataset().drop(labels, dim)
+        ds = self._to_temp_dataset().drop(labels, dim, errors=errors)
         return self._from_temp_dataset(ds)
 
     def dropna(self, dim, how='any', thresh=None):
@@ -1638,7 +1648,8 @@ class DataArray(AbstractArray, DataWithCoords):
         """
         return ops.fillna(self, other, join="outer")
 
-    def reduce(self, func, dim=None, axis=None, keep_attrs=None, **kwargs):
+    def reduce(self, func, dim=None, axis=None, keep_attrs=None,
+               keepdims=False, **kwargs):
         """Reduce this array by applying `func` along some dimension(s).
 
         Parameters
@@ -1658,6 +1669,10 @@ class DataArray(AbstractArray, DataWithCoords):
             If True, the variable's attributes (`attrs`) will be copied from
             the original object to the new one.  If False (default), the new
             object will be returned without attributes.
+        keepdims : bool, default False
+            If True, the dimensions which are reduced are left in the result
+            as dimensions of size one. Coordinates that use these dimensions
+            are removed.
         **kwargs : dict
             Additional keyword arguments passed on to `func`.
 
@@ -1668,7 +1683,8 @@ class DataArray(AbstractArray, DataWithCoords):
             summarized data and the indicated dimension(s) removed.
         """
 
-        var = self.variable.reduce(func, dim, axis, keep_attrs, **kwargs)
+        var = self.variable.reduce(func, dim, axis, keep_attrs, keepdims,
+                                   **kwargs)
         return self._replace_maybe_drop_dims(var)
 
     def to_pandas(self):
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 026be5ba..3e00640b 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -2826,7 +2826,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
             raise ValueError('One or more of the specified variables '
                              'cannot be found in this dataset')
 
-    def drop(self, labels, dim=None):
+    def drop(self, labels, dim=None, *, errors='raise'):
         """Drop variables or index labels from this dataset.
 
         Parameters
@@ -2836,33 +2836,41 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         dim : None or str, optional
             Dimension along which to drop index labels. By default (if
             ``dim is None``), drops variables rather than index labels.
+        errors: {'raise', 'ignore'}, optional
+            If 'raise' (default), raises a ValueError error if
+            any of the variable or index labels passed are not
+            in the dataset. If 'ignore', any given labels that are in the
+            dataset are dropped and no error is raised.
 
         Returns
         -------
         dropped : Dataset
         """
+        if errors not in ['raise', 'ignore']:
+            raise ValueError('errors must be either "raise" or "ignore"')
         if utils.is_scalar(labels):
             labels = [labels]
         if dim is None:
-            return self._drop_vars(labels)
+            return self._drop_vars(labels, errors=errors)
         else:
             try:
                 index = self.indexes[dim]
             except KeyError:
                 raise ValueError(
                     'dimension %r does not have coordinate labels' % dim)
-            new_index = index.drop(labels)
+            new_index = index.drop(labels, errors=errors)
             return self.loc[{dim: new_index}]
 
-    def _drop_vars(self, names):
-        self._assert_all_in_dataset(names)
+    def _drop_vars(self, names, errors='raise'):
+        if errors == 'raise':
+            self._assert_all_in_dataset(names)
         drop = set(names)
         variables = OrderedDict((k, v) for k, v in self._variables.items()
                                 if k not in drop)
         coord_names = set(k for k in self._coord_names if k in variables)
         return self._replace_vars_and_dims(variables, coord_names)
 
-    def drop_dims(self, drop_dims):
+    def drop_dims(self, drop_dims, *, errors='raise'):
         """Drop dimensions and associated variables from this dataset.
 
         Parameters
@@ -2875,14 +2883,23 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         obj : Dataset
             The dataset without the given dimensions (or any variables
             containing those dimensions)
+        errors: {'raise', 'ignore'}, optional
+            If 'raise' (default), raises a ValueError error if
+            any of the dimensions passed are not
+            in the dataset. If 'ignore', any given dimensions that are in the
+            dataset are dropped and no error is raised.
         """
+        if errors not in ['raise', 'ignore']:
+            raise ValueError('errors must be either "raise" or "ignore"')
+
         if utils.is_scalar(drop_dims):
             drop_dims = [drop_dims]
 
-        missing_dimensions = [d for d in drop_dims if d not in self.dims]
-        if missing_dimensions:
-            raise ValueError('Dataset does not contain the dimensions: %s'
-                             % missing_dimensions)
+        if errors == 'raise':
+            missing_dimensions = [d for d in drop_dims if d not in self.dims]
+            if missing_dimensions:
+                raise ValueError('Dataset does not contain the dimensions: %s'
+                                 % missing_dimensions)
 
         drop_vars = set(k for k, v in self._variables.items()
                         for d in v.dims if d in drop_dims)
@@ -3135,8 +3152,8 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         out = ops.fillna(self, other, join="outer", dataset_join="outer")
         return out
 
-    def reduce(self, func, dim=None, keep_attrs=None, numeric_only=False,
-               allow_lazy=False, **kwargs):
+    def reduce(self, func, dim=None, keep_attrs=None, keepdims=False,
+               numeric_only=False, allow_lazy=False, **kwargs):
         """Reduce this dataset by applying `func` along some dimension(s).
 
         Parameters
@@ -3152,6 +3169,10 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
             If True, the dataset's attributes (`attrs`) will be copied from
             the original object to the new one.  If False (default), the new
             object will be returned without attributes.
+        keepdims : bool, default False
+            If True, the dimensions which are reduced are left in the result
+            as dimensions of size one. Coordinates that use these dimensions
+            are removed.
         numeric_only : bool, optional
             If True, only apply ``func`` to variables with a numeric dtype.
         **kwargs : dict
@@ -3201,6 +3222,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
                         reduce_dims = None
                     variables[name] = var.reduce(func, dim=reduce_dims,
                                                  keep_attrs=keep_attrs,
+                                                 keepdims=keepdims,
                                                  allow_lazy=allow_lazy,
                                                  **kwargs)
 
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index 41f8795b..ab1be181 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -1334,7 +1334,7 @@ class Variable(common.AbstractArray, arithmetic.SupportsArithmetic,
         return ops.where_method(self, cond, other)
 
     def reduce(self, func, dim=None, axis=None,
-               keep_attrs=None, allow_lazy=False, **kwargs):
+               keep_attrs=None, keepdims=False, allow_lazy=False, **kwargs):
         """Reduce this array by applying `func` along some dimension(s).
 
         Parameters
@@ -1354,6 +1354,9 @@ class Variable(common.AbstractArray, arithmetic.SupportsArithmetic,
             If True, the variable's attributes (`attrs`) will be copied from
             the original object to the new one.  If False (default), the new
             object will be returned without attributes.
+        keepdims : bool, default False
+            If True, the dimensions which are reduced are left in the result
+            as dimensions of size one
         **kwargs : dict
             Additional keyword arguments passed on to `func`.
 
@@ -1381,8 +1384,19 @@ class Variable(common.AbstractArray, arithmetic.SupportsArithmetic,
         else:
             removed_axes = (range(self.ndim) if axis is None
                             else np.atleast_1d(axis) % self.ndim)
-            dims = [adim for n, adim in enumerate(self.dims)
-                    if n not in removed_axes]
+            if keepdims:
+                # Insert np.newaxis for removed dims
+                slices = tuple(np.newaxis if i in removed_axes else
+                               slice(None, None) for i in range(self.ndim))
+                if getattr(data, 'shape', None) is None:
+                    # Reduce has produced a scalar value, not an array-like
+                    data = np.asanyarray(data)[slices]
+                else:
+                    data = data[slices]
+                dims = self.dims
+            else:
+                dims = [adim for n, adim in enumerate(self.dims)
+                        if n not in removed_axes]
 
         if keep_attrs is None:
             keep_attrs = _get_keep_attrs(default=False)
diff --git a/xarray/plot/utils.py b/xarray/plot/utils.py
index 18215479..c9f72b17 100644
--- a/xarray/plot/utils.py
+++ b/xarray/plot/utils.py
@@ -264,8 +264,7 @@ def _determine_cmap_params(plot_data, vmin=None, vmax=None, cmap=None,
     if extend is None:
         extend = _determine_extend(calc_data, vmin, vmax)
 
-    if ((levels is not None or isinstance(norm, mpl.colors.BoundaryNorm))
-            and (not isinstance(cmap, mpl.colors.Colormap))):
+    if levels is not None or isinstance(norm, mpl.colors.BoundaryNorm):
         cmap, newnorm = _build_discrete_cmap(cmap, levels, extend, filled)
         norm = newnorm if norm is None else norm
 


Test patch:
diff --git a/xarray/tests/test_dataarray.py b/xarray/tests/test_dataarray.py
index fd9076e7..47222194 100644
--- a/xarray/tests/test_dataarray.py
+++ b/xarray/tests/test_dataarray.py
@@ -1859,19 +1859,34 @@ class TestDataArray:
         with pytest.raises(ValueError):
             arr.drop('not found')
 
+        actual = expected.drop('not found', errors='ignore')
+        assert_identical(actual, expected)
+
         with raises_regex(ValueError, 'cannot be found'):
             arr.drop(None)
 
+        actual = expected.drop(None, errors='ignore')
+        assert_identical(actual, expected)
+
         renamed = arr.rename('foo')
         with raises_regex(ValueError, 'cannot be found'):
             renamed.drop('foo')
 
+        actual = renamed.drop('foo', errors='ignore')
+        assert_identical(actual, renamed)
+
     def test_drop_index_labels(self):
         arr = DataArray(np.random.randn(2, 3), coords={'y': [0, 1, 2]},
                         dims=['x', 'y'])
         actual = arr.drop([0, 1], dim='y')
         expected = arr[:, 2:]
-        assert_identical(expected, actual)
+        assert_identical(actual, expected)
+
+        with raises_regex((KeyError, ValueError), 'not .* in axis'):
+            actual = arr.drop([0, 1, 3], dim='y')
+
+        actual = arr.drop([0, 1, 3], dim='y', errors='ignore')
+        assert_identical(actual, expected)
 
     def test_dropna(self):
         x = np.random.randn(4, 4)
@@ -1976,6 +1991,44 @@ class TestDataArray:
                              dims=['x', 'y']).mean('x')
         assert_equal(actual, expected)
 
+    def test_reduce_keepdims(self):
+        coords = {'x': [-1, -2], 'y': ['ab', 'cd', 'ef'],
+                  'lat': (['x', 'y'], [[1, 2, 3], [-1, -2, -3]]),
+                  'c': -999}
+        orig = DataArray([[-1, 0, 1], [-3, 0, 3]], coords, dims=['x', 'y'])
+
+        # Mean on all axes loses non-constant coordinates
+        actual = orig.mean(keepdims=True)
+        expected = DataArray(orig.data.mean(keepdims=True), dims=orig.dims,
+                             coords={k: v for k, v in coords.items()
+                                     if k in ['c']})
+        assert_equal(actual, expected)
+
+        assert actual.sizes['x'] == 1
+        assert actual.sizes['y'] == 1
+
+        # Mean on specific axes loses coordinates not involving that axis
+        actual = orig.mean('y', keepdims=True)
+        expected = DataArray(orig.data.mean(axis=1, keepdims=True),
+                             dims=orig.dims,
+                             coords={k: v for k, v in coords.items()
+                                     if k not in ['y', 'lat']})
+        assert_equal(actual, expected)
+
+    @requires_bottleneck
+    def test_reduce_keepdims_bottleneck(self):
+        import bottleneck
+
+        coords = {'x': [-1, -2], 'y': ['ab', 'cd', 'ef'],
+                  'lat': (['x', 'y'], [[1, 2, 3], [-1, -2, -3]]),
+                  'c': -999}
+        orig = DataArray([[-1, 0, 1], [-3, 0, 3]], coords, dims=['x', 'y'])
+
+        # Bottleneck does not have its own keepdims implementation
+        actual = orig.reduce(bottleneck.nanmean, keepdims=True)
+        expected = orig.mean(keepdims=True)
+        assert_equal(actual, expected)
+
     def test_reduce_dtype(self):
         coords = {'x': [-1, -2], 'y': ['ab', 'cd', 'ef'],
                   'lat': (['x', 'y'], [[1, 2, 3], [-1, -2, -3]]),
diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py
index 5aae5648..e3a01bbd 100644
--- a/xarray/tests/test_dataset.py
+++ b/xarray/tests/test_dataset.py
@@ -1889,6 +1889,15 @@ class TestDataset:
         with raises_regex(ValueError, 'cannot be found'):
             data.drop('not_found_here')
 
+        actual = data.drop('not_found_here', errors='ignore')
+        assert_identical(data, actual)
+
+        actual = data.drop(['not_found_here'], errors='ignore')
+        assert_identical(data, actual)
+
+        actual = data.drop(['time', 'not_found_here'], errors='ignore')
+        assert_identical(expected, actual)
+
     def test_drop_index_labels(self):
         data = Dataset({'A': (['x', 'y'], np.random.randn(2, 3)),
                         'x': ['a', 'b']})
@@ -1907,6 +1916,16 @@ class TestDataset:
             # not contained in axis
             data.drop(['c'], dim='x')
 
+        actual = data.drop(['c'], dim='x', errors='ignore')
+        assert_identical(data, actual)
+
+        with pytest.raises(ValueError):
+            data.drop(['c'], dim='x', errors='wrong_value')
+
+        actual = data.drop(['a', 'b', 'c'], 'x', errors='ignore')
+        expected = data.isel(x=slice(0, 0))
+        assert_identical(expected, actual)
+
         with raises_regex(
                 ValueError, 'does not have coordinate labels'):
             data.drop(1, 'y')
@@ -1931,6 +1950,22 @@ class TestDataset:
         with pytest.raises((ValueError, KeyError)):
             data.drop_dims('z')  # not a dimension
 
+        with pytest.raises((ValueError, KeyError)):
+            data.drop_dims(None)
+
+        actual = data.drop_dims('z', errors='ignore')
+        assert_identical(data, actual)
+
+        actual = data.drop_dims(None, errors='ignore')
+        assert_identical(data, actual)
+
+        with pytest.raises(ValueError):
+            actual = data.drop_dims('z', errors='wrong_value')
+
+        actual = data.drop_dims(['x', 'y', 'z'], errors='ignore')
+        expected = data.drop(['A', 'B', 'x'])
+        assert_identical(expected, actual)
+
     def test_copy(self):
         data = create_test_data()
         data.attrs['Test'] = [1, 2, 3]
@@ -3863,6 +3898,25 @@ class TestDataset:
         with raises_regex(TypeError, "unexpected keyword argument 'axis'"):
             ds.reduce(total_sum, dim='x')
 
+    def test_reduce_keepdims(self):
+        ds = Dataset({'a': (['x', 'y'], [[0, 1, 2, 3, 4]])},
+                     coords={'y': [0, 1, 2, 3, 4], 'x': [0],
+                             'lat': (['x', 'y'], [[0, 1, 2, 3, 4]]),
+                             'c': -999.0})
+
+        # Shape should match behaviour of numpy reductions with keepdims=True
+        # Coordinates involved in the reduction should be removed
+        actual = ds.mean(keepdims=True)
+        expected = Dataset({'a': (['x', 'y'], np.mean(ds.a, keepdims=True))},
+                           coords={'c': ds.c})
+        assert_identical(expected, actual)
+
+        actual = ds.mean('x', keepdims=True)
+        expected = Dataset({'a': (['x', 'y'],
+                                  np.mean(ds.a, axis=0, keepdims=True))},
+                           coords={'y': ds.y, 'c': ds.c})
+        assert_identical(expected, actual)
+
     def test_quantile(self):
 
         ds = create_test_data(seed=123)
diff --git a/xarray/tests/test_plot.py b/xarray/tests/test_plot.py
index a0952cac..0dc5fb32 100644
--- a/xarray/tests/test_plot.py
+++ b/xarray/tests/test_plot.py
@@ -8,6 +8,7 @@ import pytest
 import xarray as xr
 import xarray.plot as xplt
 from xarray import DataArray
+from xarray.coding.times import _import_cftime
 from xarray.plot.plot import _infer_interval_breaks
 from xarray.plot.utils import (
     _build_discrete_cmap, _color_palette, _determine_cmap_params,
@@ -537,25 +538,6 @@ class TestDetermineCmapParams:
             cmap_params = _determine_cmap_params(self.data)
             assert cmap_params['cmap'] == 'magma'
 
-    def test_do_nothing_if_provided_cmap(self):
-        cmap_list = [
-            mpl.colors.LinearSegmentedColormap.from_list('name', ['r', 'g']),
-            mpl.colors.ListedColormap(['r', 'g', 'b'])
-        ]
-
-        # can't parametrize with mpl objects when mpl is absent
-        for cmap in cmap_list:
-            cmap_params = _determine_cmap_params(self.data,
-                                                 cmap=cmap,
-                                                 levels=7)
-            assert cmap_params['cmap'] is cmap
-
-    def test_do_something_if_provided_str_cmap(self):
-        cmap = 'RdBu_r'
-        cmap_params = _determine_cmap_params(self.data, cmap=cmap, levels=7)
-        assert cmap_params['cmap'] is not cmap
-        assert isinstance(cmap_params['cmap'], mpl.colors.ListedColormap)
-
     def test_cmap_sequential_explicit_option(self):
         with xr.set_options(cmap_sequential=mpl.cm.magma):
             cmap_params = _determine_cmap_params(self.data)
@@ -775,13 +757,14 @@ class TestDiscreteColorMap:
 
     @pytest.mark.slow
     def test_discrete_colormap_int_levels(self):
-        for extend, levels, vmin, vmax in [('neither', 7, None, None),
-                                           ('neither', 7, None, 20),
-                                           ('both', 7, 4, 8),
-                                           ('min', 10, 4, 15)]:
+        for extend, levels, vmin, vmax, cmap in [
+                ('neither', 7, None, None, None),
+                ('neither', 7, None, 20, mpl.cm.RdBu),
+                ('both', 7, 4, 8, None),
+                ('min', 10, 4, 15, None)]:
             for kind in ['imshow', 'pcolormesh', 'contourf', 'contour']:
                 primitive = getattr(self.darray.plot, kind)(
-                    levels=levels, vmin=vmin, vmax=vmax)
+                    levels=levels, vmin=vmin, vmax=vmax, cmap=cmap)
                 assert levels >= \
                     len(primitive.norm.boundaries) - 1
                 if vmax is None:
diff --git a/xarray/tests/test_variable.py b/xarray/tests/test_variable.py
index 4ddd114d..5da83880 100644
--- a/xarray/tests/test_variable.py
+++ b/xarray/tests/test_variable.py
@@ -1540,6 +1540,42 @@ class TestVariable(VariableSubclassobjects):
         assert_identical(
             v.max(), Variable([], pd.Timestamp('2000-01-03')))
 
+    def test_reduce_keepdims(self):
+        v = Variable(['x', 'y'], self.d)
+
+        assert_identical(v.mean(keepdims=True),
+                         Variable(v.dims, np.mean(self.d, keepdims=True)))
+        assert_identical(v.mean(dim='x', keepdims=True),
+                         Variable(v.dims, np.mean(self.d, axis=0,
+                                  keepdims=True)))
+        assert_identical(v.mean(dim='y', keepdims=True),
+                         Variable(v.dims, np.mean(self.d, axis=1,
+                                  keepdims=True)))
+        assert_identical(v.mean(dim=['y', 'x'], keepdims=True),
+                         Variable(v.dims, np.mean(self.d, axis=(1, 0),
+                                  keepdims=True)))
+
+        v = Variable([], 1.0)
+        assert_identical(v.mean(keepdims=True),
+                         Variable([], np.mean(v.data, keepdims=True)))
+
+    @requires_dask
+    def test_reduce_keepdims_dask(self):
+        import dask.array
+        v = Variable(['x', 'y'], self.d).chunk()
+
+        actual = v.mean(keepdims=True)
+        assert isinstance(actual.data, dask.array.Array)
+
+        expected = Variable(v.dims, np.mean(self.d, keepdims=True))
+        assert_identical(actual, expected)
+
+        actual = v.mean(dim='y', keepdims=True)
+        assert isinstance(actual.data, dask.array.Array)
+
+        expected = Variable(v.dims, np.mean(self.d, axis=1, keepdims=True))
+        assert_identical(actual, expected)
+
     def test_reduce_keep_attrs(self):
         _attrs = {'units': 'test', 'long_name': 'testing'}