What happened?

When calling apply_ufunc with keep_attrs=False, the coordinate attributes are dropped only if there is more than one argument to the call.

What did you expect to happen?

I expected the behaviour to be the same, no matter the number of arguments.

I also expected the coordinate attributes to be preserved if that coordinate was appearing on only one argument.

Minimal Complete Verifiable Example

```Python import xarray as xr

def wrapper(ar1, ar2=None): return ar1.mean(axis=-1)

ds = xr.tutorial.open_dataset("air_temperature")

o1 = xr.apply_ufunc( wrapper, ds.air, ds.time, input_core_dims=[['time'], ['time']], keep_attrs=False ) print(o1.lat.attrs) # {}

o2 = xr.apply_ufunc( wrapper, ds.air, input_core_dims=[['time']], keep_attrs=False ) print(o2.lat.attrs) # {'standard_name': ... } ```

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?

The behaviour stems from this if/else:

https://github.com/pydata/xarray/blob/6d771fc82228bdaf8a4b77d0ceec1cc444ebd090/xarray/core/computation.py#L252-L260

The upper part (1 arg) doesn't touch the attributes, but in the else (more than 1 arg) , two levels deeper in merge_coordinates_without_align , we have:

https://github.com/pydata/xarray/blob/6d771fc82228bdaf8a4b77d0ceec1cc444ebd090/xarray/core/merge.py#L283-L286

When apply_ufunc is called with keep_attrs=False, the combine_attrs above is "drop". In merge_attrs even though there is only one attribute dict passed for lat, it returns an empty dict.

My preference would be for keep_attrs to only refer to the data attributes and that coordinate attributes would be preserved, or even merged if needed. This was my expectation here, as this is the behaviour in many other places of xarray. For example : python3 with xr.set_options(keep_attrs=False): o = ds.air.mean('time') This drops attributes of air, but preserves those of lat and lon.

