• [x] Closes https://github.com/xarray-contrib/xeofs/issues/148
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

The change from https://github.com/pydata/xarray/pull/8672 mostly fixed the issue of serializing a reset multiindex in the backends, but there was an additional niche issue that turned up in xeofs that was causing serialization to still fail on the zarr backend.

The issue is that zarr is the only backend that uses a custom version of encode_cf_variable called encode_zarr_variable, and the way this gets called we don't pass through the name of the variable before running ensure_not_multiindex.

As a minimal fix, this PR just passes name through as an additional arg to the general encode_variable function. See @benbovy's comment that maybe we should actually unwrap the level coordinate in reset_index and clean up the checks in ensure_not_multiindex, but I wasn't able to get that working easily.

The exact workflow this turned up in involves DataTree and looks like this: ```python import numpy as np import xarray as xr from datatree import DataTree

ND DataArray that gets stacked along a multiindex

da = xr.DataArray(np.ones((3, 3)), coords={"dim1": [1, 2, 3], "dim2": [4, 5, 6]}) da = da.stack(feature=["dim1", "dim2"])

Extract just the stacked coordinates for saving in a dataset

ds = xr.Dataset(data_vars={"feature": da.feature})

Reset the multiindex, which should make things serializable

ds = ds.reset_index("feature") dt1 = DataTree() dt2 = DataTree(name="feature", data=ds) dt1["foo"] = dt2

Somehow in this step, dt1.foo.feature.dim1.variable becomes an IndexVariable again

print(type(dt1.foo.feature.dim1.variable))

Works

dt1.to_netcdf("test.nc", mode="w")

Fails

dt1.to_zarr("test.zarr", mode="w") ```

But we can reproduce in xarray with the test added here.

Is your feature request related to a problem?

If I want to write to a region of data in a zarr, I usually have some boilerplate code like this: python ds_existing = xr.open_zarr(path) ds_new = xr.Dataset(...) # come up with some new data for a subset region of the existing zarr start_idx = (ds_existing.time == ds_new.time[0]).argmax() end_idx = (ds_existing.time == ds_new.time[-1]).argmax() ds_new.to_zarr(path, region={"time": slice(start_idx, end_idx})

Describe the solution you'd like

It would be nice to automate this within to_zarr, because having to drop into index-space always feels un-xarray-like to me.

There may be pitfalls I'm not thinking of, and I don't know exactly what the API would look like. python ds_new.to_zarr(path, region={"time": "auto"}) # ???

Describe alternatives you've considered

No response

Additional context

No response

• [x] Closes #7702, #8421
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

A quick pass at implementing these two improvements for zarr region writes:

1. allow passing region={dim: "auto"}, which opens the existing zarr store and identifies the correct slice to write to, using a variation of the approach suggested by @DahnJ here. We also check for non-matching coordinates and non-contiguous indices.
2. automatically transpose dimensions if they otherwise match the existing store but are out of order

What happened: I've been using xarray to open zarr datasets within a Flask app, and spent some time debugging a memory leak. What I found is that open_zarr() defaults to chunks='auto', rather than chunks=None which is the default for open_dataset(). The result is that open_zarr() ends up calling _maybe_chunk() on the dataset's variables by default.

For whatever reason this function is generating dask items that are not easily cleared from memory within the context of a Flask route, and memory usage continues to grow within my app, at least towards some plateau. This memory growth isn't reproducible outside of a Flask route, so it's a bit of a niche problem.

First proposal would be to simply align the default chunks argument between open_zarr() and open_dataset(). I'm happy to submit a PR there if this makes sense to others. The other more challenging piece would be to figure out whats going on in _maybe_chunk() to cause memory growth. The problem is specific to this function rather than any particular storage backend (other than the difference in default chunk args).

What you expected to happen: Memory usage should not grow when opening a zarr dataset within a Flask route.

Minimal Complete Verifiable Example:

```python from flask import Flask import xarray as xr import gc import dask.array as da

save a test dataset to zarr locally

ds_test = xr.Dataset({"foo": (["x", "y", "z"], da.random.random(size=(300,300,300)))}) ds_test.to_zarr('test.zarr', mode='w')

app = Flask(name)

ping this route repeatedly to see memory increase

@app.route('/open_zarr') def open_zarr(): # with default chunks='auto', memory grows, with chunks=None, memory is ok ds = xr.open_zarr('test.zarr', chunks='auto').compute() # Try to explicity clear memory but this doesn't help del ds gc.collect() return 'check memory'

if name == 'main': app.run(host='0.0.0.0', port=8080, debug=True) ```

Anything else we need to know?:

Environment:

What happened?

There is still something wrong with the groupby indexing changes from 2023.04.0 onward. Here are two typical ways of doing a seasonal anomaly calculation, which return the same result on 2023.03.0 but are different on 2023.04.2:

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

monthly timeseries that should return "zero anomalies" everywhere

time = xr.date_range("2023-01-01", "2023-12-31", freq="MS") data = np.linspace(-1, 1, 12) x = xr.DataArray(data, coords={"time": time}) clim = xr.DataArray(data, coords={"month": np.arange(1, 13, 1)})

seems to give the correct result if we use the full x, but not with a slice

x_slice = x.sel(time=["2023-04-01"])

two typical ways of computing anomalies

anom_gb = x_slice.groupby("time.month") - clim anom_sel = x_slice - clim.sel(month=x_slice.time.dt.month)

passes on 2023.3.0, fails on 2023.4.2

the groupby version is aligning the indexes wrong, giving us something other than 0

