xarray potential inconstistencies with cftime

[sbiner]

I am trying to use xarray with different types of calendars. I made a few tests and wonder if somebody can help me make sense of the results. In my test, I generate a DataArray da time series with a 365_day calendar using cftime. I then write that DataArray in a netCDF file and read it in another DataArray da2

Code Sample, a copy-pastable example if possible

import xarray as xr
import cftime
import numpy as np

# generate data for 365_days calendar
#
units = 'days since 2000-01-01 00:00'
time_365 = cftime.num2date(np.arange(0, 10 * 365), units, '365_day')
da = xr.DataArray(np.arange(time_365.size), coords = [time_365], dims = 'time', name = 'data')

# write dataArray in netcdf and read it in new DataArray
da.to_netcdf('data_365.nc', 'w')
da2 = xr.open_dataset('data_365.nc').data

# try resample da
try:
    mean = da.resample(time='Y').mean()
    print(mean.values)
except TypeError:
    print('got TypeError for da')

# try resample da2
mean = da2.resample(time = 'Y').mean()
print (mean.values)

Problem description

As seen in the code the resampledoes not work for da while it does for da2. The problem is related to the the type of da.time which is cftime.DatetimeNoLeap while da2.time is a datetime64. I thought that xarray is using cftime to make the conversion from time numerical values to dates but it looks to me as if it is not the case.

I wonder if this makes sense or if it is something that should eventually be corrected.

INSTALLED VERSIONS
In [6]: print (cftime.version)
1.0.1

------------------ commit: None python: 3.6.5.final.0 python-bits: 64 OS: Darwin OS-release: 17.7.0 machine: x86_64 processor: i386 byteorder: little LC_ALL: None LANG: fr_CA.UTF-8 LOCALE: fr_CA.UTF-8

xarray: 0.10.8
pandas: 0.23.0
numpy: 1.14.3
scipy: 1.1.0
netCDF4: 1.4.1
h5netcdf: None
h5py: 2.7.1
Nio: None
zarr: None
bottleneck: 1.2.1
cyordereddict: None
dask: 0.17.5
distributed: 1.21.8
matplotlib: 2.2.2
cartopy: None
seaborn: 0.8.1
setuptools: 39.1.0
pip: 10.0.1
conda: 4.5.11
pytest: 3.5.1
IPython: 6.4.0
sphinx: 1.7.4

[spencerkclark]

@sbiner these are all reasonable points of confusion. The current behavior in xarray regarding non-standard calendars is complicated, and we are working toward improving the situation. I've tried to provide a recommended solution based on your example as well as some historical/future context. Apologies for the long-winded answer!

Recommendation

For accurate round-tripping of date types, I would recommend that you run your code to open the dataset with the xarray option enable_cftimeindex set to True (by default it is currently set to False; for why this is the case, see the note regarding the default behavior). In the case of your example this would look like:

In [1]: import cftime

In [2]: import numpy as np

In [3]: import xarray as xr

In [4]: units = 'days since 2000-02-25'

In [5]: times = cftime.num2date(np.arange(7), units=units, calendar='365_day')

In [6]: da = xr.DataArray(np.arange(7), coords=[times], dims=['time'], name='a')

In [7]: da.to_netcdf('data-noleap.nc')

In [8]: with xr.set_options(enable_cftimeindex=True):
    ...:     cftimeindex_enabled = xr.open_dataset('data-noleap.nc')
    ...:

Here we can see that the time index is a CFTimeIndex, and that the time coordinate contains instances of cftime.DatetimeNoLeap objects (as they were in the original DataArray we created):

In [9]: cftimeindex_enabled.indexes['time']
Out[9]:
CFTimeIndex([2000-02-25 00:00:00, 2000-02-26 00:00:00, 2000-02-27 00:00:00,
             2000-02-28 00:00:00, 2000-03-01 00:00:00, 2000-03-02 00:00:00,
             2000-03-03 00:00:00],
            dtype='object', name=u'time')

In [10]: cftimeindex_enabled.time[0]
Out[10]:
<xarray.DataArray 'time' ()>
array(cftime._cftime.DatetimeNoLeap(2000, 2, 25, 0, 0, 0, 0, 6, 56), dtype=object)
Coordinates:
    time     object 2000-02-25 00:00:00

