xarray mean and numpy mean yield very different result since v2025.1.0

[sanghyukmoon]

What happened?

For the particular example given below, DataArray.mean() returned a factor of ~4 smaller value compared to the corresponding numpy mean, DataArray.data.mean() (see the minimal working example below). This behavior is not observed with v2024.11.0, but appeared using v2025.1.0 and later. Also, the problem depends on the size of the array and the data precision (float32 or float64)

Interestingly, DataArray.mean(skipna=False) produced correct result (i.e., identical with numpy mean), even though the DataArray does not contain any NaNs.

What did you expect to happen?

I expect that xarray and numpy mean yield identical results to machine precision.

Minimal Complete Verifiable Example

import numpy as np
import xarray as xr

print(f"xarray version = {xr.__version__}")

N = 512
np.random.seed(0)
arr = np.random.random([N, N, N]).astype(np.float32)
da = xr.DataArray(arr)
print(f"xarray mean = {da.mean().data[()]}")
print(f"numpy mean = {da.data.mean()}")

# output
# ------
# xarray version = 2025.1.0
# xarray mean = 0.125         # <--- wrong value!
# numpy mean = 0.49997711181640625
# ------
# xarray version = 2024.11.0
# xarray mean = 0.49997711181640625
# numpy mean = 0.49997711181640625


# Alternative example with skipna=False
print(f"xarray mean = {da.mean(skipna=False).data[()]}")
print(f"numpy mean = {da.data.mean()}")

# output
# ------
# xarray version = 2025.1.0
# xarray mean = 0.49997711181640625   # Note that with skipna=False, this becomes correct.
# numpy mean = 0.49997711181640625
# ------
# xarray version = 2024.11.0
# xarray mean = 0.49997711181640625
# numpy mean = 0.49997711181640625

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.
• Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

Anything else we need to know?

No response

Environment

INSTALLED VERSIONS

commit: None
python: 3.13.1 | packaged by conda-forge | (main, Jan 13 2025, 09:53:10) [GCC 13.3.0]
python-bits: 64
OS: Linux
OS-release: 4.18.0-553.34.1.el8_10.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: en_US.UTF-8
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.14.4
libnetcdf: 4.9.2

xarray: 2025.1.0
pandas: 2.2.3
numpy: 2.2.2
scipy: 1.15.1
netCDF4: 1.7.2
pydap: None
h5netcdf: None
h5py: 3.12.1
zarr: None
cftime: 1.6.4
nc_time_axis: None
iris: None
bottleneck: 1.4.2
dask: 2025.1.0
distributed: 2025.1.0
matplotlib: 3.10.0
cartopy: None
seaborn: None
numbagg: None
fsspec: 2025.2.0
cupy: None
pint: None
sparse: None
flox: 0.10.0
numpy_groupies: 0.11.2
setuptools: 75.8.0
pip: 25.0
conda: None
pytest: None
mypy: None
IPython: 8.32.0
sphinx: None

[welcome]

Thanks for opening your first issue here at xarray! Be sure to follow the issue template!
If you have an idea for a solution, we would really welcome a Pull Request with proposed changes.
See the Contributing Guide for more.
It may take us a while to respond here, but we really value your contribution. Contributors like you help make xarray better.
Thank you!

[kmuehlbauer]

Thanks @sanghyukmoon for raising this. This is related to #1346. The issue is with bottleneck which has issues with float32 data. Please use something like:

use_bottleneck = True
with xr.set_options(use_bottleneck=use_bottleneck):
    print(f"xarray mean = {da.mean().data[()]}")
    print(f"numpy mean = {da.data.mean()}")

Close this with this comment.