I tried this: diff --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -4701,7 +4701,9 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords): if not reduce_dims: variables[name] = var else: - if ( + if not reduce_dims: + variables[name] = var + elif ( not numeric_only or np.issubdtype(var.dtype, np.number) or (var.dtype == np.bool_) which works great for mean - "var" stays an integer, as expected.

However, that breaks ds.std - which should be zero for "var", but isn't. I guess that is ok for coords - as the assumption is that on coordinates the calculation is not done, but for data variables this is probably not ok.

```diff --- a/xarray/core/duck_array_ops.py +++ b/xarray/core/duck_array_ops.py @@ -537,6 +537,11 @@ def mean(array, axis=None, skipna=None, **kwargs): dtypes""" from .common import _contains_cftime_datetimes

• The mean over an empty axis shouldn't change the data

• See https://github.com/pydata/xarray/issues/4885

• if not axis:
• return array + array = asarray(array) if array.dtype.kind in "Mm": offset = _datetime_nanmin(array) ``` I think it is best to change mean - which would work also for dataArrays. This implies that mean does not convert to float64 - as the numpy version does, but I guess that should be fine.

Should I open a PR?