Climatology calculation for daily or higher temporal frequency input

[lee1043]

Following #303 (comment), calculation of climatology for daily or higher temporal frequency time series will need to remove data on dates of 29th February in leap years from the input time series.

Calendars having normal leap year in CF-convention are "gregorian", "standard", and "proleptic_gregorian". (Calendars in CF-convention: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch04s04.html)

We may want to implement a kind of pre-processing function for input data, which could be something like below. I am yet to have a good idea where this function should be placed. @tomvothecoder would you share your insight?

# Get daily climatology (ds: dataset with daily time series)
ds_daily_climo = ds.temporal.climatology("tas", freq="day", weighted=True)

# Below is somewhere in the temporal.climatology, in advance to calculate climatology, to pre-process input
def function(ds, freq):
    if freq in ["day", "hour"]:
        if calendar in ["gregorian", "standard", "proleptic_gregorian"]:
            # Remove leap year prior to daily climatology
            ds2 = ds.sel(time=~((ds.time.dt.month==2)&(ds.time.dt.day==29)))
            ds = copy.deepcopy(ds2)
    return ds

[jypeter]

This could be referenced in a calendar page

[tomvothecoder]

Thank you @lee1043 and nice starting point with your pseudo code.

I'll take a closer look soon, maybe this week!

[tomvothecoder]

We can also convert this to an issue and link #303 to it.

[lee1043]

@tomvothecoder I don't mind converting this to issue. Thanks for taking a look.

[msahn]

@lee1043 I have a question about this method. It seems like you calculate 365-day climatology excluding all 2/29. If so, how can you calculate the anomaly of 2/29? In my method, I calculate 366-day climatology including 2/29.

[lee1043]

@msahn In this method 2/29 is not being considered. Presume an example case of getting daily climatology for 10 years, then 2/29 is average of 2 years while other days are average of 10 years. I won't prefer having such inconsistency, so I would prefer dropping 2/29. It looks like different groups have done this differently, and I don't think there is the one right way though.

One possible path could be allowing some flexibility for how to calculate daily climatology by enabling the following options, such as:

1. Drop leap year days (maybe as a default, as like in the PR: Drop leap days based on CF calendar type to calculate daily climatologies and departures  #350. I think Xarray takes a similar approach when converting calendar)
2. Insert leap year day for every year filled with NaN (or equivalent) value -- however this may cause error from cftime when writing calculated climatology to netCDF file, as documented here and here. The calculated daily climatology will be written as year 1 (no leap year), not year 0 (yes leap year). cftime will return error because year 1 has leap year day.
3. Use julian day for climatology calculation using groupby(day of year) (I think NCL takes this approach: related link1, link2) -- however the same issue occurs for the last day of the year (365th in normal year and 366th in leap year)

But I also think the code could be more complicated (thus less maintainable) to enable all options -- this question may go to @tomvothecoder

I know you were taking the approach 2 in his precipitation variability code in PMP, would you like to give a quick check if you take approach 1 how much impact would be made to your result? I think it can be checked with just one or two observations, no need to use the entire suite of CMIP models.

[pochedls]

One thought: could you average 2/28, 2/29 (when it exists), and 3/1, then assign that average value to 2/29 for the purposes of creating the climatology (or something like this). If 2/29 is populated in the climatology, then I think it should be possible to calculate departures.

I think this option and/or dropping 2/29 could be good options; I can see scenarios where you would want 2/29 in the time series.

[msahn]

Thank you for your detailed answer. I see your point. I think the precipitation variability results will be almost the same for methods 1 and 2 because only 5 data from 7305 data (20 years) are different.

[lee1043]

@pochedls I think that is a good idea if dealing variables that have temporal continuity (like temperature), but could mislead analysis for variables like precipitation that has intermittency. I think the averaging could be one option, but we may not want to consider that as a default option.

@msahn I guess so, would be interested to see comparison if you were planning so.

[tomvothecoder]

Thank you everyone for coming up with possible solutions for handling leap days. From what I understand, how this is handled is opinionated and can be influenced by the variable types and analysis scenarios.

For now, I believe we should move forward with method 1 since it provides the most value with the least amount of work ((it's already done in #350). This allows us xcdat to calculate daily climatologies/departures, which means PMP can integrate and thoroughly test the temporal averaging APIs.

But I also think the code could be more complicated (thus less maintainable) to enable all options -- this question may go to @tomvothecoder

I generally agree with this assessment. Usually if an API has many configurable options, the code is less maintainable and usage/behavior is less predictable because of how many possible side-effects can occur.

I propose we revisit how to handle leap days if there is high demand/value in also keeping leap days.

How does this sound?

[lee1043]

@tomvothecoder I agree on that, thanks for sharing your thoughts.