What is your issue?

I'm trying to perform 2D interpolation on a large 3D array that is heavily chunked along the interpolation dimensions and not the third dimension. The application could be extracting a timeseries from a reanalysis dataset chunked in space but not time, to compare to observed station data with more precise coordinates.

I use the advanced interpolation method as described in the documentation, with the interpolation coordinates specified by DataArray's with a shared dimension like so:

```python %load_ext memory_profiler import numpy as np import dask.array as da import xarray as xr

Synthetic dataset chunked in the two interpolation dimensions

nt = 40000 nx = 200 ny = 200 ds = xr.Dataset( data_vars = { 'foo':( ('t', 'x', 'y'), da.random.random(size=(nt, nx, ny), chunks=(-1, 10, 10))), }, coords = { 't': np.linspace(0, 1, nt), 'x': np.linspace(0, 1, nx), 'y': np.linspace(0, 1, ny), } )

Interpolate to some random 2D locations

ni = 10 xx = xr.DataArray(np.random.random(ni), dims='z', name='x') yy = xr.DataArray(np.random.random(ni), dims='z', name='y') interpolated = ds.foo.interp(x=xx, y=yy) %memit interpolated.compute() ```

With just 10 interpolation points, this example calculation uses about 1.5 * ds.nbytes of memory, and saturates around 2 * ds.nbytes by about 100 interpolation points.

This could definitely work better, as each interpolated point usually only requires a single chunk of the input dataset, and at most 4 if it is right on the corner of a chunk. For example we can instead do it in a loop and get very reasonable memory usage, but this isn't very scalable:

python interpolated = [] for n in range(ni): interpolated.append(ds.foo.interp(x=xx.isel(z=n), y=yy.isel(z=n))) interpolated = xr.concat(interpolated, dim='z') %memit interpolated.compute()

I tried adding a .chunk({'z':1}) to the interpolation coordinates but this doesn't help. We can also do .sel(x=xx, y=yy, method='nearest') with very good performance.

Any tips to make this calculation work better with existing options, or otherwise ways we might improve the interp method to handle this case? Given the performance behavior I'm guessing we may be doing sequntial interpolation for the dimensions, basically an interp1d call for all the xx points and from there another to the yy points, which for even a small number of points would require nearly all chunks to be loaded in. But I haven't explored the code enough yet to understand the details.

What is your issue?

```python import dask.array import xarray as xr

da1 = xr.DataArray(dask.array.ones((3,3), chunks=(1, 1)), dims=["x", "y"]) da2 = xr.DataArray(dask.array.ones((3,), chunks=(1,)), dims=["x"])

da2.broadcast_like(da1).chunksizes

Frozen({'x': (1, 1, 1), 'y': (3,)}) ```

Was surprised to not find any other issues around this. Feels like a major limitation of the method for a lot of use cases. Is there an easy hack around this?

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

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:

Is your feature request related to a problem?

As far as I know, there is no straightforward API in xarray to do what np.sort or pandas.sort_values does. We have DataArray.sortby("x"), which will sort the array according to the coordinate itself. But if instead you want to sort the values of the array to be monotonic, you're on your own. There are probably a lot of ways we could do this, but I ended up with the couple line solution below after a little trial and error.

Describe the solution you'd like

Would there be interest in implementing a Dataset/DataArray.sort_values(dim="x") method?

Note: this 1D example is not really relevant, see the 2D version and more obvious implementation in comments below for what I really want.

python def sort_values(self, dim: str): sort_idx = self.argsort(axis=self.get_axis_num(dim)).drop_vars(dim) return self.isel({dim: sort_idx}).drop_vars(dim).assign_coords({dim: self[dim]})

The goal is to handle arrays that we want to monotize like so:

```python da = xr.DataArray([1, 3, 2, 4], coords={"x": [1, 2, 3, 4]}) da.sort_values("x")

<xarray.DataArray (x: 4)> array([1, 2, 3, 4]) Coordinates: * x (x) int64 1 2 3 4 ```

In addition to sortby which can deal with an array that is just unordered according to the coordinate:

```python da = xr.DataArray([1, 3, 2, 4], coords={"x": [1, 3, 2, 4]}) da.sortby("x")

<xarray.DataArray (x: 4)> array([1, 2, 3, 4]) Coordinates: * x (x) int64 1 2 3 4 ```

Describe alternatives you've considered

I don't know if argsort is dask-enabled (the docs just point to the numpy function). Is there a more intelligent way to implement this with apply_ufunc and something else? I assume chunking in the sort dimension would be problematic.

Additional context

Some past related threads on this topic: https://github.com/pydata/xarray/issues/3957 https://stackoverflow.com/questions/64518239/sorting-dataset-along-axis-with-dask

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

What is your issue?

Follow on from #6891 and #6978 to discuss how we could homogenize the passing of keyword arguments to wrapped external functions across xarray methods.

There are quite a few methods like this where we are ultimately passing data to numpy, scipy, or some other library and want the option to send variable length kwargs to that underlying function. There are two different ways of doing this today:

1. xarray method accepts flexible **kwargs so these can be written directly in the xarray call
2. xarray method accepts a single dict kwargs (sometimes named differently) and passes these along in expanded form via **kwargs

I could only find a few examples of the latter:

• Dataset.interp, which takes kwargs
• Dataset.curvefit, which takes kwargs (although the docstring is wrong here)
• xr.apply_ufunc, which takes kwargs passed to func and dask_gufunc_kwargs passed to dask.array.apply_gufunc
• xr.open_dataset, which takes either **kwargs or backend_kwargs and merges the two

Allowing direct passage with **kwargs seems nice from a user perspective. But, this could occasionally be problematic, for example in the Dataset.interp case where this method also accepts the kwarg form of coords with **coords_kwargs. There are many methods like this that use **indexers_kwargs or **chunks_kwargs with either_dict_or_kwargs but don't happen to wrap external functions.

What is your issue?

I make use of fsspec to quickly open netcdf files in the cloud and pull out slices of data without needing to read the entire file. Quick and dirty is just ds = xr.open_dataset(fs.open("gs://...")).

This works great, in that a many GB file can be lazy-loaded as a dataset in a few hundred milliseconds, by only parsing the netcdf headers with under-the-hood byte range requests. But, only if the netcdf is written from dask-backed arrays. Somehow, writing from numpy-backed arrays produces a different netcdf that requires reading deeper into the file to parse as a dataset.

I spent some time digging into the backends and see xarray is ultimately passing off the store write to dask.array here. A look at ncdump and Dataset.encoding didn't reveal any obvious differences between these files, but there is clearly something. Anyone know why the straight xarray store methods would produce a different netcdf structure, despite the underlying data and encoding being identical?

This should work as an MCVE: ```python import os import string import fsspec import numpy as np import xarray as xr

fs = fsspec.filesystem("gs") bucket = "gs://<your-bucket>"

create a ~160MB dataset with 20 variables

variables = {v: (["x", "y"], np.random.random(size=(1000, 1000))) for v in string.ascii_letters[:20]} ds = xr.Dataset(variables)

Save one version from numpy backed arrays and one from dask backed arrays

ds.compute().to_netcdf("numpy.nc") ds.chunk().to_netcdf("dask.nc")

Copy these to a bucket of your choice

fs.put("numpy.nc", bucket) fs.put("dask.nc", bucket) ```

Then time reading in these files as datasets with fsspec: ```python %timeit xr.open_dataset(fs.open(os.path.join(bucket, "numpy.nc")))

2.15 s ± 40.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

```

```python %timeit xr.open_dataset(fs.open(os.path.join(bucket, "dask.nc")))

187 ms ± 26.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

```

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

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

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: