[Enhancement]: use CF information to compute seasonal means

[oliviermarti]

Is your feature request related to a problem?

Short story

Custom seasonnal means doesn't work with my model outputs :-(

Long story

IPSL Earth System Model raw outputs have one time axis time_counter, but two time variables : time_counter and time_centered. (outputs on ESGF works fine)

So if I run the following code :

import xcdat as xc

input_data = "https://thredds-su.ipsl.fr/thredds/dodsC/tgcc_thredds/work/p86mart/IPSLCM6/PROD/Holocene/TR6kCM6AS-Sr01/ATM/Analyse/TS_MO/TR6kCM6AS-Sr01_20000101_29991231_1M_psol.nc"

ds = xc.open_dataset (input_data, use_cftime=True, decode_times=True, decode_cf=True).isel({'time_counter':slice(0,120)})
ds = ds.bounds.add_missing_bounds()

# Custom seasons in a three month format:
custom_seasons = [ ['Dec', 'Jan', 'Feb'], ]

season_config = {'custom_seasons': custom_seasons, 'dec_mode': 'DJF', 'drop_incomplete_seasons':True}

ds.temporal.group_average ("psol", "season", season_config=season_config)

I get :

 File "/Users/marti/Unix/XCDAT/xcdat/temporal.py", line 1163, in
 _subset_coords_for_custom_seasons
     coords_by_month = ds.time.groupby(f"{self.dim}.month").groups
                      ^^^^^^^

A time variable is not identified.

Describe the solution you'd like

Proposed correction

I don't know what could be the side effects, but I've succesfully tested the following change in xcdat/temporal.py :

diff --git a/xcdat/temporal.py b/xcdat/temporal.py
index 2681d1d..2353e16 100644
--- a/xcdat/temporal.py
+++ b/xcdat/temporal.py
@@ -1160,7 +1160,8 @@ class TemporalAccessor:
         """
         month_ints = sorted([MONTH_STR_TO_INT[month] for month in months])

-        coords_by_month = ds.time.groupby(f"{self.dim}.month").groups
+        coords_by_month = ds.cf[['T']].groupby(f"{self.dim}.month").groups
+
         month_to_time_idx = {
             k: coords_by_month[k] for k in month_ints if k in coords_by_month
         }

cf_xarray correctly finds the CF time variable (probably from the attribute axis: T), and the xcdat code runs smoothly.

Describe alternatives you've considered

No response

Additional context

Versions

Version python : 3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:45:13) [Clang 16.0.6 ] Version xcdat : 0.7.1 GitHub branch : feature/416-custom-season-span

More details about time variables in the file

>>> ds.time_counter
<xarray.DataArray 'time_counter' (time_counter: 120)> Size: 960B
array([cftime.DatetimeGregorian(2000, 1, 16, 12, 0, 0, 0, has_year_zero=False),
....
     dtype=object)
Coordinates:
    time_centered  (time_counter) object 960B ...
  * time_counter   (time_counter) object 960B 2000-01-16 12:00:00 ... 2009-12...
Attributes:
    axis:           T
    standard_name:  time
    long_name:      Time axis
    time_origin:    2000-01-01 00:00:00
    bounds:         time_counter_bnds
    _ChunkSizes:    1
    
>>> ds.time_centered
<xarray.DataArray 'time_centered' (time_counter: 120)> Size: 960B
[120 values with dtype=object]
Coordinates:
    time_centered  (time_counter) object 960B ...
  * time_counter   (time_counter) object 960B 2000-01-16 12:00:00 ... 2009-12...
Attributes:
    standard_name:  time
    long_name:      Time axis
    time_origin:    2000-01-01 00:00:00
    bounds:         time_centered_bounds
    _ChunkSizes:    1