This makes complete sense to me, thanks! I'm quite surprised I didn't get this issue before in my scripts, sorry about that!

This is happenning to me too. For the record, I need to use engine='pydap' because it supports authentification to access private thredds data programmatically (by passing session), while netCDF4 does not.

I pointed to this issue in the related (but currently badly named and stale) issue upstream.

Hmm, I don't think so! Thanks for the heads up.

I found another problem that this issue causes : the dataset can't be written to netCDF and the error message suggests exactly what has already been done: ```python import xarray as xr ds = xr.tutorial.open_dataset('air_temperature')

ds = ds.stack(spatial=['lon', 'lat']) ds = ds.reset_index('spatial', drop=True)

ds.to_netcdf('test.nc') raises: NotImplementedError: variable 'lat' is a MultiIndex, which cannot yet be serialized to netCDF files (https://github.com/pydata/xarray/issues/1077). Use reset_index() to convert MultiIndex levels into coordinate variables instead. ```

@benbovy I'm guessing your PR will fix this, but I'm raising the issue to be sure.

I have the same problem in xarray 2022.3.0. The issue is that this creates unnecessary dask tasks in the graph and some operations acting on the coordinates unexpectedly trigger dask computations. "Unexpected" because the coordinates at the beginning of the process where not chunked. So computation that was expected to happen in the main thread (or not happen at all) is now happenning in the dask workers.

An example: ```python3 import numpy as np import xarray as xr from dask.diagnostics import ProgressBar

A 2D variable

da = xr.DataArray( np.ones((12, 10)), dims=('x', 'y'), coords={'x': np.arange(12), 'y': np.arange(10)} )

A 1D variable sharing a dim with da

db = xr.DataArray( np.ones((12,)), dims=('x'), coords={'x': np.arange(12)} )

A non-dimension coordinate

cx = xr.DataArray(np.zeros((12,)), dims=('x',), coords={'x': np.arange(12)})

Assign it to da and db

da = da.assign_coords(cx=cx) db = db.assign_coords(cx=cx)

We need to chunk along y

da = da.chunk({'y': 1})

Notice how cx is now a dask array, even if it is a 1D coordinate and does not have 'Y' as a dimension.

print(da)

This triggers a dask computation

with ProgressBar(): da - db ``` The reason my example triggers dask is that xarray ensure the coordinates are aligned and equal (I think?). Anyway, I didn't expect it.

Personally, I think the chunk method shouldn't apply to the coordinates at all, no matter their dimensions. They're coordinate so we expect to be able to read them easily when aligning/comparing dataset. Dask is to be used with the "real" data only. Does this vision fit the one from the devs? I feel this "skip" could be easily implemented.

To add details, I suspect this commit as the culprit: https://github.com/headtr1ck/xarray/commit/cd2f2b247cf6e1c68c7db6e92100898272fb0e1a

It removes a direct import of the core.resample submodule in core.dataarray and moves it to the "Type checking" imports. Thus, importing xarray (import xarray as xr), does not trigger an explicit import of resample anymore.

Presumably you haven't found any bugs?

Not for now,, but I haven't done much and it wasn't very complex.

You can pass method as .mean(..., method=...)

Yes! I saw that, but I did search a bit. In fact, my idea was triggered partly because I'm lazy. My suggestion was to do something alike what xarray does with keep_attrs : a global/context option that changes the default value of method, so it doesn't need to be explicitly added everywhere.

Note you can only blockwise if all members of a group are in a single block.

Oh! I didn't catch that. My test data was already appropriately chunked.

Right now, I'd just like this to get merged :)

Haha, I understand. My idea, if pertinent, could be done in another PR.

Hi @dcherian and @andersy005, thanks a lot for this PR and for flox. I'm currently trying this with our climate indices library (xclim) and it works well! I get a clear improvement in performance,

In xclim's indicators, all groupby operations are time resampling ops and the best performance is achieved with method='blockwise', by far. Is there a plan to let users control the default method through something like xr.set_options ?

Like, use_flox could accept False, True or any method name? (The default still being map-reduce) ? What do you think?

I agree that it kinda crowds the xarray API. My first idea was to implement these calendar utilities in cf-xarray, but I was redirected here. Seems to me that if that move is made, it's the most logical place for all CF[time] related things.

Here we are! Also, I don't think the test failure comes from my changes.

@spencerkclark Done!

Note that we have this issue : Ouranosinc/xclim#841 open on our side. It's for a third method for converting 360_day calendars. Based on something used in the litterature, it randomizes the days to insert/drop when converting. My plan is to implement it in xclim soonish, then if it works I could bring it up here and open a new PR.

@spencerkclark After more than 2 weeks, I had time to address your comments!

@spencerkclark I came back to this today, merged upstream/main and now I can't commit the merge because xarray/core/_typed_ops.pyi is not passing the mypy pre-commit hook. Is this a problem on my side? Have you ever seen this? It's strange cause in setup.cfg only *.py files are supposed to be checked by mypy, but I guess pre-commit doesn't respect that? pre-commit run --all does work though...

May be @max-sixty has an idea, since he recently pushed a PR playing with this file and typing...

The errors :

@aulemahal what do you think about how easy it might be to change polyfit to use apply_ufunc?

Should be easy! The code was written without apply_ufunc for a few reasons, one being that using apply_ufunc+vectorize was slower than the apply_along_axis trick it currently uses. But, with the issue raised here, that might not be true...

EDIT: I hit enter a bit too soon... I think I found the very source of the problem, maybe this could be fixed in cftime:

```python import pickle

ds.time # same No-Leap coordinate as in the MWE above

Items are instances of cftime.DatetimeNoLeap

pickle.loads(pickle.dumps(ds.time))

Items are instances of cftime.datetime, calendar information is now only an attribute, get_clean_interp_index will fail

```

However, as xarray has the possiblity to contain cftime coordinates where the calendar is an attribute, I guess get_clean_interp_index needs to be fixed?

Agreed that this is limited and there's need for a more general solution!

However, this problem is quite tricky... The best I see is what I think you suggested in the other thread : special workaround with re-casting when the operation involves "<m8[ns]" and "O" (and that the "O" is cftime datetimes). That clears the instability problem of my solution. It loads a value, so dask-compat is not optimal, but that's already what dt does.

The encoding problem remains, but I made some tests and I realized numpy automatically casts timedeltas to the smallest unit, but accepts operations between timedeltas of different units. Could we accept timedelta64 arrays with another units than "ns"? Then, there could be a check when converting python timedeltas : if nanosecond information exist : go with [ns], if not, go with [ms]. This bumps the limit from 292 years to 292471. And, with this relaxation of the dtype strictness, would allow a user to use even coarser frequencies if need be.

Throwing ideas. There might be other issues that I did not see!

FYI, I opened SciTools/nc-time-axis#59 to update nc-time-axis. Waiting for maintainers to review.

I'm not sure I understand the "dimensions" thing, isn't the MWE 1D? What is the case where it should work?

Thanks everyone. Just to be clear, I am using map_blocks in a rather complicated way (way more than this MWE) and at the moment of the map_blocks call, the code doesn't know the dimension order of the output, only the names. Moreover, in this code, the wrapped functions are too complex for map_blocks to infer the template. The idea is that all wrapped funcs reduce some dims or add some other (but with known sizes), thus I can construct a template, but the dim order is uncertain.

However, I understand that it is expected from map_blocks that the result fits the template, so my expectation was not correct. I'll leave the "workaround" in my code and close this issue. Thanks!

Haha, sorry, I didn't explain this properly.

This MWE is just MWE. The use case I have of map_blocks is in a decorator. At the moment I write the map_blocks call, I do not know if func will transpose the result or not. In my case, the template is generated from a list of dimensions to reduce and new dimensions.

I realize this might be an edge case that is too narrow to merit a fix. I expected dimension order not to matter, as it is the case in most of xarray. Right now, my code wraps func in a another function that performs this tranpose when needed.

Hi @max-sixty, just edited the top comment to include the traceback.

The MWE does work with return d.

@spencerkclark Applied suggestions from the review, added Dataset.convert_calendar and Dataset.interp_calendar, and same on DataArray.

Wasn't sure what to do with the tests. They already done in test_calendar_ops.py...

Updated the code following first wave of suggestions!

I see the strange errors here : https://github.com/aulemahal/xarray/actions. The failure triggered an email and this is why I didn't look here. Well, if it's something in the configuration of my fork, then we don't care!

Will fix the doc errors...

I am excessively confused at the errors in the builds... They are triggered by tests not touching my work and also it says the installed version of xarray is 0.11???

Edited the comment above as I realized that numpy /pandas / xarray use nanoseconds by default, so the earliest date possible is 1677-09-21 00:12:43.145225. Thus, I suggest we return "standard" as the calendar of numpy-backed datetime indexes.

Cool! I started a branch, will push a PR soon.

I understand the "default" issue and using use_cftime=None makes sense to me!

~For dt.calendar and date_range, there remains the question on how we name numpy's calendar: Python uses what CF conventions call proleptic_gregorian, but the default and most common calendar we see and use is CF's "standard". May be users would expect "standard"? A solution would be to check if the earliest value in the array is before 1582-10-15. If yes, return "proleptic_gregorian", if not, return "standard".~

Took a look and it seems to originate from the stacking part and someting in dask.

In polyfit, we rearrange the DataArrays to 2D arrays, so we can run the least squares with np/dsa.apply_along_axis. But I checked and the chunking problem seems to appear before any call of the sort. MWE:

```python3 import xarray as xr import dask.array as dsa nz, ny, nx = (10, 20, 30) data = dsa.ones((nz, ny, nx), chunks=(1, 5, nx)) da = xr.DataArray(data, dims=['z', 'y', 'x']) da.chunks

((1, 1, 1, 1, 1, 1, 1, 1, 1, 1), (5, 5, 5, 5), (30,))

stk = da.stack(zy=['z', 'y']) print(stk.dims, stk.chunks)

('x', 'zy') ((30,), (20, 20, 20, 20, 20, 20, 20, 20, 20, 20))

Merged chunks!

```

And then I went down the rabbit hole (ok it's not that deep) and is all goes down here: https://github.com/pydata/xarray/blob/e0358e586079c12525ce60c4a51b591dc280713b/xarray/core/variable.py#L1507

In Variable._stack_one the stacking is performed and Variable.data.reshape is called. Dask itself is rechunking the output, merging the chunks. There is a merge_chunks kwarg for reshape, but I think it has a bug:

```python

Let's stack as xarray does: x, z, y -> x, zy

data_t = data.transpose(2, 0, 1) # Dask array with shape (30, 10, 20), the same as reordered in Variable._stack_once.

new_data = data_t.reshape((30, -1), merge_chunks=True) # True is the default, this is the same call as in xarray new_data.chunks

((30,), (20, 20, 20, 20, 20, 20, 20, 20, 20, 20))

new_data = data_t.reshape((30, -1), merge_chunks=False) new_data.shape # I'm printing shape because chunks is too large, but see the bug:

(30, 6000) # instead of (30, 200)!!!

Doesn't happen when we do not transpose. So let's reshape data as z, y, x -> zy, x

new_data = data.reshape((-1, 30), merge_chunks=True) new_data.chunks

((5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5), (30,))

Chunks were not merged? But this is the output expected by paigem.

new_data = data.reshape((-1, 30), merge_chunks=False) new_data.chunks

((5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5), (30,))

That's what I expected with merge_chunks=False.

```

For polyfit itself, the apply_along_axis call could be changed to a apply_ufunc with vectorize=True, I think this would avoid the problem and behave the same on the user's side. Would need some refactoring.

@aaronspring Oh, I didn't think of that trick for 1, thanks! But again, this fails with numpy-backed time coordinates. With the definition of a "default" calendar, there could be a more general way.

Dang. How did I miss that?

Woups, the bug seems more important than I first thought. It happens when reading files: MWE: ```python import xarray as xr

da = xr.DataArray( [0, 1, 2], dims=('time',), coords={'time': xr.cftime_range('2000-01-01', freq='D', periods=3, calendar='noleap')}, name='test' ) da.to_netcdf('test.nc')

da2 = xr.open_dataset('test.nc') da2.indexes['time'] Fails with an error similar to the example above. Also failing on :python da2.sel(time='2000') ```

Oh! I should have thought of this. The error message is indeed unclear... I'll close the issue as the problem is solved, but I'll note that: python3 da.sel(y=dy, x=dx, method='nearest') does work, without the need to explicitly pass the shared dimension t. This is why I expected that interp would work the same.

It doesn't work with 0.16.0, but does with later versions (just tried with the current master). Thanks!

I did encounter this when writing up polyfit! I decided to discard that for the time being, but I believe changing get_clean_interp_index would the best way to fix this. One would also need to modify polyfit sligthly, though. Here: https://github.com/pydata/xarray/blob/9c85dd5f792805bea319f01f08ee51b83bde0f3b/xarray/core/dataset.py#L6064-L6075 When the non-fitted dimensions are stacked together, dim should be renamed to the real dimension of the variable. I will not have time to make a PR, but it should be quite simple, if anyone has time!

@mathause Fixed the tests to catch the warning and added a case for the specific of the issue : too many nans causing a rank deficiency.

So while adding tests, I realized there were more bugs concerning deficient rank matrices and full=True output. I believe I fixed all I could find. However, I was not able to reduce the number of warnings for the case using dask : the computation occurs outside of the catch_warnings() context... Also, there is a bug in the dask.array.linalg.lstsq output (see issue dask/dask#6516), so if full=True we must use skipna=True with dask. Same when rank != order (deficient rank), dask.array.linalg.lstsq will fail, so we work around by using our (slower) nan-skipping method.

My comments Q.1 : For now xr.apply_ufunc does not accept core dimensions to be chunked along, which would be kinda a sad for curve fitting. However, dask's least square method does, which is part of the reason why I used it in polyfit. On the other hand, scipy's least-squares procedure is not simple. Curve fitting is quite complex and rewriting all the code to use dask might be a project too ambitious, and surely out-of-scope for xarray...

Q.3 : For simple directly declared function, inspect does a good job, but it can get tricky with wrapped functions, which might arise in more complex workflows. Could we have a params arg that takes in a list of names?

@spencerkclark Here you go! Tests added for all cases you mentioned.

@spencerkclark Thanks for the feedback!

I rewrote the _CFTimeFrequencyInferer class independently of pandas private objects. I removed outputs that wouldn't match offsets defined in xarray: weekly frequencies, week-of-month (WOM) and all business-related anchors. It wouldn't be difficult to add them later on if need be and if cftime_offsets.py defines them. Note that infer_freq is quite independent of the rest of xarray, this choice was purely so that the behavior matches the rest of xr. I did leave the sub-seconds frequencies (L for ms and U for µs).

I expect the performance of xr.infer_freq to be slightly worst than pandas' since they use optimized and compiled code whereas I only used common numpy and xarray stuff. However, I assume that this part of the workflow is far from being the bottleneck of most xarray usecases...

Will write up some docs soon.

Yay! Many thanks, this will greatly improve our projects.

Ping. It be green.

I pushed a new PR trying to implement polyfit in xarray, #3733. It is still work in progress, but I would like the opinion on those who participated in this thread.

Considering all options discussed in the thread, I chose an implementation that seemed to give the best performance and generality (skipping NaN values), but it duplicates a lot of code from numpy.polyfit.

Main question:

• Should xarray's implementation really replicate the behaviour of numpy's?

A lot of extra code could be removed if we'd say we only want to compute and return the residuals and the coefficients. All the other variables are a few lines of code away for the user that really wants them, and they don't need the power of xarray and dask anyway.

I'm guessing @huard @dcherian @rabernat and @shoyer might have comments.

Oh dask="allowed", that's true! For the coeficients and residuals, I can easily do this stacking and unstacking and return only one variable. So, apply_ufunc could be used. However, it would not be very useful right now, considering they way I implemented the least-squares fit directly. Meaning: the left hand side computation and scaling and the not-fitted-dimensions stacking would be done on each call of numpy.polyfit if I did call it directly. Right now, we save some of this overhead by doing it once before iterating on 1d axes (which is anyway needed with skipna=True).