One other tricky edge case that occurs to me is one where an extreme fill value (e.g. 1e30) is used for floating point fields. If we decode the times first, it might appear that the dates cannot be represented as nanosecond-precision values, but in reality they would be. We may need to think more about how to handle this edge case in addition to #7817.

Great, yeah, that's a nice example without writing to disk. Indeed I saw those warnings too, but omitted them in my earlier message to focus on the encoding issue (sorry about that). I agree that these are something we should address.

For a little more discussion of this issue see #7493. As the example at the top of the issue notes, there is at least one place where non-nanosecond times can slip in (this is not intentional), but for most code pathways xarray should currently convert things automatically.

Xarray will convert any non-nanosecond precision times to nanosecond precision (or an error will be raised if this is not possible).

Thanks @kmuehlbauer -- I just wanted to give you a heads up that I'm pretty busy this week. Hopefully I'll get a free moment to look at this more closely next week.

For the time being xarray should behave as it always has, converting any non-nanosecond precision times to nanosecond-precision before being used internally.

Thanks for the ping @dcherian -- I just gave #7098 a review. I think it's close to ready to merge.

The cftime.Datetime360Day objects are still expected here. We have not switched over to the universal cftime.datetime class yet within xarray, though my sense is this would likely still be an issue regardless (the traceback indicates the error comes up in xarray's indexing logic). cftime.Datetime360Day instances are instances of cftime.datetime objects:

```

import cftime isinstance(cftime.Datetime360Day(2000, 1, 1), cftime.datetime) True ```

This is the full traceback for reference: Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/Users/spencer/software/xarray/xarray/core/dataset.py", line 1917, in to_netcdf return to_netcdf( # type: ignore # mypy cannot resolve the overloads:( ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/software/xarray/xarray/backends/api.py", line 1216, in to_netcdf dump_to_store( File "/Users/spencer/software/xarray/xarray/backends/api.py", line 1263, in dump_to_store store.store(variables, attrs, check_encoding, writer, unlimited_dims=unlimited_dims) File "/Users/spencer/software/xarray/xarray/backends/common.py", line 269, in store variables, attributes = self.encode(variables, attributes) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/software/xarray/xarray/backends/common.py", line 358, in encode variables, attributes = cf_encoder(variables, attributes) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/software/xarray/xarray/conventions.py", line 773, in cf_encoder _update_bounds_encoding(variables) File "/Users/spencer/software/xarray/xarray/conventions.py", line 347, in _update_bounds_encoding ) or contains_cftime_datetimes(v) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/software/xarray/xarray/core/common.py", line 1818, in contains_cftime_datetimes return _contains_cftime_datetimes(var._data) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/software/xarray/xarray/core/common.py", line 1811, in _contains_cftime_datetimes return isinstance(np.asarray(sample).item(), cftime.datetime) ^^^^^^^^^^^^^^^^^^ File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/array/core.py", line 1700, in __array__ x = self.compute() ^^^^^^^^^^^^^^ File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/base.py", line 314, in compute (result,) = compute(self, traverse=False, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/base.py", line 599, in compute results = schedule(dsk, keys, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/threaded.py", line 89, in get results = get_async( ^^^^^^^^^^ File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/local.py", line 511, in get_async raise_exception(exc, tb) File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/local.py", line 319, in reraise raise exc File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/local.py", line 224, in execute_task result = _execute_task(task, data) ^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/core.py", line 119, in _execute_task return func(*(_execute_task(a, cache) for a in args)) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/Software/miniconda3/envs/xarray-tests-py311/lib/python3.11/site-packages/dask/array/core.py", line 120, in getter c = a[b] ~^^^ File "/Users/spencer/software/xarray/xarray/core/indexing.py", line 490, in __getitem__ result = self.array[self.indexer_cls(key)] ~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/software/xarray/xarray/core/indexing.py", line 699, in __getitem__ return type(self)(_wrap_numpy_scalars(self.array[key])) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/software/xarray/xarray/core/indexing.py", line 686, in __init__ self.array = _wrap_numpy_scalars(as_indexable(array)) ^^^^^^^^^^^^^^^^^^^ File "/Users/spencer/software/xarray/xarray/core/indexing.py", line 727, in as_indexable raise TypeError(f"Invalid array type: {type(array)}") TypeError: Invalid array type: <class 'cftime._cftime.Datetime360Day'>

Thanks for noting that @dcherian -- I think I got to all of them now.

Thanks all for the help! Fingers crossed things should be all green now. Happy to address any more review comments.

I'm fine waiting until #7724 is merged to let our main CI cover this. Indeed the upstream tests are flaky. Locally I just installed pandas 2 via pip to do testing during development.

Thanks @keewis -- I went ahead and added some more detailed thoughts about how to address the cftime-related failures in #7707. I can try to address those soon.

I think it is fine that CFTimeIndex.to_datetimeindex() no longer raises an error for out-of-nanosecond-precision range dates, so we can simply relax that test for pandas versions greater than or equal to 2 and update its docstring.

In practice I don't think this will come up very often (we can address later if we want probably), but one subtle issue there is that prior to October 15th, 1582, the proleptic Gregorian calendar and the "standard" calendar are not equivalent, so we may want to update how we warn when converting between calendars. The "standard" calendar according to the CF Conventions is a mixed Julian/Gregorian calendar, which uses a Julian calendar prior to 1582-10-15 and a Gregorian calendar after.

In cftime the DatetimeGregorian object conforms to this definition, and is what is created if you provide "standard" as the calendar argument to num2date: ```

cftime.num2date([0, 1], units="days since 1582-10-04", calendar="standard") array([cftime.DatetimeGregorian(1582, 10, 4, 0, 0, 0, 0, has_year_zero=False), cftime.DatetimeGregorian(1582, 10, 15, 0, 0, 0, 0, has_year_zero=False)], dtype=object) cftime.num2date([0, 1], units="days since 1582-10-04", calendar="proleptic_gregorian") array([cftime.DatetimeProlepticGregorian(1582, 10, 4, 0, 0, 0, 0, has_year_zero=True), cftime.DatetimeProlepticGregorian(1582, 10, 5, 0, 0, 0, 0, has_year_zero=True)], dtype=object) ```

I'll need to think more about how to handle the test_should_cftime_be_used_source_outside_range failure; I'm not sure if we're ready to handle changing the behavior of this until we fully address #7493.

The cftime-related ones are still pandas related, but hopefully should be straightforward to address.

Thanks for cleaning up the merge error I must have introduced; I agree this should be ready to go.

Thanks @keewis for fixing this upstream (https://github.com/pandas-dev/pandas/pull/52220)!

Still mulling this over a bit, but one other thing that occurs to me is that if we go with a strftime solution we should be careful (if it exists) to preserve any sub-second information of the Timestamp as well.

Answered in #7610.

Thanks @keewis -- I pushed a fix. The documentation build succeeded. The other build failures look unrelated.

Hmm I guess I should have run the tests before saying anything. We could probably work around the NaT issue fairly easily, but I forgot about timezones. The ability to include a colon in the UTC offset, which we expect in the failing test, with strftime was only recently added to Python (https://github.com/python/cpython/pull/95983), which makes things a little messier. I'll think about this a little more.

Thanks @keewis--great suggestion--I think this should be ready for review now too!

@keewis and @dcherian indeed I just wanted to give it another once over. I think it should be ready for review now. Thanks for your patience on this!

Indeed currently we do not support indexing a CFTimeIndex-backed array using a list of strings, but that's something I think we would be happy to change (e.g. we do accept a list of strings to interp for CFTimeIndex-backed arrays).

For the time being you should be able to use cftime.DatetimeAllLeap values themselves: python ds.sel(time=[cftime.DatetimeAllLeap(2023, 1, 1), cftime.DatetimeAllLeap(2023, 1, 2)])

After debugging locally it looks like the remaining doc build failure was the result of a line I inadvertently deleted further down in the cftime_range docstring. It should be good now I think.

Currently in xarray we make the choice based on the calendar attribute associated with the data on disk (following the CF conventions). If the data has a non-standard calendar (or cannot be represented with nanosecond-precision datetime values) then we use cftime; otherwise we use NumPy. Which kind of calendar do PMIP simulations typically use?

For some background -- my initial need in this realm came mainly from idealized climate model simulations (e.g. configured to start on 0001-01-01 with a no-leap calendar), so I do not have a ton of experience with paleoclimate research. I would be happy to learn more about your application, however!

Thanks for joining the meeting today @khider. Some potentially relevant places in the code that come to my mind are: - Automatic casting to nanosecond precision - Decoding times via pandas - Encoding times via pandas - datetime_to_numeric

Though as @shoyer says, searching for datetime64[ns] or timedelta64[ns] will probably go a long way toward finding most of these issues.

Some design questions that come to my mind are (but you don't need an answer to these immediately to start working): - How do we decide which precision to decode times to? Would it be the finest precision that enables decoding without overflow? - This is admittedly in the weeds, but how do we decide when to use cftime and when not to? It seems obvious that in the long term we should use NumPy values for proleptic Gregorian dates of all precisions, but what about dates from the Gregorian calendar (where we may no longer have the luxury that the proleptic Gregorian and Gregorian calendars are equivalent for all representable times)? - Not a blocker (since this is an existing issue) but are there ways we could make working with mixed precision datetime values friendlier with regard to overflow (https://github.com/numpy/numpy/issues/16352)? I worry about examples like this:

```
>>> np.seterr(over="raise")
>>> np.datetime64("1970-01-01", "ns") - np.datetime64("0001-01-01", "D")
numpy.timedelta64(6795364578871345152,'ns')
```

Great -- I'll plan on joining. That's correct. It is at 8:30 AM PT (https://github.com/pydata/xarray/issues/4001).

I can block out time to join today's meeting or an upcoming one if it would be helpful.

@dcherian +1. I'm happy to engage with others if they are motivated to start on this earlier.

Indeed it would be nice if this "just worked" but it may take some time to sort out (sorry that this example initially got your hopes up!). Here what I mean by "address" is continuing to prevent non-nanosecond-precision datetime values from entering xarray through casting to nanosecond precision and raising an informative error if that is not possible. This of course would be temporary until we work through the kinks of enabling such support. In the big picture it is exciting that pandas is doing this in part due to your grant.

Thanks for posting this general issue @khider. This is something that has been on my radar for several months and I'm on board with it being great to support (eventually this will likely help cftime support as well).

I might hesitate to say that I'm actively working on it yet 😬. Right now, in the time I have available, I'm mostly trying to make sure that xarray's existing functionality does not break under pandas 2.0. Once things are a little more stable in pandas with regard to this new feature my plan is to take a deeper dive into what it will take to adopt in xarray (some aspects might need to be handled delicately). We can plan on using this issue for more discussion.

As @keewis notes, xarray currently will cast any non-nanosecond precision datetime64 or timedelta64 values that are introduced to nanosecond-precision versions. This casting machinery goes through pandas, however, and I haven't looked carefully into how this is behaving/is expected to behave under pandas 2.0. @khider based on your nice example it seems that it is possible for non-nanosecond-precision values to slip through, which is something we may need to think about addressing for the time being.

Thanks for the help @keewis and @dcherian. I think this should be ready for a full review.

(Obviously the above comment only applies if cftime installed; technically cftime is an optional dependency)

@yurivict could you show the output of xarray.show_versions()? I suspect you are using an old version of cftime. Xarray requires a minimum version of 1.5.

@jhamman nice -- thanks for noting that. I think you can replace timedeltas with unique_timedeltas here:

https://github.com/pydata/xarray/blob/b21f62ee37eea3650a58e9ffa3a7c9f4ae83006b/xarray/coding/times.py#L365

and then remove any reference to timedeltas in the rest of the function.

Issue reported here: https://github.com/pandas-dev/pandas/issues/50867. Thanks again for noting that @keewis. I was too narrowly focused on fixing this in xarray.

We may still want to be careful about non-nanosecond-precision Timestamp objects leaking into our code for the time being, but that's a different conversation.

Thanks @keewis.

I think this is a bug in pandas: pd.Timestamp("1-01-01 00:00:00") returns a date in 2001.

I suspect this may be related to the comment here.

In any case, t.isoformat(sep=" ") should return the year with 4 digits so maybe we should use that instead? (once again, it doesn't, but that I think is also a bug, maybe the same one?)

Though indeed the isoformat behavior is definitely a bug and it would be fair to at least expect this kind of roundtrip to hold, but it doesn't: ```

import pandas as pd pd.Timestamp(str(pd.Timestamp("0001-01-01"))) Timestamp('2001-01-01 00:00:00') ```

I can report that to pandas. It's not necessarily surprising considering that it previously was not possible to write pd.Timestamp("0001-01-01").

After a git bisect I traced the units encoding errors to https://github.com/pandas-dev/pandas/pull/49737. I made a PR to address this in xarray (#7441), but it turns out this change may be reverted in pandas (https://github.com/pandas-dev/pandas/issues/50704) so I would hold off on reviewing it until a decision is made there.

7373 should take care of the date_range errors and https://github.com/pydata/xarray/issues/7266#issuecomment-1336172648 is my thinking on how to address the base deprecation. I'll have to look into what's causing the other datetime-related errors.

Apologies I haven't much time to work on these issues lately. Hopefully some time should free up for me towards the end of next week.

Yes, I think that's what we'll end up going with @mroeschke. Thanks for your patience as we worked it out.

Technically I think we would consider get_loc to be public (both due to how it is named and in that it is at least primitively documented). That said, I don't know how many people use the lower level methods of CFTimeIndexes like get_loc outside of xarray data structures.

My inclination is that it would not require a deprecation cycle, given the warning had already been being raised as in the example above. It's true though that this warning was not being emitted in the case that key was a string, so perhaps it's a little borderline. Maybe @dcherian or @mathause, do you have any thoughts?

As I think about this more, it wouldn't be too hard for us to support the base argument even after pandas removes it, so perhaps this isn't so urgent (at least as far as deprecation is concerned; we still need to make updates for compatibility, however). The code to translate a base argument to an offset argument can be found here, and is all possible with public API functionality. I already did something similar for the CFTimeIndex resampling code in #7284.

Maybe you were already thinking along those lines @dcherian.

Sorry I didn't get to adding the warning today. I'll try and put something together over the weekend so that it gets into the release after today's. I'm not sure exactly when pandas 2.0 will be out, but regardless I guess at least it could still be valuable for anyone who doesn't upgrade xarray and pandas at the same time.

Should we add some sort of deprecation warning regarding the use of the base argument with future versions of pandas before the next release?

(I did not end up restoring the pandas warning in #7284)

Perhaps we can at least restore the warning in #7284 in case our next release happens to take place before the next pandas release to give users somewhat of a heads up. Apologies for being a bit out of the loop of #4292 at the time.

I went ahead and actually implemented the origin and offset options for the CFTimeIndex version of resample as part of #7284. It might be good to finish that and then we can decide how we would like to handle the deprecation.

Maybe we should try to catch Deprecation warnings in the nightly builds and raise an error / Automatic issue, so we can fix things before they break.

I agree -- something like that would be useful in general. In this particular case it seems like we were aware of it at one point, but just lost track after silencing it initially for compatibility reasons (https://github.com/pydata/xarray/pull/4292#issuecomment-691665611). Unfortunately that means that this was silenced in user code as well.

Thanks for pinging this issue again in #7266 -- apologies for missing it earlier. I think we should adapt to the new inclusive argument and implement it for the cftime_range function as well.

Yes, I think so too. I can look into it more this weekend. Clearly we need to change the code that relies on pandas immediately. For resampling with a CFTimeIndex I may create a separate issue for implementing these new arguments (we can probably get those tests passing in the meantime, however).

It looks like I've stumbled upon some new upstream test failures related to https://github.com/pandas-dev/pandas/pull/49101.

Nice find @mathause! That seems to work. I added a version of that function to utils. The warning from my example above now looks like: test_warning.py:5: UserWarning: Converting non-nanosecond precision datetime values to nanosecond precision. This behavior can eventually be relaxed in xarray, as it is an artifact from pandas which is now beginning to support non-nanosecond precision values. This warning is caused by passing non-nanosecond np.datetime64 or np.timedelta64 values to the DataArray or Variable constructor; it can be silenced by converting the values to nanosecond precision ahead of time. da = xr.DataArray([1], dims=["time"], coords=[times])

Hmm...I was kind of hoping we could avoid something like adding a _stacklevel_increment argument, but hardcoding the stacklevel obviously doesn't work in the case that fields are converted to variables within xarray, e.g. in the case of passing coordinates to the DataArray constructor:

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

times = [np.datetime64("2000-01-01", "us")] da = xr.DataArray([1], dims=["time"], coords=[times]) ```

$ python test_warning.py /Users/spencer/software/xarray/xarray/core/dataarray.py:148: UserWarning: Converting non-nanosecond precision datetime values to nanosecond precision. This behavior can eventually be relaxed in xarray, as it is an artifact from pandas which is now beginning to support non-nanosecond precision values. This warning is caused by passing non-nanosecond np.datetime64 or np.timedelta64 values to the DataArray or Variable constructor; it can be silenced by converting the values to nanosecond precision ahead of time. var = as_variable(coord, name=dim)

For good measure I also added some assertions that the dtypes after conversion were all of nanosecond precision.

...which turned up a case that was not valid under the existing pandas version, so I added another separate test.

Thanks @dcherian -- I did my best to try and cover our bases. For good measure I also added some assertions that the dtypes after conversion were all of nanosecond precision.

No worries @djhoese. I should be able to submit a PR implementing https://github.com/pydata/xarray/issues/7175#issuecomment-1281254828 later today.

Thanks @dcherian — yes it probably would be good to create an issue to track further discussion about this. Encoding times with float128 values might be of interest to some. I’ll create one tomorrow.

Yes, I think that’s the safest approach on our side to start out.

This is neat to be honest. I think this is related to more progress on enabling different precision datetimes and timedeltas in pandas (issues and PRs related to the "Non-nano" effort). We'll have to read more about how they plan to roll this out and figure out how we will adapt in xarray.

@durack1 I guess the answer is today thanks to @dcherian :).

Thanks for taking a look @dcherian. This should at least handle any exceptions raised by decode_cf_variable (which could be more than just time-related).

I was able to reproduce this issue in a Docker container using the s390x Debian image. After a little experimentation I narrowed it down to the following minimal example: ```

import numpy as np; import pandas as pd np.version '1.23.3' pd.version '1.4.4' pd.Series(np.array([1]).astype("<M8[h]")) Traceback (most recent call last): File "\<stdin>", line 1, in \<module> File "/usr/local/lib/python3.9/dist-packages/pandas/core/series.py", line 451, in __init__ data = sanitize_array(data, index, dtype, copy) File "/usr/local/lib/python3.9/dist-packages/pandas/core/construction.py", line 570, in sanitize_array subarr = _try_cast(data, dtype, copy, raise_cast_failure) File "/usr/local/lib/python3.9/dist-packages/pandas/core/construction.py", line 729, in _try_cast return sanitize_to_nanoseconds(arr, copy=copy) File "/usr/local/lib/python3.9/dist-packages/pandas/core/dtypes/cast.py", line 1717, in sanitize_to_nanoseconds values = conversion.ensure_datetime64ns(values) File "pandas/_libs/tslibs/conversion.pyx", line 257, in pandas._libs.tslibs.conversion.ensure_datetime64ns File "pandas/_libs/tslibs/np_datetime.pyx", line 120, in pandas._libs.tslibs.np_datetime.check_dts_bounds pandas._libs.tslibs.np_datetime.OutOfBoundsDatetime: Out of bounds nanosecond timestamp: 8220291319602-05-05 16:00:00 This confirms it is an upstream issue. Interestingly if we use the native byte order (big-endian on this architecture) for the dtype, this example works: pd.Series(np.array([1]).astype("M8[h]")) 0 1970-01-01 01:00:00 dtype: datetime64[ns] or more explicitly pd.Series(np.array([1]).astype(">M8[h]")) 0 1970-01-01 01:00:00 dtype: datetime64[ns] ``` It appears the inverse of this issue (big-endian dtype leading to a failure on a little-endian system) came up once in pandas: https://github.com/pandas-dev/pandas/issues/29684. @amckinstry I'm not sure what it will take to fix this issue in pandas, but you are welcome to open an issue there. They may also have a difficult time reproducing and testing this, however (https://github.com/pandas-dev/pandas/pull/30976#issuecomment-573989082).

Oops sorry I missed the merge conflict. Thanks @dcherian and @mathause.

Thanks @mathause for reviewing.

Interesting. Thanks for checking that #6988 indeed solves this. I went ahead and merged it, but when I get a chance I’ll keep trying to track down the root cause of this issue.

I will go ahead and merge this. I think it’s a positive change regardless, and it’s nice it also helps with #6906.

Thanks @amckinstry. I guess my last try to produce a pandas minimal example might be: ```

import numpy as np; import pandas as pd pd.Series(np.array([np.int64(1000000).astype("<M8[h]")])) 0 2084-01-29 16:00:00 dtype: datetime64[ns] or potentially more simply: import numpy as np; import pandas as pd pd.Series(np.int64(1000000).astype("<M8[h]")) 0 2084-01-29 16:00:00 dtype: datetime64[ns] `` Somewhere something is going wrong in converting a non-nanosecond-precision datetime value to a nanosecond-precision one (maybe the cast to apd.Timestamp` in my earlier example was short-circuiting this).

I think #6988 should likely work around this issue at least on the xarray side, since it passes datetime64[ns] values into the DataArray constructor immediately. It also seems like the function where the error occurs (ensure_datetime64ns) was recently eliminated in favor of an updated implementation in pandas, so I wonder if this will be an issue there going forward.

I encountered the promote_attrs / merge issue recently in an example similar to @lanougue's above. I'm with @keewis that I would also tentatively support setting promote_attrs=False in merge, but I'm also not that familiar with that part of the codebase and the decisions that went into it.

Thanks for trying that. Maybe it has to do with casting to a pd.Series. Could you maybe also try: ```

import numpy as np; import pandas as pd pd.Series(pd.Timestamp(np.int64(1000000).astype("<M8[h]"))) 0 2084-01-29 16:00:00 dtype: datetime64[ns] ```

It looks like we can bump to 2021.08.0 instead of dealing with this compatibility code.

Thanks for noting this @dcherian. I wasn't sure if this was worth bumping our minimum version for, but since our policy allows for it, I'm happy to make that change.

Ah, this still didn't quite work for older dask versions:

This appears to be the commit that fixed this issue in dask, meaning it should be fixed as of dask version 2021.07.0: ``` git bisect broken 626cc724ffa1d14cd74f48ff4f21ef443d832afa is the first fixed commit commit 626cc724ffa1d14cd74f48ff4f21ef443d832afa Author: Julia Signell jsignell@gmail.com Date: Tue Jun 29 11:50:36 2021 -0400

Allow mixing dask and numpy arrays in @guvectorize (#6863)

dask/array/core.py | 7 ++++++- dask/array/tests/test_gufunc.py | 11 +++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) ```

Apologies for taking a while to look into this. I have not been able to set up an environment to reproduce these test failures, which makes it tricky. It seems like the tests are failing in the setup step, where a DataArray of some random times is generated: data = xr.DataArray(np.random.randint(1, 1000000, size=(4, 5)).astype("<M8[h]"), dims=("x", "y")) In principle the NumPy code should not generate any times larger than 1,000,000 hours since 1970-01-01, i.e. 2084-01-29T16:00:00, which in theory should be should be representable with a nanosecond-precision pandas Timestamp.

Trying to narrow things down, I guess my first question would be: does the following fail in this environment? Is this maybe a pandas issue? ```

import numpy as np; import pandas as pd pd.Timestamp(np.int64(1000000).astype("<M8[h]")) Timestamp('2084-01-29 16:00:00') ``` I think these tests could be simplified some to remove the randomness, but that's probably a separate issue.

6717 should fix the time decoding errors.

I will investigate these cftime failures more in a bit, but I suspect they may be related to these recent PRs in pandas. It appears they are working towards adding support for non-nanosecond-precision datetimes, which is quite an interesting development.

Many thanks for restoring this functionality @headtr1ck!

Thanks @sappjw -- this is a distillation of the bug derived from your example:

```

import numpy as np import xarray as xr xr.coding.times.decode_cf_datetime(np.uint32(50), "seconds since 2018-08-22T03:23:03Z") array('2018-08-22T03:23:05.755359744', dtype='datetime64[ns]') ```

I believe the solution is to also cast all unsigned integer values -- anything with dtype.kind == "u" -- to np.int64 values here: https://github.com/pydata/xarray/blob/770e878663b03bd83d2c28af0643770bdd43c3da/xarray/coding/times.py#L220-L224 Ordinarily we might worry about overflow in this context -- i.e. some np.uint64 values cannot be represented by np.int64 values -- but I believe since we already verify that the minimum and maximum value of the input array can be represented by nanosecond-precision timedelta values, we can safely do this.

This issue was addressed as part of https://github.com/pydata/xarray/pull/4684.

After a little more investigating I think I've narrowed this down to the following example. This particular example happens to be problematic even with the development version of xarray (i.e. it is an existing bug): ``` In [1]: import pandas as pd ...: import xarray as xr

In [2]: times = pd.date_range("2000", periods=2)

In [3]: variable = xr.Variable(["time"], times)

In [4]: variable.encoding Out[4]: {}

In [5]: xr.conventions.decode_cf_variable(variable).encoding Out[5]: {'dtype': dtype('<M8[ns]')} `` The issue is that callingdecode_cf_variableon a variable that already has adatetime64[ns]dtype adds an encoding attribute that leads to problems downstream ('<M8[ns]'` is not a valid dtype in a netCDF file).

In this case it seems like we'd prefer not to add or modify any encoding attributes at all -- we are perfectly capable of handling empty encoding dictionaries for datetime-like data when writing it out to disk.

Thanks for the report! This is indeed odd behavior. I was able to reproduce this issue with the latest release of xarray, 2022.3.0, but interestingly when I try running the example with the development version of xarray, everything works:

``` DATASET1: {'long_name': 'time since initialization', 'field': 'time, scalar, series'} {'source': 'https://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/CBOFS/MODELS/2022/04/09/nos.cbofs.regulargrid.n001.20220409.t00z.nc', 'original_shape': (1,), 'dtype': dtype('float64'), 'units': 'seconds since 2016-01-01 00:00:00', 'calendar': 'gregorian'} datetime64[ns] {'source': 'https://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/CBOFS/MODELS/2022/04/09/nos.cbofs.regulargrid.n001.20220409.t00z.nc', 'original_shape': (1,), 'dtype': dtype('float64'), 'units': 'seconds since 2016-01-01 00:00:00', 'calendar': 'gregorian'}

DATASET2: {'long_name': 'time since initialization', 'units': 'seconds since 2016-01-01 00:00:00', 'calendar': 'gregorian', 'field': 'time, scalar, series'} {'source': 'https://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/CBOFS/MODELS/2022/04/09/nos.cbofs.regulargrid.n001.20220409.t00z.nc', 'original_shape': (1,), 'dtype': dtype('float64')} float64 {'source': 'https://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/CBOFS/MODELS/2022/04/09/nos.cbofs.regulargrid.n001.20220409.t00z.nc', 'original_shape': (1,), 'dtype': dtype('float64'), 'units': 'seconds since 2016-01-01 00:00:00', 'calendar': 'gregorian'}

DATASET3: {'long_name': 'time since initialization', 'field': 'time, scalar, series'} {'source': 'https://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/CBOFS/MODELS/2022/04/09/nos.cbofs.regulargrid.n001.20220409.t00z.nc', 'original_shape': (1,), 'dtype': dtype('float64'), 'units': 'seconds since 2016-01-01 00:00:00', 'calendar': 'gregorian'} datetime64[ns] {'source': 'https://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/CBOFS/MODELS/2022/04/09/nos.cbofs.regulargrid.n001.20220409.t00z.nc', 'original_shape': (1,), 'dtype': dtype('float64'), 'units': 'seconds since 2016-01-01 00:00:00', 'calendar': 'gregorian'} ```

It's good that this appears to be fixed, but it perhaps suggests we have some missing test coverage. I might see if I can develop a simple test we can add to our test suite.

I agree this should be better documented.

Here when I want to retrieve the year of datetime, instead of casting back to an array of object & using datetime.year, it's handy to use built-in numpy datetime64 conversion.

The recommended way to extract datetime components, like the year, is to use the DatetimeAccessor. For example if you have a DataArray of datetime-like values (whether they are of type np.datetime64[ns] or cftime.datetime) you can do something like this:

da.dt.year

and it will return a DataArray containing the year of each datetime (for more information see the "Datetime components" section of the documentation).

This is a common point of confusion, but is in fact expected. Xarray intentionally converts any np.datetime64 type to np.datetime64[ns]. The primary motivation is compatibility with pandas, which xarray relies on for time indexing and other time-related operations through things like pandas.DatetimeIndex or the pandas.Series.dt accessor (see, e.g., discussion in https://github.com/pydata/xarray/issues/789).

May I ask what your reason is for requiring a lower-precision datetime type? In xarray we have tried to provide alternatives (e.g. through cftime) for use-cases like longer date ranges.

Thanks again @okhoma -- welcome to xarray!

I did a little more digging. I'm not a backend expert, but I think the issue can be distilled to the following. When passed an ordinary NumPy array as an input, xarray.Variable will automatically cast any np.datetime64 values to nanosecond precision: ``` In [3]: arr = np.array([np.datetime64("2000-01-01", "us")])

In [4]: xarray.Variable(["time"], arr) Out[4]: <xarray.Variable (time: 1)> array(['2000-01-01T00:00:00.000000000'], dtype='datetime64[ns]') However, if passed a `LazilyIndexedArray`, this casting will not occur (note the dtype is still `'datetime64[us]'`): In [5]: lazily_indexed_arr = xarray.core.indexing.LazilyIndexedArray(arr)

In [6]: xarray.Variable(["time"], lazily_indexed_arr) Out[6]: <xarray.Variable (time: 1)> array(['2000-01-01T00:00:00.000000'], dtype='datetime64[us]') `` This is based on [this code](https://github.com/pydata/xarray/blob/f42ac28629b7b2047f859f291e1d755c36f2e834/xarray/backends/zarr.py#L404-L419) withinxarray.backends.zarr.ZarrStore`.

The casting does not occur in the Variable constructor, because it requires the type of the array be np.ndarray -- see here. Regardless, even if we relaxed that, _possibly_convert_objects would raise an error, because it is not compatible with LazilyIndexedArray objects.

Is this something that we maybe need to address within the xarray zarr backend? As I understand it, zarr is a bit unusual compared to other storage formats we deal with in that it can store and open np.datetime64 data directly -- normally datetime data starts as some numeric type and is converted to datetime data after going through our decoding pipeline.

To be honest I didn't know it was possible to open a Dataset and maintain np.datetime64[us] values. I feel like casting maybe should occur automatically -- we do this already in other contexts, e.g.: ``` In [3]: da = xr.DataArray([np.datetime64("2000-01-01", "us")], dims=["time"], name="time")

In [4]: da Out[4]: <xarray.DataArray 'time' (time: 1)> array(['2000-01-01T00:00:00.000000000'], dtype='datetime64[ns]') Dimensions without coordinates: time ``` I'll try and dig deeper into this in the next few days, but @keewis's workaround should be good in the meantime.

Awesome, thanks @antarcticrainforest -- yes, it will always have an object dtype.

Thanks @antarcticrainforest -- I think that's exactly what @mathause is getting at. It seems fairly safe to add a new keyword argument to CFTimeIndex.__new__, and the example @mathause uses would make a nice test. Would either of you be up to make a PR?

Great! I'll go ahead and close this issue.

Thanks @pierreloicq!

@dcherian @mathause @andersy005 given a DataArray of times with missing values, do you have any thoughts on what the preferred result of da.dt.season would be?

One option would be to return np.nan in place of the missing time values:

``` In [3]: times = [np.datetime64("NaT"), np.datetime64("2000-01-01")]

In [4]: da = xr.DataArray(times, dims=["x"])

In [5]: da.dt.season Out[5]: <xarray.DataArray 'season' (x: 2)> array([nan, 'DJF'], dtype=object) Dimensions without coordinates: x ```

This would be somewhat in line with how pandas handles this in other contexts (e.g. https://github.com/pydata/xarray/pull/5876#discussion_r734447120). But this sort of awkwardly returns a DataArray of mixed types. Another option, and this is how @pierreloicq implemented things originally, would be to simply return a string label for missing values, e.g.

In [5]: da.dt.season Out[5]: <xarray.DataArray 'season' (x: 2)> array(['nan', 'DJF'], dtype='<U32') Dimensions without coordinates: x

As I've thought about this more, this feels nicer, because the types are consistent and "season" is just a category label, and "nan" categorizes these values just as well as np.nan. Do you agree?

Something that might make this easier would be to implement a now method on cftime.datetime objects. I would imagine that could have use in more than just xarray.

This seems reasonable. We would probably only enable it for real-world calendars, i.e. "standard", "gregorian", "proleptic_gregorian", or "julian". "julian" is a little trickier since it is not merely casting from datetime.datetime.now() to a cftime date (see Conversion between Julian and Gregorian calendars).

For shift intervals that can be represented as timedeltas this seems reasonably straightforward to add. I would hold off for monthly or annual intervals -- even for 360-day calendars, I don't think that non-integer shift factors are very well-defined in that context, since those frequencies involve rounding, e.g. to the beginnings or ends of months: ``` In [2]: times = xr.cftime_range("2000", freq="7D", periods=7)

In [3]: times Out[3]: CFTimeIndex([2000-01-01 00:00:00, 2000-01-08 00:00:00, 2000-01-15 00:00:00, 2000-01-22 00:00:00, 2000-01-29 00:00:00, 2000-02-05 00:00:00, 2000-02-12 00:00:00], dtype='object', length=7, calendar='gregorian', freq='7D')

In [4]: times.shift(2, "M") Out[4]: CFTimeIndex([2000-02-29 00:00:00, 2000-02-29 00:00:00, 2000-02-29 00:00:00, 2000-02-29 00:00:00, 2000-02-29 00:00:00, 2000-03-31 00:00:00, 2000-03-31 00:00:00], dtype='object', length=7, calendar='gregorian', freq='None') ```

Sorry for not responding to this issue earlier -- I think this is related to https://github.com/pydata/xarray/issues/5686 (see discussion and links there for more details). I can reproduce your issue with cftime version 1.5.0, and I tested things with cftime version 1.5.1 and it was fixed (i.e. cftime.DatetimeJulian objects are returned in all cases).

Sounds good -- thanks @mathause!

I recall seeing discussion regarding something like this in https://github.com/pydata/xarray/issues/5028 and https://github.com/pydata/xarray/discussions/5275; it seems like upgrading to zarr version 2.8.3 might help.

No worries about the cftime tests by the way -- how to best handle missing values in that context is still somewhat of an open question (https://github.com/Unidata/cftime/issues/145).

@mathause I think the issue is that non-periodic interpolation introduces NaNs into the dataset, which get plotted as a white line.

Xarray's built-in interpolation method does not know how to handle periodic coordinates, which is why you see the blank center line. I would recommend using xESMF, which is an xarray wrapper of some code specifically written to interpolate from one geospatial grid to another, and knows how to handle periodic coordinates (see the periodic flag in the constructor of the Regridder object). It also provides multiple different interpolating methods which, depending on your application, might be more applicable than xarray's options.

This logic has been around in xarray for a long time (I think it dates back to https://github.com/pydata/xarray/pull/12!), so it predates me. If I had to guess though, it would have to do with the fact that back then, a form of cftime.date2num was used to encode all times, even those that started as np.datetime64 values. I think that's significant for two reasons: 1. In the old days, date2num would only return floating point values, even if the times could in principle be encoded with integers. For that reason, for accuracy reasons, it was best to keep the encoded values as small as possible to avoid roundoff error. 2. Even if (1) was not the case back then, date2num did not -- and still does not -- support nanosecond units, because it relies on microsecond-precision datetimes.

This of course is not true anymore. We no longer use date2num to encode np.datetime64 values, and we no longer encode dates with floating point values by default (#4045); we use integers for optimal round-tripping accuracy, and are capable of encoding dates with nanosecond units.

To be honest, currently it seems the only remaining advantage to choosing a larger time encoding unit and proximate reference date is that it makes the raw encoded values a little more human-readable. However, encoding dates with units of "nanoseconds since 1970-01-01" is objectively optimal for np.datetime64[ns] values, as it guarantees the maximum range of possible encoded times, and maximum round-trip accuracy, so it could be worth revisiting our approach in light of the fact that it makes appending somewhat dangerous.

Given the increasing scope of cftime-related utilities, I wonder if it would make sense to try to eventually move the cftime-related utilities outside of Xarray into a separate package? After the current ongoing "explicit indexes" refactor, I believe we will have public APIs for all of the essential functionality.

Yeah, I think that could be something good to strive for, and it's exciting that the explicit indexes refactor could facilitate that. To be clear, are you ok with merging #5233 and waiting until after the refactor is complete to discuss future plans there?

It would be cool if we could merge #5233 before the next release as well. I think it is ready, and @aulemahal has been diligently working on it for the last several months (and been very patient with my sporadic reviewing frequency).

@mathause I think this should be ready now. I added some version-dependent logic to the tests.