Note that resample along a CFTimeIndex has not been implemented yet (#2191). Attempting to do so will raise an error. If you are interested in computing something as simple as a time series of annual means using resample, then you could work around that in the meantime by using groupby, for example:

In [11]: cftimeindex_enabled.groupby('time.year').mean('time')
Out[11]:
<xarray.Dataset>
Dimensions:  (year: 1)
Coordinates:
  * year     (year) int64 2000
Data variables:
    a        (year) float64 3.0

For more information on what is enabled and what is not enabled when using a CFTimeIndex for indexing, see this section in the documentation.

Default behavior

The default behavior can be traced back to the early days of xarray (see the original discussion in #118, #121, and #126). It boils down to coercing any dates decoded into cftime.datetime objects (formerly netCDF4.datetime) into np.datetime64[ns] objects whenever possible. If this coercion is not possible (e.g. a date in the file is not allowed in the standard calendar like 2000-02-30 in the case of a 360-day calendar, or a date has a year outside the range 1678-2262) then cftime.datetime objects are allowed to remain. In other words, by default xarray indeed does use cftime to decode the dates; however, after decoding, it will try its hardest to convert those dates into a friendly type for pandas.

The advantage of the default approach is that, when possible, it allows you to take advantage of all the nice features that a time coordinate indexed by a pandas.DatetimeIndex provides (like resample). A disadvantage is that the dates in memory may not have the same calendar type as those encoded in the file (e.g. if the dates in the file are from a non-standard calendar, like no leap). For operations that rely on computing differences between dates (e.g.differentiate orinterp involving a time coordinate), this can lead to subtle (and silent) errors. Therefore when using dates coerced into a DatetimeIndex from a non-standard calendar, one should use caution to only do operations that are independent of the calendar type (one notable exception here is that xarray does make an effort to encode these dates accurately when writing out to a netCDF file).

Connecting back to your example, we can see that if we don't open the dataset with enable_cftimeindex=True, the dates are coerced to np.datetime64 objects and a DatetimeIndex is used:

In [12]: default = xr.open_dataset('data-noleap.nc')

In [13]: default.indexes['time']
Out[13]:
DatetimeIndex(['2000-02-25', '2000-02-26', '2000-02-27', '2000-02-28',
               '2000-03-01', '2000-03-02', '2000-03-03'],
              dtype='datetime64[ns]', name=u'time', freq=None)

In [14]: default.time[0]
Out[14]:
<xarray.DataArray 'time' ()>
array(951436800000000000L, dtype='datetime64[ns]')
Coordinates:
    time     datetime64[ns] 2000-02-25

In this case, as noted above, resample works:

In [15]: default.resample(time='Y').mean('time')
Out[15]:
<xarray.Dataset>
Dimensions:  (time: 1)
Coordinates:
  * time     (time) datetime64[ns] 2000-12-31
Data variables:
    a        (time) float64 3.0

Future behavior

In xarray we are slowly working towards better support for operations involving cftime.datetime objects (see #789, #1084, #1252, #2008, #2142, #2301, #2434). Eventually we would like to switch to using enable_cftimeindex=True as the default: in that case the behavior would be to use np.datetime64 objects (associated with DatetimeIndexes) only for standard calendars, and cftime.datetime objects (associated with CFTimeIndexes) for any other calendar types.

The two major outstanding issues on this front are probably:

• Adding resample functionality to CFTimeIndex (Adding resample functionality to CFTimeIndex #2191)
• Plotting data with cftime.datetime coordinate axes in matplotlib or holoviews (holoviews / bokeh doesn't like cftime coords #2164).

Once those two remaining issues are addressed, one should be able to do most of the significant things one can do with np.datetime64 dates with cftime.datetime dates (and therefore changing the default behavior would be justified).

[shoyer]

@spencerkclark do you think it would make sense to enable cftimeindex by default in the next major xarray release? On the whole I think it's probably a win for usability at this point...

[sbiner]

@spencerkclark I made tests with enable_cftimeindex=True and noticed I got an error. I did not fit my need but on second thought it's probably better to have this than unexpected resample behaviour, escpecially if we use it to upsample (i.e. from lower frequency toward higher).

Thanks for the complete answer.

[spencerkclark]

@shoyer I agree that seems like a good idea at this stage. Now that there are a number of functions in xarray that do depend differences in dates (as @sbiner notes upsampling with resample, interp, and differentiate), which did not exist in the past, it is perhaps better that things error by default, rather than silently return potentially incorrect results if they have not yet been implemented for dates from non-standard calendars. Users can explicitly opt in to the old workaround if they feel it would be safe in their use cases.

[spencerkclark]

With #2516 already in released versions of xarray, and #2593 and #2665 recently merged, this situation has been significantly improved. I think it is safe now to close this general issue. @sbiner thanks for starting this conversation; feel free to post other issues related to cftime if they come up.

[hafez-ahmad]

How can I canvert julian to dmy index or datetime in pandas? , I have following dataset
Dimensions:
id: 170time: 1560
Coordinates:
time
(time)
int64
67935 67966 67994 ... 115355 115385
array([ 67935, 67966, 67994, ..., 115324, 115355, 115385], dtype=int64)

[spencerkclark]

@hafez-ahmad could you provide more detail about your dataset? Does the "time" coordinate have associated "calendar" and "units" attributes?

[hafez-ahmad]

Here screenshot of my data

here dataset link https://portal.grdc.bafg.de/grdcdownload/external/d94b2ca6-90c3-4220-b0ae-03250f809afe/2021-02-26_08-30.zip
I like to convert julian date to normal calendar datetime like (22-01-2020)

thank you
Hafez

[spencerkclark]

Could you show me what the output of ds.info() looks like for the dataset?

[keewis]

@spencerkclark: if the dataset linked to above is a good example, that would be

ds.info()

xarray.Dataset {
dimensions:
	id = 37 ;
	time = 1380 ;

variables:
	float32 runoff_mean(time, id) ;
		runoff_mean:units = m3/s ;
		runoff_mean:long_name = GRDC calculated from daily data ;
	float32 flag(time, id) ;
		flag:long_name = percentage of valid values used for calculation from daily data ;
	float32 area(id) ;
		area:units = km2 ;
		area:long_name = catchment area ;
	object country(id) ;
		country:long_name = country name ;
		country:iso2 = ISO 3166-1 alpha-2 - two-letter country code ;
	float32 geo_x(id) ;
		geo_x:units = degree_east ;
		geo_x:long_name = station longitude (WGS84) ;
	float32 geo_y(id) ;
		geo_y:units = degree_north ;
		geo_y:long_name = station latitude (WGS84) ;
	float32 geo_z(id) ;
		geo_z:units = m ;
		geo_z:long_name = station altitude (m above sea level) ;
	object owneroforiginaldata(id) ;
		owneroforiginaldata:long_name = Owner of original data ;
	object river_name(id) ;
		river_name:long_name = river name ;
	object station_name(id) ;
		station_name:long_name = station name ;
	float32 timezone(id) ;
		timezone:units = 00:00 ;
		timezone:long_name = utc offset, in relation to the national capital ;
	int64 time(time) ;
		time:long_name = time ;
	int32 id(id) ;
		id:long_name = grdc number ;

// global attributes:
	:title = Mean daily discharge (Q) ;
	:Conventions = CF-1.7 ;
	:references = grdc.bafg.de ;
	:institution = GRDC ;
	:history = Download from GRDC Database, 26/02/2021 ;
	:missing_value = -999.000 ;
}

so no units or calendar attributes on time

[spencerkclark]

Thanks @keewis.

@hafez-ahmad by Julian date do you mean that the time coordinate represents "days since -4713-01-01T12:00:00" in a Julian calendar?

Once we know the units (expressed as "{time_unit} since {reference_date}") and the calendar of the time coordinate, we can convert it to datetime objects via something like the following:

units = "days since -4713-01-01T12:00:00"
calendar = "julian"
ds["time"] = ds.time.assign_attrs(units=units, calendar=calendar)
ds = xr.decode_cf(ds)

I'll admit though, with the values in your dataset, this assumption produces dates like cftime.DatetimeJulian(-4527, 1, 30, 12, 0, 0, 0), which feel unlikely to be correct. Perhaps you are using a different reference date?

[hafez-ahmad]

@spencerkclark
I like to convert datetime. My dataset time is not familiar. It looks like 456852,85993, is there anyvway to convert 456852 to dmy [01-01-2020]?

[spencerkclark]

@hafez-ahmad yes, I'm trying to help, but in order to do that I need more information. What does 456852 represent?

[hafez-ahmad]

@spencerkclark
67935 67966 67994 ... 115355 115385 are actually time. But they are stored aa integer. I like to convert all integer to date. My attached data has exact same time.

[keewis]

the issue is that without more information (units, reference date, calendar), xarray can't decode these integers correctly. For example, if the units are days your data would span more than 100 years, but if we assume milliseconds that would be about 2 hours. Same for reference date and calendar: if they are chosen incorrectly, the decoded data will be completely wrong.

You will have to check either the data provider's website or contact them and ask for help (it seems they started offering their data as netcdf files about two weeks ago so this might actually be a issue with their conversion code). Once you have that information we can help you with the code necessary to convert using xarray.

[jsta]

I believe the dates assocsiated with this particular dataset are days since "1700-01-01"

86287 == 1936-04-01

ds["time"] = ds.time.assign_attrs(units="days since 1700-01-01")
ds = xr.decode_cf(ds)