What happened?

With the 2024.1 release, I am unable to write objects to netCDF after having stacked dimensions with .stack() and called .reset_index() to get rid of the multi-index

What did you expect to happen?

No response

Minimal Complete Verifiable Example

Python import numpy as np import xarray as xr da = xr.DataArray(np.zeros([2, 3]), dims=["x", "y"]) da = da.stack(point=("x", "y")) da = da.reset_index("point") da.to_netcdf("test.nc")

MVCE confirmation

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

Relevant log output

```Python 86 def ensure_not_multiindex(var: Variable, name: T_Name = None) -> None: 87 if isinstance(var._data, indexing.PandasMultiIndexingAdapter): ---> 88 raise NotImplementedError( 89 f"variable {name!r} is a MultiIndex, which cannot yet be " 90 "serialized. Instead, either use reset_index() " 91 "to convert MultiIndex levels into coordinate variables instead " 92 "or use https://cf-xarray.readthedocs.io/en/latest/coding.html." 93 )

NotImplementedError: variable 'x' is a MultiIndex, which cannot yet be serialized. Instead, either use reset_index() to convert MultiIndex levels into coordinate variables instead or use https://cf-xarray.readthedocs.io/en/latest/coding.html. ```

Anything else we need to know?

Creating the stacked object from scratch and saving it to netCDF works fine. The difference is that type(da.x.variable._data) is xarray.core.indexing.PandasMultiIndexingAdapter if it was stacked and reset and numpy.ndarray if it's created from scratch

Environment

Code Sample

```python import numpy as np import xarray as xr

np.random.seed(42)

dims = ('a', 'b', 'c', 'd') shape = (10, 10, 500, 500)

coords = {d: np.arange(s) for d, s in zip(dims, shape)}

Using data with non-normal distribution

data = np.random.lognormal(size=shape) data = data.astype(np.float32)

da = xr.DataArray(data, coords=coords, dims=dims)

Numpy method gives the correct value

print(da.values.mean())

Explicitly specifying all axis gives the correct value

print(da.mean(axis=(0, 1, 2, 3)))

Default DataArray mean method gives incorrect value

print(da.mean()) # <- Problem arise here

float64 arrays produce the correct value

print(da.astype(np.float64).mean()) This is the output I see: 1.6489075

<xarray.DataArray ()> array(1.648908, dtype=float32)

<xarray.DataArray ()> array(1.517693)

<xarray.DataArray ()> array(1.648907) ```

Problem description

Wrong mean value calculated by DataArray.mean() method with default arguments. I have only observed the problem for float32 arrays. It appears to be sensitive to the shape of the array, e.g. a shape of (10, 10, 10, 10) seems to be fine.

Expected Output

This is the output I expect for the sample above: ``` 1.6489075

<xarray.DataArray ()> array(1.648908, dtype=float32)

<xarray.DataArray ()> array(1.648908, dtype=float32)

<xarray.DataArray ()> array(1.648907) ```