I see no easy way out here, except by handling it explicitly somewhere in apply_ufunc ? If the decision is that "untouched" coordinate attribute preservation is not ensured by xarray, I think it would be worth noting somewhere (but I don't know where). And I would change my codes to "manually" preserve those where appropriate.

Environment

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

Simply add a encodings argument to save_mfdataset. As for the other args, it expects a list of dictionaries, with encoding information to be passed to to_netcdf for each dataset. Added a minimal test, simply to see if the argument was taken into account.

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

When playing with cftime objects a problem I encountered many times is that I can sub two arrays and them add it back to another. Subtracting to cftime datetime arrays result in an array of np.timedelta64. And when trying to add it back to another cftime array, we get a UFuncTypeError because the two arrays have incompatible dtypes : '<m8[ns]' and 'O'.

Example: ```python import xarray as xr da = xr.DataArray(xr.cftime_range('1900-01-01', freq='D', periods=10), dims=('time',))

An array of timedelta64[ns]

dt = da - da[0]

da[-1] + dt # Fails ```

However, if the two arrays were of 'O' dtype, then the subtraction would be made by cftime which supports datetime.timedelta objects.

This solution here adds a to_pytimedelta to the TimedeltaAccessor, mirroring the name of the similar function on pd.Series.dt. It uses a monkeypatching workaround to prevent xarray to case the array back into numpy objects.

The user still has to check if the data is in cftime or numpy to adapt the operation (calling dt.to_pytimedelta or not), but custom workaround were always overly complicated for such a simple problem, this helps.

Also, this doesn't work with dask arrays because loading a dask array triggers the variable constructor and thus recasts the array of datetime.timedelta to numpy.timedelta[64].

I realize I maybe should have opened an issue before, but I had this idea and it all rushed along.

What happened: When using map_blocks on an object that is dask-backed and has a CFTimeIndex coordinate, the construction step (not computation done) is very slow. I've seen up to 100x slower than an equivalent object with a numpy datetime index.

What you expected to happen: I would understand a performance difference since numpy/pandas objects are usually more optimized than cftime/xarray objects, but the difference is quite large here.

Minimal Complete Verifiable Example:

Here is a MCVE that I ran in a jupyter notebook. Performance is basically measured by execution time (wall time). I included the current workaround I have for my usecase.

```python import numpy as np import pandas as pd import xarray as xr import dask.array as da from dask.distributed import Client

c = Client(n_workers=1, threads_per_worker=8)

Test Data

Nt = 10_000 Nx = Ny = 100 chks = (Nt, 10, 10)

A = xr.DataArray( da.zeros((Nt, Ny, Nx), chunks=chks), dims=('time', 'y', 'x'), coords={'time': pd.date_range('1900-01-01', freq='D', periods=Nt), 'x': np.arange(Nx), 'y': np.arange(Ny) }, name='data' )

Copy of a, but with a cftime coordinate

B = A.copy() B['time'] = xr.cftime_range('1900-01-01', freq='D', periods=Nt, calendar='noleap')

A dumb function to apply

def func(data): return data + data

Test 1 : numpy-backed time coordinate

%time outA = A.map_blocks(func, template=A) # %time outA.load();

Res on my machine:

CPU times: user 130 ms, sys: 6.87 ms, total: 136 ms

Wall time: 127 ms

CPU times: user 3.01 s, sys: 8.09 s, total: 11.1 s

Wall time: 13.4 s

Test 2 : cftime-backed time coordinate

%time outB = B.map_blocks(func, template=B) %time outB.load();

Res on my machine

CPU times: user 4.42 s, sys: 219 ms, total: 4.64 s

Wall time: 4.48 s

CPU times: user 13.2 s, sys: 3.1 s, total: 16.3 s

Wall time: 26 s

Workaround in my code

def func_cf(data): data['time'] = xr.decode_cf(data.coords.to_dataset()).time return data + data

def map_blocks_cf(func, data): data2 = data.copy() data2['time'] = xr.conventions.encode_cf_variable(data.time) return data2.map_blocks(func, template=data)

Test 3 : cftime time coordinate with encoding-decoding

%time outB2 = map_blocks_cf(func_cf, B) %time outB2.load();

Res

CPU times: user 536 ms, sys: 10.5 ms, total: 546 ms

Wall time: 528 ms

CPU times: user 9.57 s, sys: 2.23 s, total: 11.8 s

Wall time: 21.7 s

```

Anything else we need to know?: After exploration I found 2 culprits for this slowness. I used %%prun to profile the construction phase of map_blocks and found that in the second case (cftime time coordinate):

1. In map_blocks calls to dask.base.tokenize take the most time. Precisely, tokenizing a numpy ndarray of O dtype goes through the pickling process of the array. This is already quite slow and cftime objects take even more time to pickle. See Unidata/cftime#253 for the corresponding issue. Most of the construction phase execution time is spent pickling the same datetime array at least once per chunk.
2. Second, but only significant when the time coordinate is very large (55000 in my use case). CFTimeIndex.__new__ is called more than twice as many times as there are chunks. And within the object creation there is this line : https://github.com/pydata/xarray/blob/3956b73a7792f41e4410349f2c40b9a9a80decd2/xarray/coding/cftimeindex.py#L228 The larger the array, the more time is spent in this iteration. Changing the example above to use Nt = 50_000, the code spent a total of 25 s in dask.base.tokenize calls and 5 s in CFTimeIndex.__new__ calls.

My workaround is not the best, but it was easy to code without touching xarray. The encoding of the time coordinate changes it to an integer array, which is super fast to tokenize. And the speed up of the construction phase is because there is only one call to encode_cf_variable compared to N_chunks calls to the pickling, As shown above, I have not seen a slowdown in the computation phase. I think this is mostly because the added decode_cf calls are done in parallel, but there might be other reason I do not understand.

I do not know for sure how/why this tokenization works, but I guess the best improvment in xarray could be to: - Look into the inputs of map_blocks and spot cftime-backed coordinates - Convert those coordinates to a ndarray of a basic dtype. - At the moment of tokenization of the time coordinates, do a switheroo and pass the converted arrays instead.

I have no idea if that would work, but if it does that would be the best speed-up I think.

Environment:

What happened: Given a dataset that has a variable with cftime objects along dimension A, averaging (mean) leads to buggy behaviour:

1. Averaging over 'A' drops the variable instead of averaging it.
2. Averaging over any other dimension will fail if that variable is on the dask backend.

What you expected to happen:

1. I expected the average to fail in the case of a dask-backed cftime variable, given that this code exists: https://github.com/pydata/xarray/blob/fdabf3bea5c750939a4a2ae60f80ed34a6aebd58/xarray/core/duck_array_ops.py#L562-L572

And I expected the average to work (not drop the var) in the case of the numpy backend.

1. I expected the fact that dask is used to be irrelevant to the result. I expected the mean to conserve the cftime variable as-is since it doesn't include the averaged dimension.

Minimal Complete Verifiable Example:

```python

Put your MCVE code here

import xarray as xr

ds = xr.Dataset({ 'var1': (('time',), xr.cftime_range('2021-10-31', periods=10, freq='D')), 'var2': (('x',), list(range(10))) })

var1 contains cftime objects

var2 contains integers

They do not share dims

ds.mean('time') # var1 has disappeared instead of being averaged

ds.mean('x') # Everything ok

dsc = ds.chunk({})

dsc.mean('time') # var1 has disappeared. I would expected this line to fail.

dsc.mean('x') # Raises NotImplementedError. I would expect this line to run flawlessly. ```

Anything else we need to know?: A culprit is #5393, but maybe the bug is older? I think the change introduced there causes the issue (2) above.

In duck_array_ops.py the mean operation is declared numeric_only, which is kinda incoherent with the implementation allowing means of datetime objects. This setting causes my (1) above.

Environment:

What happened:

Quite a specific bug! Using map_blocks to wrap a polyfit computation, using a dask client (not the local scheduler) and a time axis with a cftime calendar, I got the error : TypeError: cannot compute the time difference between dates with different calendars.

What you expected to happen:

No bug.

Minimal Complete Verifiable Example:

```python ds = xr.tutorial.open_dataset('rasm').chunk({'x': 25, 'y': 25})

templ = ds.Tair

def func(ds, verbose=False): # Dummy function that call get_clean_interp_index function # Return the Tair as-is just for the test. if verbose: print(ds.time) print(type(ds.time[0].item())) x = xr.core.missing.get_clean_interp_index(ds, 'time')

return ds.Tair

This works (time is a coordinate, so it is already loaded

x = xr.core.missing.get_clean_interp_index(ds, 'time')

This works too. The local scheduler is used.

out = ds.map_blocks( func, template=templ, kwargs={'verbose': False} ) out.load()

This fails!

with Client(n_workers=1, threads_per_worker=8, dashboard_address=8786, memory_limit='7GB') as c: out = ds.map_blocks( func, template=templ, kwargs={'verbose': True} ) out.load() ```

The full traceback is here:

The printout to the console. I am calling this in a jupyter notebook so the prints from within workers are in the console, not in the cell's output. I removed useless lines to shorten it.

``` <xarray.DataArray 'time' (time: 36)> array([cftime.datetime(1980, 9, 16, 12, 0, 0, 0, calendar='noleap', has_year_zero=True), cftime.datetime(1980, 10, 17, 0, 0, 0, 0, calendar='noleap', has_year_zero=True), cftime.datetime(1980, 11, 16, 12, 0, 0, 0, calendar='noleap', has_year_zero=True), .... )], dtype=object) Coordinates: * time (time) object 1980-09-16 12:00:00 ... 1983-08-17 00:00:00 Attributes: long_name: time type_preferred: int

<class 'cftime._cftime.datetime'> ```

And for reference: ```python

ds.time array([cftime.DatetimeNoLeap(1980, 9, 16, 12, 0, 0, 0, has_year_zero=True), cftime.DatetimeNoLeap(1980, 10, 17, 0, 0, 0, 0, has_year_zero=True), cftime.DatetimeNoLeap(1980, 11, 16, 12, 0, 0, 0, has_year_zero=True), ... type(ds.time[0].item()) cftime._cftime.DatetimeNoLeap ``` Anything else we need to know?:

I'm not sure where the exact breaking change lies (dask or cftime?), but this worked with dask 2021.5 and cftime <= 1.4.1. The problem lies in get_clean_interp_index, specifically these lines:

https://github.com/pydata/xarray/blob/5ccb06951cecd59b890c1457e36ee3c2030a67aa/xarray/core/missing.py#L274-L280

On the original dataset, the class of the time values is DatetimeNoLeap whereas the time coordinates received by func are of class datetime, the calendar is only a kwargs. Thus, in get_clean_interp_index the offset is created with the default "standard" calendar and becomes incompatible with the array itself. Which makes datetime_to_numeric fail.

Environment:

What happened: I have a computation that output dates but that sometimes also outputs missing data. (It computes the start date of a run in a timeseries, if there is no run, it outputs NaN). Afterwards, I'd like to convert those dates to dayofyear, thus I call out.dt.dayofyear. In a case where the first value of out is missing, it fails.

What you expected to happen: I expected out.dt.dayofyear to return an array where the first value would be NaN.

Minimal Complete Verifiable Example:

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

da = xr.DataArray( [[np.nan, np.nan], [1, 2]], dims=('x', 'time'), coords={'x': [1, 2], 'time': xr.cftime_range('2000-01-01', periods=2)}, )

out is a "object" array, where the first element is NaN

out = da.idxmin('time')

out.dt.dayofyear

Expected : [nan, 1.]

Got:

TypeError Traceback (most recent call last) <ipython-input-56-06aa9bdfd6b8> in <module> ----> 1 da.idxmin('time').dt.dayofyear

~/.conda/envs/xclim-dev/lib/python3.8/site-packages/xarray/core/utils.py in get(self, obj, cls) 917 return self._accessor 918 --> 919 return self._accessor(obj) 920 921

~/.conda/envs/xclim-dev/lib/python3.8/site-packages/xarray/core/accessor_dt.py in new(cls, obj) 514 # do all the validation here. 515 if not _contains_datetime_like_objects(obj): --> 516 raise TypeError( 517 "'.dt' accessor only available for " 518 "DataArray with datetime64 timedelta64 dtype or "

TypeError: '.dt' accessor only available for DataArray with datetime64 timedelta64 dtype or for arrays containing cftime datetime objects. ```

Anything else we need to know?: This also triggers computation when da is lazy. A lazy .dt accessor would be useful.

The laziness of it aside, would it be meaningful to change: https://github.com/pydata/xarray/blob/d4b7a608bab0e7c140937b0b59ca45115d205145/xarray/core/common.py#L1822 to cycle on the array while np.isnan(sample) ?

Environment:

What happened: With cftime 1.2.0, one can create datetime object with cftime.datetime(*args, calendar='calendar'), instead of using one of the subclasses (ex cftime.DatetimeNoLeap(*args)). In the latest release (1.4.0, yesterday), the subclasses have been deprecated, but kept as legacy. While all xr code still works (it is using the legacy subclasses), the CFTimeIndex object relies on the type of the datetime object in order to infer the calendar. If the datetime was created outside xarray, using the now default constructor, the returned type is not understood and CFTimeIndexbreaks.

What you expected to happen: I expected CFTimeIndex to be independent of the way the datetime object is created.

Minimal Complete Verifiable Example:

```python3 import cftime import numpy as np import xarray as xr

A datetime array, not created in xarray

time = cftime.num2date(np.arange(365), "days since 2000-01-01", calendar="noleap") a = xr.DataArray(np.zeros(365), dims=('time',), coords={'time': time})

a.indexes['time'] Fails with :python3 Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/home/phobos/Python/xclim/.tox/py38/lib/python3.8/site-packages/xarray/coding/cftimeindex.py", line 342, in repr attrs_str = format_attrs(self) File "/home/phobos/Python/xclim/.tox/py38/lib/python3.8/site-packages/xarray/coding/cftimeindex.py", line 264, in format_attrs attrs["freq"] = f"'{index.freq}'" if len(index) >= 3 else None File "/home/phobos/Python/xclim/.tox/py38/lib/python3.8/site-packages/xarray/coding/cftimeindex.py", line 692, in freq return infer_freq(self) File "/home/phobos/Python/xclim/.tox/py38/lib/python3.8/site-packages/xarray/coding/frequencies.py", line 96, in infer_freq inferer = _CFTimeFrequencyInferer(index) File "/home/phobos/Python/xclim/.tox/py38/lib/python3.8/site-packages/xarray/coding/frequencies.py", line 105, in init self.values = index.asi8 File "/home/phobos/Python/xclim/.tox/py38/lib/python3.8/site-packages/xarray/coding/cftimeindex.py", line 673, in asi8 [ File "/home/phobos/Python/xclim/.tox/py38/lib/python3.8/site-packages/xarray/coding/cftimeindex.py", line 674, in <listcomp> _total_microseconds(exact_cftime_datetime_difference(epoch, date)) File "/home/phobos/Python/xclim/.tox/py38/lib/python3.8/site-packages/xarray/core/resample_cftime.py", line 370, in exact_cftime_datetime_difference seconds = b.replace(microsecond=0) - a.replace(microsecond=0) File "src/cftime/_cftime.pyx", line 1153, in cftime._cftime.datetime.sub ValueError: cannot compute the time difference between dates with different calendars ``` Anything else we need to know?:

Environment:

What happened: When trying to interpolate a N-D array with 2 other arrays sharing a common (new) dimension and with one (at least) being multidimensional fails. Kinda a complex edge case I agree. Here's a MWE: python3 da = xr.DataArray([[[1, 2, 3], [2, 3, 4]], [[1, 2, 3], [2, 3, 4]]], dims=('t', 'x', 'y'), coords={'x': [1, 2], 'y': [1, 2, 3], 't': [10, 12]}) dy = xr.DataArray([1.5, 2.5], dims=('u',), coords={'u': [45, 55]}) dx = xr.DataArray([[1.5, 1.5], [1.5, 1.5]], dims=('t', 'u'), coords={'u': [45, 55], 't': [10, 12]}) So we have da a 3D array with dims (t, x, y). We have dy, containing the values of y along new dimension u. And dx containing the values of x along both u and t. We want to interpolate with: python3 out = da.interp(y=dy, x=dx, method='linear') As so to have a new array over dims t and u.

What you expected to happen: I expected (with the dummy data I gave): python3 xr.DataArray([[2, 3], [2, 3]], dims=('t', 'u'), coords={'u': [45, 55], 't': [10, 12]})

But instead it fails with ValueError: axes don't match array.

Full traceback:

Anything else we need to know?: It works if dx doesn't vary along t. I .e.: da.interp(y=dy, x=dx.isel(t=0, drop=True), method='linear') works.

Environment: