I see, they are not the same - the slow one is still a dask array, the other one is not: Sn (r, theta, phi, sampling) float64 dask.array<chunksize=(14, 52, 2, 10), meta=np.ndarray>, Sn (r, theta, phi, sampling) float64 nan nan nan nan ... nan nan nan Otherwise they are the same, so this might be dask related ...

I just changed theisel = ds[k].isel(**slc, missing_dims="ignore") to: slcp = [slc[d] if d in slc else slice(None) for d in ds[k].dims] theisel = ds[k].values[tuple(slcp)] And that changed the runtime of my code from (unknown, still running after 3 hours) to around 10 seconds.

ds[k] is a 3 dimensional array slc[d] are 7-d numpy array of integers

The change does matter - but deep copies are still much more expensive than they used to be (as to be expected, I guess)

Indeed, it does help. In 6 hours the CI completed 50% of the tests, compared to 17%.

This however still very much slower than before, where we finished in around half an hour - so around 24x slower ...

A call graph was posted in the referenced thread:

https://github.com/boutproject/xBOUT/pull/252#issuecomment-1282222985 https://user-images.githubusercontent.com/1486942/196415148-ca7ea730-34f6-4622-8f0c-1e98d8b06e26.svg

I tried it with master, and it failed. Trying with main worked :-D

I think there will be cases where the data is not suitable for plotting, where the new error will be less clear than the old one, but I still think that would be overall an improvement.

I think there are several options:

1. remove the offending tests (simple)
2. reintroduce _ensure_plottable but only check for multiindex. This has still the issue that it might at some point be possible, and the tests need to be changed and _ensure_plottable needs to be changed
3. Change the test to check for any error. This basically just means we are checking it isn't working, which I see only limited value in.

Which of these options would you prefer?

I have hacked something that does support the reading and writing of sparse arrays to a netcdf file, however I didn't know how and where to put this within xarray.

``` def ds_to_netcdf(ds, fn): dsorg = ds ds = dsorg.copy() for v in ds: if hasattr(ds[v].data, "nnz") and ( hasattr(ds[v].data, "to_coo") or hasattr(ds[v].data, "linear_loc") ): coord = f"{v}_xarray_index" assert coord not in ds data = ds[v].data if hasattr(data, "to_coo"): data = data.to_coo() ds[coord] = coord, data.linear_loc() dims = ds[v].dims ds[coord].attrs["compress"] = " ".join(dims) at = ds[v].attrs ds[v] = coord, data.data ds[v].attrs = at ds[v].attrs["fill_value"] = str(data.fill_value) for d in dims: if d not in ds: ds[f"_len{d}"] = len(dsorg[d])

print(ds)
ds.to_netcdf(fn)

```

``` def xr_open_dataset(fn): ds = xr.open_dataset(fn)

def fromflat(shape, i):
    index = []
    for fac in shape[::-1]:
        index.append(i % fac)
        i //= fac
    return tuple(index[::-1])

for c in ds.coords:
    if "compress" in ds[c].attrs:
        vs = c.split("_")
        if len(vs) < 5:
            continue
        if vs[-1] != "" or vs[-2] != "index" or vs[-3] != "xarray":
            continue
        v = "_".join(vs[1:-3])
        at = ds[v].attrs
        dat = ds[v].data
        fill = ds[v].attrs.pop("_fill_value", None)
        if fill:
            knownfails = {"nan": np.nan, "False": False, "True": True}
            if fill in knownfails:
                fill = knownfails[fill]
            else:
                fill = np.fromstring(fill, dtype=dat.dtype)
        dims = ds[c].attrs["compress"].split()
        shape = []
        for d in dims:
            try:
                shape.append(len(ds[d]))
            except KeyError:
                shape.append(int(ds[f"_len_{d}"].data))
                ds = ds.drop_vars(f"_len_{d}")

        locs = fromflat(shape, ds[c].data)
        data = sparse.COO(locs, ds[v].data, shape, fill_value=fill)
        ds[v] = dims, data, ds[v].attrs, ds[v].encoding
print(ds)
return ds

```

Has there been any progress since last year?

I have now changed so that several wrapped functions preserve the data.

It is more generic, and hopefully still readable. The flag might also be called identity_0d because the function is the identity function for 0d data, while invariant_0d means that any 0d data is invariant under this function.

Thanks, it is indeed fixed in 070d815 :+1: Should I close #4978 ?

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?

Very useful :+1: I would add: try: c.attrs["units"] = a.attrs["units"] + '*' + b.attrs["units"] except KeyError: pass to preserve units - but I am not sure that is in scope for xarray.