DISCUSSION: Datetime Alignment Standards

[bradyrx]

CC @kpegion, @aaronspring, @ahuang11, @jmunroe, @Thomas-Moore-Creative, @dougiesquire, @jukent. Looking for insight on time stuff here if/when folks have time to think about this!

With the soon-to-be-merged https://github.com/bradyrx/climpred/pull/294, it's time to figure out some assumptions and standards for dealing with pesky datetime in aligning initialized predictions. I want to work out a system here so we can slowly introduce PRs to handle some of these ideas.

I think we should also add a dedicated docs page on all this datetime stuff. We could draw from some of the thoughts from @jukent in https://ncar.github.io/xdev/posts/time/.

FYI I have some code snippets at https://github.com/bradyrx/esmtools/blob/master/esmtools/temporal.py that could be useful for the below cases.

Here are some thoughts/issues I think are important to discuss (I am going to break them up into separate comments):

1. Do we assume that "monthly" units in the lead attribute aligns with the start of the next month? Currently we handle monthly units by using .shift() with cftime of nlags 'MS' (month-start).

import cftime
import xarray as xr

inits = [cftime.DatetimeProlepticGregorian(y, m, 1) for y in [1995, 1996] for m in [1, 2, 3]]
inits = xr.CFTimeIndex(inits)  # Includes leap year in 1996

print(inits)
>>> CFTimeIndex([1995-01-01 00:00:00, 1995-02-01 00:00:00, 1995-03-01 00:00:00,
                 1996-01-01 00:00:00, 1996-02-01 00:00:00, 1996-03-01 00:00:00],
                 dtype='object')

inits.shift(1, 'MS')  # Current practice
>>> CFTimeIndex([1995-02-01 00:00:00, 1995-03-01 00:00:00, 1995-04-01 00:00:00,
                 1996-02-01 00:00:00, 1996-03-01 00:00:00, 1996-04-01 00:00:00],
                 dtype='object')

inits.shift(1, 'M')  # Month-end?
>>> CFTimeIndex([1995-01-31 00:00:00, 1995-02-28 00:00:00, 1995-03-31 00:00:00,
                 1996-01-31 00:00:00, 1996-02-29 00:00:00, 1996-03-31 00:00:00],
                 dtype='object')

So currently we use the MS practice which means we don't have to worry about leap years or leap-year calendars. We just align with a verification product on the first of the next month for a monthly forecast. Of course there are issues.

• What if the verification product represents the e.g. January average with January 31st, as do many climate models? How do we know when to then use .shift(1, 'M') instead of .shift(1, 'MS')? Do we just ask that users convert their verification product to MS?
• What if they represent with month-middle? Do you have a lookup table that relates ndays in month including leap years so you can shift in units of days?
• Do we accommodate time bounds? The two above decisions are dependent on time bounds. Is the January 31st index on a verification product the January 15th to February 15th average or the January 1st to January 31st average? (This is a common issue for any time resolution)

[bradyrx]

2. Similarly to above, how do we handle yearly forecasts in a non noleap calendar? Currently the practice is to use YS similar to MS above.

import cftime
import xarray as xr

inits = xr.cftime_range('1994', '1996', freq='YS')
inits = xr.CFTimeIndex(inits)
print(inits)
>>> CFTimeIndex([1994-01-01 00:00:00, 1995-01-01 00:00:00, 1996-01-01 00:00:00], dtype='object')

# Align with verification for next year.
inits.shift(1, 'YS')  # Current way we handle it 
>>> CFTimeIndex([1995-01-01 00:00:00, 1996-01-01 00:00:00, 1997-01-01 00:00:00], dtype='object')

# Initialization on a specific day with year forecasts.
inits = [cftime.DatetimeProlepticGregorian(1995, 3, 1),
         cftime.DatetimeProlepticGregorian(1996, 3, 1)]
inits = xr.CFTimeIndex(inits)
print(inits)
>>> CFTimeIndex([1995-03-01 00:00:00, 1996-03-01 00:00:00], dtype='object')
inits.shift(365, 'D')  # Perhaps their verification is aligned one year out.
>>> CFTimeIndex([1996-02-29 00:00:00, 1997-03-01 00:00:00], dtype='object')

Similarly we use YS here to just assume that the 1995 average is written as '1995-01-01' in the verification product. Same issues as above here.

• How do we deal with time bounds?
• How do we deal with leap years? -- See the second example. Maybe they initialize Nov. 1st and want to align with the Nov 1st entry from the verification product. Using 365 gets out of sync here due to leap years (although I think we could just have a variable 365 or 366 addition depending on if it's a leap year and they're using a leap year calendar).

[bradyrx]

3. How do we deal with integer decimal time?

Currently, if an init coordinate comes in with integer or floats, we just assume that it's annual resolution initializations and convert to year-start for that integer.

E.g.

import numpy as np
import pandas as pd
from climpred.utils import convert_time_index

init = np.arange(1990, 1995)
da = xr.DataArray(np.random.rand(len(init)), dims='time', coords=[init])

print(convert_time_index(da, 'time', '').time.to_index())
>>> CFTimeIndex([1990-01-01 00:00:00, 1991-01-01 00:00:00, 1992-01-01 00:00:00,
                 1993-01-01 00:00:00, 1994-01-01 00:00:00],
                 dtype='object', name='time')

So here we just take floats/ints and essentially append on -01-01. But what if the values coming in are decimal time (e.g. 1990.3) or netcdf style integers that correspond to "days since ...". Should we have a try statement that tries to decode_cf on integer/float times? What if it comes in as 19950101 and could just be converted that way?

[bradyrx]

Alright that's all I've got for now! Just wanted to open the discussion with people who have thought about these things so we can find a standard and start chipping away at them with PRs.

[aaronspring]

Regarding 3.: is there any data center that provides fractional float inits? We could simply not allow this... easiest way.

[kpegion]

I'm thinking these issues are big enough that it would be helpful to collect information from the various hindcast datasets currently out there and then see what is common and what is different between them and what is clearly defined and what is ambiguous wrt to handling of time and time-related metadata. I think this will help us see how to answer these questions robustly. At the very least, it will help me get a big picture view of the issues across the various communities and envision a robust set of solutions. Where is the best place to collect all this information?

Once we agree on a location, I can start populating information for NMME, SubX, International S2S database, and the NCEP/UFS hindcasts (with the idea that we keep weather timescales in mind for the future as well). I think we have the potential to help define a more common and consistent set of metadata and definitions across the entire Earth System Prediction community with respect to time. We are planning for the future for SubX and NMME databases so there is also an opportunity here to influence how they will be setup in the future.

[aaronspring]

I think having an overview of the different product names, eg. NMME, SubX, CESM-DPLE, MPI-ESM MiKlip, ... and the details and what this means in the context of climpred might be nice and more transparent for all users.
however, this would take same of effort for all of us to document.

[aaronspring]

#332 requires Leap or NoLeap. it cannot handle the other calendars yet. for now it checks dayofyear in resampling uninitialized.

[bradyrx]

This goal is worked at in https://github.com/bradyrx/climpred/pull/332 but not resolved.

[aaronspring]

a comment when implementing .plot():

need to align with discussion about init times: I propose:

• initialization: as the exact date of the start of the forecast
• observation/hist time: as the middle of time spanned
• lead: as ints but when converted to real time as middle of time