assert anom_sel.equals(anom_gb) ```

Related: #7759 #7766 cc @dcherian

What did you expect to happen?

No response

Minimal Complete Verifiable Example

No response

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.

Relevant log output

No response

Anything else we need to know?

No response

Environment

• [x] Closes #7362
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

What happened?

6461 had some unintended consequences for xr.where(..., keep_attrs=True), where coordinate attributes are getting overwritten by variable attributes. I guess this has been broken since 2022.06.0.

What did you expect to happen?

Coordinate attributes should be preserved.

Minimal Complete Verifiable Example

Python import xarray as xr ds = xr.tutorial.load_dataset("air_temperature") xr.where(True, ds.air, ds.air, keep_attrs=True).time.attrs

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.

Relevant log output

```Python

New time attributes are:

{'long_name': '4xDaily Air temperature at sigma level 995', 'units': 'degK', 'precision': 2, 'GRIB_id': 11, 'GRIB_name': 'TMP', 'var_desc': 'Air temperature', 'dataset': 'NMC Reanalysis', 'level_desc': 'Surface', 'statistic': 'Individual Obs', 'parent_stat': 'Other', 'actual_range': array([185.16, 322.1 ], dtype=float32)}

Instead of:

{'standard_name': 'time', 'long_name': 'Time'} ```

Anything else we need to know?

I'm struggling to figure out how the simple lambda change in #6461 brought this about. I tried tracing my way through the various merge functions but there are a lot of layers. Happy to submit a PR if someone has an idea for an obvious fix.

Environment

• [x] Closes #7220
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

Reverts the getattr method used in xr.where(..., keep_attrs=True) from #6461, but keeps handling for scalar inputs. Adds some test cases to ensure consistent attribute handling.

• [x] Closes #6444
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

Fixes a bug introduced by #4687 where passing a non-xarray object to x in xr.where(cond, x, y, keep_attrs=True) caused a failure. The keep_attrs callable passed to merge_attrs() tries to access the attributes of x which do not exist in this case. This fix just checks to make sure x has attributes, and if not will pass through keep_attrs=True.

What happened?

Calling rolling().mean() on a dask-backed array sometimes outputs a different dtype than with a numpy backed array, for example with a float32 input. This is due to the optimized _mean function introduced in #4915.

What did you expect to happen?

This is a simple enough operation that if you start with float32 I would expect to get float32 back.

Minimal Complete Verifiable Example

```Python

import xarray as xr da = xr.DataArray([1,2,3], coords={'x':[1,2,3]}).astype('float32') da.rolling(x=3, min_periods=1).mean().dtype dtype('float32') da.chunk({'x':1}).rolling(x=3, min_periods=1).mean().dtype dtype('float64') ```

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.

Relevant log output

No response

Anything else we need to know?

5877 is somewhat related.

Environment

• [x] Closes #7062
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

This just tests to make sure we at least get the same dtype whether we have a numpy or dask array.

• [x] Closes #6484
• [x] Tests added

I've often been tripped up by the stack trace noted in #6484. This PR changes two things:

1. Handles the zarr GroupNotFoundError error with a more informative FileNotFoundError, displaying the path where we didn't find a zarr store.
2. Moves the consolidated metadata warning to after the step of successfully opening the zarr with non-consolidated metadata. This way the warning isn't shown if we are actually trying to open a non-existent zarr store, in which case we only get the error above and no warning.

What happened: With the latest dask release 2021.4.0, there seems to be a regression in memory management that has broken some of my standard climate science workflows like the simple anomaly calculation below. Rather than intelligently handling chunks for independent time slices, the code below will now load the entire ~30GB x array into memory before writing to zarr.

What you expected to happen: Dask would intelligently manage chunks and not fill up memory. This works fine in 2021.3.0.

Minimal Complete Verifiable Example: Generate a synthetic dataset with time/lat/lon variable and associated climatology stored to disk, then calculate the anomaly: ```python import xarray as xr import pandas as pd import numpy as np import dask.array as da

dates = pd.date_range('1980-01-01', '2019-12-31', freq='D') ds = xr.Dataset( data_vars = { 'x':( ('time', 'lat', 'lon'), da.random.random(size=(dates.size, 360, 720), chunks=(1, -1, -1))), 'clim':( ('dayofyear', 'lat', 'lon'), da.random.random(size=(366, 360, 720), chunks=(1, -1, -1))), }, coords = { 'time': dates, 'dayofyear': np.arange(1, 367, 1), 'lat': np.arange(-90, 90, .5), 'lon': np.arange(-180, 180, .5), } )

My original use case was pulling this data from disk, but it doesn't actually seem to matter

ds.to_zarr('test-data', mode='w') ds = xr.open_zarr('test-data') ds['anom'] = ds.x.groupby('time.dayofyear') - ds.clim ds[['anom']].to_zarr('test-anom', mode='w') ```

Anything else we need to know?: Distributed vs local scheduler and file backend e.g. zarr vs netcdf don't seem to affect this.

Dask graphs look the same for both 2021.3.0: and 2021.4.0:

Environment:

• [x] Closes #4300
• [x] Tests added
• [x] Passes pre-commit run --all-files
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst
• [x] New functions/methods are listed in api.rst

This is a simple implementation of a more general curve-fitting API as discussed in #4300, using the existing scipy curve_fit functionality wrapped with apply_ufunc. It works for arbitrary user-supplied 1D functions that ingest numpy arrays. Formatting and nomenclature of the outputs was largely copied from .polyfit, but could probably be improved.