Repository: xCDAT/xcdat

Issue:
[Feature]: custom seasons that span calendar years
### Is your feature request related to a problem?

Following #393, it seems that it would be useful to expand [custom seasons](https://xcdat.readthedocs.io/en/stable/generated/xarray.Dataset.temporal.group_average.html) functionality across calendar years. Examples include: taking the water year from October to September, or taking a boreal winter average from December to March.

### Describe the solution you'd like

I envision the main change would be so that the order of months listed in [custom seasons](https://xcdat.readthedocs.io/en/stable/generated/xarray.Dataset.temporal.group_average.html) **does** matter, and could span across calendar years.

For example, this configuation would create Apr-Nov and Dec-March averages, the latter which extends into the following year:
```
custom_seasons = [
    ["Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov"],
    ["Dec", "Jan", "Feb", "Mar"],
]
```
Associated with this change I think it would also make sense to generalize the [season_config](https://xcdat.readthedocs.io/en/stable/generated/xarray.Dataset.temporal.group_average.html) parameters beyond the existing options for DJF, e.g. the flag `“drop_incomplete_djf”` could become something like `"drop_incomplete_season."` In addition to seasons that cross the calendar year, this could apply to datasets that end in the middle of a season (for example, the last Apr-June season for a dataset that ends in May).

### Describe alternatives you've considered

_No response_

### Additional context

As a possible aside, I was also wondering if it is necessary to keep the requirement to include all 12 months in the `custom_seasons` list. I imagine that often users are only interested in one season.

Discussion:
tomvothecoder: Example result of `_drop_incomplete_seasons()`:

```python

# Before dropping
# -----------------
# 2000-1, 2000-2, and 2001-12 months in incomplete "DJF" seasons" so they are dropped
ds.time
<xarray.DataArray 'time' (time: 15)>
array([cftime.DatetimeGregorian(2000, 1, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 2, 15, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 3, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 4, 16, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 5, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 6, 16, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 7, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 8, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 9, 16, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 10, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 11, 16, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 12, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2001, 1, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2001, 2, 15, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2001, 12, 16, 12, 0, 0, 0, has_year_zero=False)],
      dtype=object)
Coordinates:
  * time     (time) object 2000-01-16 12:00:00 ... 2001-12-16 12:00:00
Attributes:
    axis:           T
    long_name:      time
    standard_name:  time
    bounds:         time_bnds

# After dropping
# -----------------
ds_new.time
<xarray.DataArray 'time' (time: 12)>
array([cftime.DatetimeGregorian(2000, 3, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 4, 16, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 5, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 6, 16, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 7, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 8, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 9, 16, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 10, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 11, 16, 0, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2000, 12, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2001, 1, 16, 12, 0, 0, 0, has_year_zero=False),
       cftime.DatetimeGregorian(2001, 2, 15, 0, 0, 0, 0, has_year_zero=False)],
      dtype=object)
Coordinates:
  * time     (time) object 2000-03-16 12:00:00 ... 2001-02-15 00:00:00
Attributes:
    axis:           T
    long_name:      time
    standard_name:  time
    bounds:         time_bnds

```
tomvothecoder: Hey @lee1043, this PR seemed to be mostly done when I stopped working on it last year. I just had to fix a few things and update the tests. 

Would you like to check out this branch to test it out on real data? Also a code review would be appreciated.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/xCDAT/xcdat/pull/423?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) Report
Attention: Patch coverage is `98.41270%` with `1 lines` in your changes are missing coverage. Please review.
> Project coverage is 99.93%. Comparing base [(`472c04b`)](https://app.codecov.io/gh/xCDAT/xcdat/commit/472c04b36aa89eb7f9d6d49cff0b0851dfb8edce?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) to head [(`74c136c`)](https://app.codecov.io/gh/xCDAT/xcdat/pull/423?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).

| [Files](https://app.codecov.io/gh/xCDAT/xcdat/pull/423?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) | Patch % | Lines |
|---|---|---|
| [xcdat/temporal.py](https://app.codecov.io/gh/xCDAT/xcdat/pull/423?src=pr&el=tree&filepath=xcdat%2Ftemporal.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT#diff-eGNkYXQvdGVtcG9yYWwucHk=) | 98.41% | [1 Missing :warning: ](https://app.codecov.io/gh/xCDAT/xcdat/pull/423?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) |

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##              main     #423      +/-   ##
===========================================
- Coverage   100.00%   99.93%   -0.07%     
===========================================
  Files           15       15              
  Lines         1521     1552      +31     
===========================================
+ Hits          1521     1551      +30     
- Misses           0        1       +1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/xCDAT/xcdat/pull/423?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).

lee1043: @tomvothecoder sure, I will test it out and review. Thank you for the update!
lee1043: @tomvothecoder Can this be considered for v0.7.0?
tomvothecoder: > @tomvothecoder Can this be considered for v0.7.0?

This PR still needs thorough review before I'm confident in releasing it. I'll probably tag Steve at some point.

The release after v0.7.0 is more realistic and reasonable. We can always initiate a new release for this feature whenever it is merged.
lee1043: > > @tomvothecoder Can this be considered for v0.7.0?
> 
> This PR still needs thorough review before I'm confident in merging it. I'll probably tag Steve at some point.
> 
> The release after v0.7.0 is more realistic and reasonable. We can always initiate a new release for this feature whenever it is merged.

@tomvothecoder no problem. Thank you for consideration. 
lee1043: @tomvothecoder it looks like when custom season go beyond calendar year (Nov, Dec, Jan) there is error as follows.

```python
import os
import xcdat as xc

input_data = os.path.join(
    "/p/css03/esgf_publish/CMIP6/CMIP/AWI/AWI-CM-1-1-MR/historical/r1i1p1f1/Amon/psl/gn/v20181218/",
    "psl_Amon_AWI-CM-1-1-MR_historical_r1i1p1f1_gn_201301-201312.nc")

ds = xc.open_mfdataset(input_data)

# Example of custom seasons in a three month format:
custom_seasons = [
    ['Dec', 'Jan'],
]

season_config = {'custom_seasons': custom_seasons, 'dec_mode': 'DJF', 'drop_incomplete_djf': True}

ds.temporal.group_average("psl", "season", season_config=season_config)
```


```
CPU times: user 448 ms, sys: 32.8 ms, total: 481 ms
Wall time: 471 ms
xarray.Dataset
Dimensions:
lat: 192bnds: 2lon: 384time: 4
Coordinates:
lat
(lat)
float64
-89.28 -88.36 ... 88.36 89.28
lon
(lon)
float64
0.0 0.9375 1.875 ... 358.1 359.1
time
(time)
object
2013-02-01 00:00:00 ... 2013-11-...
Data variables:
lat_bnds
(lat, bnds)
float64
dask.array<chunksize=(192, 2), meta=np.ndarray>
lon_bnds
(lon, bnds)
float64
dask.array<chunksize=(384, 2), meta=np.ndarray>
psl
(time, lat, lon)
float64
dask.array<chunksize=(1, 192, 384), meta=np.ndarray>
Indexes: (3)
Attributes: (44)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xarray/core/dataarray.py:862, in DataArray._getitem_coord(self, key)
    861 try:
--> 862     var = self._coords[key]
    863 except KeyError:

KeyError: 'time.year'

During handling of the above exception, another exception occurred:

AttributeError                            Traceback (most recent call last)
<timed exec> in ?()

~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, data_var, freq, weighted, keep_weights, season_config)
    400         }
    401         """
    402         self._set_data_var_attrs(data_var)
    403 
--> 404         return self._averager(
    405             data_var,
    406             "group_average",
    407             freq,

~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, data_var, mode, freq, weighted, keep_weights, reference_period, season_config)
    870 
    871         if self._mode == "average":
    872             dv_avg = self._average(ds, data_var)
    873         elif self._mode in ["group_average", "climatology", "departures"]:
--> 874             dv_avg = self._group_average(ds, data_var)
    875 
    876         # The original time dimension is dropped from the dataset because
    877         # it becomes obsolete after the data variable is averaged. When the

~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, ds, data_var)
   1295         dv = _get_data_var(ds, data_var)
   1296 
   1297         # Label the time coordinates for grouping weights and the data variable
   1298         # values.
-> 1299         self._labeled_time = self._label_time_coords(dv[self.dim])
   1300 
   1301         if self._weighted:
   1302             time_bounds = ds.bounds.get_bounds("T", var_key=data_var)

~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, time_coords)
   1470         >>>       dtype='datetime64[ns]')
   1471         >>> Coordinates:
   1472         >>> * time     (time) datetime64[ns] 2000-01-01T00:00:00 ... 2000-04-01T00:00:00
   1473         """
-> 1474         df_dt_components: pd.DataFrame = self._get_df_dt_components(
   1475             time_coords, drop_obsolete_cols=True
   1476         )
   1477         dt_objects = self._convert_df_to_dt(df_dt_components)

~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, time_coords, drop_obsolete_cols)
   1534 
   1535         # Use the TIME_GROUPS dictionary to determine which components
   1536         # are needed to form the labeled time coordinates.
   1537         for component in TIME_GROUPS[self._mode][self._freq]:
-> 1538             df[component] = time_coords[f"{self.dim}.{component}"].values
   1539 
   1540         # The season frequency requires additional datetime components for
   1541         # processing, which are later removed before time coordinates are

~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xarray/core/dataarray.py in ?(self, key)
    869     def __getitem__(self, key: Any) -> Self:
    870         if isinstance(key, str):
--> 871             return self._getitem_coord(key)
    872         else:
    873             # xarray-style array indexing
    874             return self.isel(indexers=self._item_key_to_dict(key))

~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xarray/core/dataarray.py in ?(self, key)
    861         try:
    862             var = self._coords[key]
    863         except KeyError:
    864             dim_sizes = dict(zip(self.dims, self.shape))
--> 865             _, key, var = _get_virtual_variable(self._coords, key, dim_sizes)
    866 
    867         return self._replace_maybe_drop_dims(var, name=key)

~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xarray/core/dataset.py in ?(variables, key, dim_sizes)
    212     if _contains_datetime_like_objects(ref_var):
    213         ref_var = DataArray(ref_var)
    214         data = getattr(ref_var.dt, var_name).data
    215     else:
--> 216         data = getattr(ref_var, var_name).data
    217     virtual_var = Variable(ref_var.dims, data)
    218 
    219     return ref_name, var_name, virtual_var

AttributeError: 'IndexVariable' object has no attribute 'year'
CPU times: user 240 ms, sys: 13.2 ms, total: 253 ms
Wall time: 249 ms
xarray.Dataset
Dimensions:
lat: 192bnds: 2lon: 384time: 1
Coordinates:
lat
(lat)
float64
-89.28 -88.36 ... 88.36 89.28
lon
(lon)
float64
0.0 0.9375 1.875 ... 358.1 359.1
time
(time)
object
2013-02-01 00:00:00
Data variables:
lat_bnds
(lat, bnds)
float64
dask.array<chunksize=(192, 2), meta=np.ndarray>
lon_bnds
(lon, bnds)
float64
dask.array<chunksize=(384, 2), meta=np.ndarray>
psl
(time, lat, lon)
float64
dask.array<chunksize=(1, 192, 384), meta=np.ndarray>
Indexes: (3)
Attributes: (44)
CPU times: user 106 ms, sys: 6.96 ms, total: 113 ms
Wall time: 110 ms
xarray.Dataset
Dimensions:
lat: 192bnds: 2lon: 384time: 1
Coordinates:
lat
(lat)
float64
-89.28 -88.36 ... 88.36 89.28
lon
(lon)
float64
0.0 0.9375 1.875 ... 358.1 359.1
time
(time)
object
2013-02-01 00:00:00
Data variables:
lat_bnds
(lat, bnds)
float64
dask.array<chunksize=(192, 2), meta=np.ndarray>
lon_bnds
(lon, bnds)
float64
dask.array<chunksize=(384, 2), meta=np.ndarray>
psl
(time, lat, lon)
float64
dask.array<chunksize=(1, 192, 384), meta=np.ndarray>
Indexes: (3)
Attributes: (44)
```
tomvothecoder: @lee1043 Thanks for trying to this out and providing an example script! I'll debug the stack trace.
tomvothecoder: Hi @oliviermarti, thanks for your interest in this feature. I've had minimal time to work on it and plan to ramp up development soon. I don't have a set timeline on when it will be done, but I'm hoping within the next few months.
tomvothecoder: > @tomvothecoder it looks like when custom season go beyond calendar year (Nov, Dec, Jan) there is error as follows.
> 
> ```python
> import os
> import xcdat as xc
> 
> input_data = os.path.join(
>     "/p/css03/esgf_publish/CMIP6/CMIP/AWI/AWI-CM-1-1-MR/historical/r1i1p1f1/Amon/psl/gn/v20181218/",
>     "psl_Amon_AWI-CM-1-1-MR_historical_r1i1p1f1_gn_201301-201312.nc")
> 
> ds = xc.open_mfdataset(input_data)
> 
> # Example of custom seasons in a three month format:
> custom_seasons = [
>     ['Dec', 'Jan'],
> ]
> 
> season_config = {'custom_seasons': custom_seasons, 'dec_mode': 'DJF', 'drop_incomplete_djf': True}
> 
> ds.temporal.group_average("psl", "season", season_config=season_config)
> ```
> 
> ```
> CPU times: user 448 ms, sys: 32.8 ms, total: 481 ms
> Wall time: 471 ms
> xarray.Dataset
> Dimensions:
> lat: 192bnds: 2lon: 384time: 4
> Coordinates:
> lat
> (lat)
> float64
> -89.28 -88.36 ... 88.36 89.28
> lon
> (lon)
> float64
> 0.0 0.9375 1.875 ... 358.1 359.1
> time
> (time)
> object
> 2013-02-01 00:00:00 ... 2013-11-...
> Data variables:
> lat_bnds
> (lat, bnds)
> float64
> dask.array<chunksize=(192, 2), meta=np.ndarray>
> lon_bnds
> (lon, bnds)
> float64
> dask.array<chunksize=(384, 2), meta=np.ndarray>
> psl
> (time, lat, lon)
> float64
> dask.array<chunksize=(1, 192, 384), meta=np.ndarray>
> Indexes: (3)
> Attributes: (44)
> ---------------------------------------------------------------------------
> KeyError                                  Traceback (most recent call last)
> File ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xarray/core/dataarray.py:862, in DataArray._getitem_coord(self, key)
>     861 try:
> --> 862     var = self._coords[key]
>     863 except KeyError:
> 
> KeyError: 'time.year'
> 
> During handling of the above exception, another exception occurred:
> 
> AttributeError                            Traceback (most recent call last)
> <timed exec> in ?()
> 
> ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, data_var, freq, weighted, keep_weights, season_config)
>     400         }
>     401         """
>     402         self._set_data_var_attrs(data_var)
>     403 
> --> 404         return self._averager(
>     405             data_var,
>     406             "group_average",
>     407             freq,
> 
> ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, data_var, mode, freq, weighted, keep_weights, reference_period, season_config)
>     870 
>     871         if self._mode == "average":
>     872             dv_avg = self._average(ds, data_var)
>     873         elif self._mode in ["group_average", "climatology", "departures"]:
> --> 874             dv_avg = self._group_average(ds, data_var)
>     875 
>     876         # The original time dimension is dropped from the dataset because
>     877         # it becomes obsolete after the data variable is averaged. When the
> 
> ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, ds, data_var)
>    1295         dv = _get_data_var(ds, data_var)
>    1296 
>    1297         # Label the time coordinates for grouping weights and the data variable
>    1298         # values.
> -> 1299         self._labeled_time = self._label_time_coords(dv[self.dim])
>    1300 
>    1301         if self._weighted:
>    1302             time_bounds = ds.bounds.get_bounds("T", var_key=data_var)
> 
> ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, time_coords)
>    1470         >>>       dtype='datetime64[ns]')
>    1471         >>> Coordinates:
>    1472         >>> * time     (time) datetime64[ns] 2000-01-01T00:00:00 ... 2000-04-01T00:00:00
>    1473         """
> -> 1474         df_dt_components: pd.DataFrame = self._get_df_dt_components(
>    1475             time_coords, drop_obsolete_cols=True
>    1476         )
>    1477         dt_objects = self._convert_df_to_dt(df_dt_components)
> 
> ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xcdat/temporal.py in ?(self, time_coords, drop_obsolete_cols)
>    1534 
>    1535         # Use the TIME_GROUPS dictionary to determine which components
>    1536         # are needed to form the labeled time coordinates.
>    1537         for component in TIME_GROUPS[self._mode][self._freq]:
> -> 1538             df[component] = time_coords[f"{self.dim}.{component}"].values
>    1539 
>    1540         # The season frequency requires additional datetime components for
>    1541         # processing, which are later removed before time coordinates are
> 
> ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xarray/core/dataarray.py in ?(self, key)
>     869     def __getitem__(self, key: Any) -> Self:
>     870         if isinstance(key, str):
> --> 871             return self._getitem_coord(key)
>     872         else:
>     873             # xarray-style array indexing
>     874             return self.isel(indexers=self._item_key_to_dict(key))
> 
> ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xarray/core/dataarray.py in ?(self, key)
>     861         try:
>     862             var = self._coords[key]
>     863         except KeyError:
>     864             dim_sizes = dict(zip(self.dims, self.shape))
> --> 865             _, key, var = _get_virtual_variable(self._coords, key, dim_sizes)
>     866 
>     867         return self._replace_maybe_drop_dims(var, name=key)
> 
> ~/.conda/envs/xcdat_dev_20240401/lib/python3.11/site-packages/xarray/core/dataset.py in ?(variables, key, dim_sizes)
>     212     if _contains_datetime_like_objects(ref_var):
>     213         ref_var = DataArray(ref_var)
>     214         data = getattr(ref_var.dt, var_name).data
>     215     else:
> --> 216         data = getattr(ref_var, var_name).data
>     217     virtual_var = Variable(ref_var.dims, data)
>     218 
>     219     return ref_name, var_name, virtual_var
> 
> AttributeError: 'IndexVariable' object has no attribute 'year'
> CPU times: user 240 ms, sys: 13.2 ms, total: 253 ms
> Wall time: 249 ms
> xarray.Dataset
> Dimensions:
> lat: 192bnds: 2lon: 384time: 1
> Coordinates:
> lat
> (lat)
> float64
> -89.28 -88.36 ... 88.36 89.28
> lon
> (lon)
> float64
> 0.0 0.9375 1.875 ... 358.1 359.1
> time
> (time)
> object
> 2013-02-01 00:00:00
> Data variables:
> lat_bnds
> (lat, bnds)
> float64
> dask.array<chunksize=(192, 2), meta=np.ndarray>
> lon_bnds
> (lon, bnds)
> float64
> dask.array<chunksize=(384, 2), meta=np.ndarray>
> psl
> (time, lat, lon)
> float64
> dask.array<chunksize=(1, 192, 384), meta=np.ndarray>
> Indexes: (3)
> Attributes: (44)
> CPU times: user 106 ms, sys: 6.96 ms, total: 113 ms
> Wall time: 110 ms
> xarray.Dataset
> Dimensions:
> lat: 192bnds: 2lon: 384time: 1
> Coordinates:
> lat
> (lat)
> float64
> -89.28 -88.36 ... 88.36 89.28
> lon
> (lon)
> float64
> 0.0 0.9375 1.875 ... 358.1 359.1
> time
> (time)
> object
> 2013-02-01 00:00:00
> Data variables:
> lat_bnds
> (lat, bnds)
> float64
> dask.array<chunksize=(192, 2), meta=np.ndarray>
> lon_bnds
> (lon, bnds)
> float64
> dask.array<chunksize=(384, 2), meta=np.ndarray>
> psl
> (time, lat, lon)
> float64
> dask.array<chunksize=(1, 192, 384), meta=np.ndarray>
> Indexes: (3)
> Attributes: (44)
> ```


### The issue is that this dataset only has a single year of data and does not have a complete `["Dec", "Jan"]` season with `drop_incomplete_djf=True`. 
```
<xarray.DataArray 'time' (time: 12)> Size: 96B
array([cftime.DatetimeProlepticGregorian(2013, 1, 16, 12, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 2, 15, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 3, 16, 12, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 4, 16, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 5, 16, 12, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 6, 16, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 7, 16, 12, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 8, 16, 12, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 9, 16, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 10, 16, 12, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 11, 16, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeProlepticGregorian(2013, 12, 16, 12, 0, 0, 0, has_year_zero=True)],
      dtype=object)
Coordinates:
  * time     (time) object 96B 2013-01-16 12:00:00 ... 2013-12-16 12:00:00
Attributes:
    bounds:         time_bnds
    axis:           T
    long_name:      time
    standard_name:  time
```

* `dec_mode="DJF"` shifts Decembers to the next year (for proper DJF grouping instead of JFD)
* `drop_incomplete_djf=True` (now `drop_incomplete_seasons=True`) drops all of the time coordinates since there are no complete seasons using `["Dec", "Jan"]`.
 
### Resulting time coordinates after subsetting for Dec and Jan coordinates and applying above options:
```python
   year  season  month  expected_months  actual_months
0  2013  DecJan      1                2              1
1  2014  DecJan     12                2              1
```

Since expected and actual months don't align, the coordinates are dropped. This results in no time coordinates remaining leading to `> AttributeError: 'IndexVariable' object has no attribute 'year'`

Workaround -- specify `drop_incomplete_seasons=False`

I also added a clear error message no complete seasons are found ([commit](https://github.com/xCDAT/xcdat/pull/423/commits/f6858ff8dd2d7e17114b66e85642bfdbc22f87bc)).


tomvothecoder: After thorough testing, I actually think this PR is close to complete.
tomvothecoder: Hey @oliviermarti, @arfriedman, and @DamienIrving, would any of you be interested in checking out this branch to test the custom season functionality? If so, please check out our [contributing guide](https://xcdat.readthedocs.io/en/stable/contributing.html#version-control-git-and-github). Thanks!
oliviermarti: @tomvothecoder : this version seems perfect to me :-)

(I've got the branch feature/416-custom-season-span : correct ? )

Thanks !!

I join my very simple test case, which runs well : 
```python
import os, sys
import numpy as np, xarray as xr, xcdat as xc

# Version check
print ( f'{sys.version = }' )

print ( f'{np.__version__ = }' )
print ( f'{xr.__version__ = }' )
print ( f'{xc.__version__ = }' )
print ( f'{xc.__file__ = }' )

# Creates time axis
nt=36

time = cftime.num2date ( np.arange(36,dtype=float)*30.0+15, units="day since 2000-01-01", calendar='360_day', only_use_cftime_datetimes=True,
                only_use_python_datetimes=False, has_year_zero=None)

time = xr.DataArray (time, dims=('time',), coords=(time,) )

# Creates two tests variables
nmonth = np.arange (nt,dtype=float) + 1
nmod   = np.mod (np.arange (nt,dtype=float), 12) +1

nmonth = xr.DataArray (nmonth, dims=('time',), coords=(time,) )
nmod   = xr.DataArray (nmod  , dims=('time',), coords=(time,) )

# Creates dataset
dd = xr.Dataset ( {'nmonth': nmonth, 'nmod':nmod})

# Writes dataset and open it with xcdat
if os.path.exists ('test.nc') : os.remove ('test.nc')
dd.to_netcdf ( 'test.nc', mode='w')
dx = xc.open_dataset ('test.nc',  use_cftime=True, decode_times=True).bounds.add_missing_bounds()

print ( 'Check variables')
print (dx.nmonth.values)
print (dx.nmod.values)

print ( 'Year means are correct')
print (dx.temporal.group_average("nmonth", freq="year", weighted=True)['nmonth'].values)
print (dx.temporal.group_average("nmod"  , freq="year", weighted=True)['nmod'  ].values)


print ( ' Standard seasonnal means are correct')
print (dx.temporal.group_average("nmonth", freq="season", weighted=True,
                                season_config={'drop_incomplete_djf':False, 'dec_mode': 'DJF',})['nmonth'].values)
print (dx.temporal.group_average("nmod"  , freq="season", weighted=True,
                                season_config={'drop_incomplete_djf':False, 'dec_mode': 'DJF',})['nmod'  ].values)

print (dx.temporal.group_average("nmonth", freq="season", weighted=True, 
                                 season_config={'drop_incomplete_djf':True, 'dec_mode': 'DJF',})['nmonth'].values)
print (dx.temporal.group_average("nmod"  , freq="season", weighted=True, 
                                 season_config={'drop_incomplete_djf':True, 'dec_mode': 'DJF',})['nmod'  ].values)


print ( '3 months custom seasons are correct' )

custom_seasons = [  ["Dec", "Jan", "Feb"], ]

print (dx.temporal.group_average("nmonth", freq="season", 
                                 season_config={'drop_incomplete_djf':False, 'custom_seasons':custom_seasons, 'dec_mode': 'DJF',}, 
                                 weighted=True)['nmonth'].values)
print (dx.temporal.group_average("nmod"  , freq="season", 
                                 season_config={'drop_incomplete_djf':False, 'custom_seasons':custom_seasons, 'dec_mode': 'DJF',}, 
                                 weighted=True)['nmod'  ].values)

print (dx.temporal.group_average("nmonth", freq="season", 
                                 season_config={'drop_incomplete_djf':True, 'custom_seasons':custom_seasons, 'dec_mode': 'DJF',}, 
                                 weighted=True)['nmonth'].values)
print (dx.temporal.group_average("nmod"  , freq="season", 
                                 season_config={'drop_incomplete_djf':True, 'custom_seasons':custom_seasons, 'dec_mode': 'DJF',}, 
                                 weighted=True)['nmod'  ].values)


print ( '4 months custom seasons are correct' )

custom_seasons = [  ["Dec", "Jan", "Feb", "Mar"], ]

print (dx.temporal.group_average("nmonth", freq="season", 
                                 season_config={'drop_incomplete_djf':False, 'custom_seasons':custom_seasons, 'dec_mode': 'DJF',}, 
                                 weighted=True)['nmonth'].values)
print (dx.temporal.group_average("nmod"  , freq="season", 
                                 season_config={'drop_incomplete_djf':False, 'custom_seasons':custom_seasons, 'dec_mode': 'DJF',}, 
                                 weighted=True)['nmod'  ].values)


print (dx.temporal.group_average("nmonth", freq="season", 
                                 season_config={'drop_incomplete_djf':True, 'custom_seasons':custom_seasons, 'dec_mode': 'DJF',}, 
                                 weighted=True)['nmonth'].values)
print (dx.temporal.group_average("nmod"  , freq="season", 
                                 season_config={'drop_incomplete_djf':True, 'custom_seasons':custom_seasons, 'dec_mode': 'DJF',}, 
                                 weighted=True)['nmod'  ].values)
```

Results : 

```
sys.version = '3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:45:13) [Clang 16.0.6 ]'
np.__version__ = '1.26.4'
xr.__version__ = '2024.7.0'
xc.__version__ = '0.7.1'
xc.__file__ = '/Users/marti/Unix/XCDAT/xcdat/__init__.py'
Check variables
[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16. 17. 18.
 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36.]
[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12.  1.  2.  3.  4.  5.  6.
  7.  8.  9. 10. 11. 12.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12.]
Year means are correct
[ 6.5 18.5 30.5]
[6.5 6.5 6.5]
 Standard seasonnal means are correct
[ 1.5  4.   7.  10.  13.  16.  19.  22.  25.  28.  31.  34.  36. ]
[ 1.5  4.   7.  10.   5.   4.   7.  10.   5.   4.   7.  10.  12. ]
[ 4.  7. 10. 13. 16. 19. 22. 25. 28. 31. 34.]
[ 4.  7. 10.  5.  4.  7. 10.  5.  4.  7. 10.]
3 months custom seasons are correct
[ 1.5 13.  25.  36. ]
[ 1.5  5.   5.  12. ]
[13. 25.]
[5. 5.]
4 months custom seasons are correct
[ 2.  13.5 25.5 36. ]
[ 2.   4.5  4.5 12. ]
[13.5 25.5]
[4.5 4.5]
```

tomvothecoder: > @tomvothecoder : this version seems perfect to me :-)
> 
> (I've got the branch feature/416-custom-season-span : correct ? )
> 
> Thanks !!
> 
> I join my very simple test case, which runs well :

Yes that is the correct branch. I'm happy to see your test cases pass! If you're interested, it'd be helpful to test on some real data too. 
oliviermarti: > > @tomvothecoder : this version seems perfect to me :-)
> > (I've got the branch feature/416-custom-season-span : correct ? )
> > Thanks !!
> > I join my very simple test case, which runs well :
> 
> Yes that is the correct branch. I'm happy to see your test cases pass! If you're interested, it'd be helpful to test on some real data too.

Tom,

On real data, version 0.7.1 returns a result, but a wrong result. It's hards to have an easy check on real data. Hence this academic test case.

I've made a fews test on real data, from ESGF, and raw data from my IPSL model. It seems OK.

I've got problems with the IPSL model, but it's an another story outside the present scope : IPSL raw outputs have one time axis 'time_counter', but two time variables : `time_counter` and `time_centered`. `add_bounds`works nicely on both, but in the end I got : 
`
 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
`
I have to rename `time_counter` to `time` variable to make it work.

IPSL outputs processed and available on ESGF are OK.
tomvothecoder: > On real data, version 0.7.1 returns a result, but a wrong result. It's hards to have an easy check on real data. Hence this academic test case.
> 
> I've made a fews test on real data, from ESGF, and raw data from my IPSL model. It seems OK.

Thank you. @lee1043 is also going to test on CMIP data. 

> I've got problems with the IPSL model, but it's an another story outside the present scope : IPSL raw outputs have one time axis 'time_counter', but two time variables : `time_counter` and `time_centered`. `add_bounds`works nicely on both, but in the end I got : `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` I have to rename `time_counter` to `time` variable to make it work.
> 
> IPSL outputs processed and available on ESGF are OK.

I pushed [`388f4da` (#423)](https://github.com/xCDAT/xcdat/pull/423/commits/388f4dac82a1330197228ff02e1551438703a5b2) to fix #684. Your example code in that issue now works. 

arfriedman: > Hey @oliviermarti, @arfriedman, and @DamienIrving, would any of you be interested in checking out this branch to test the custom season functionality? If so, please check out our [contributing guide](https://xcdat.readthedocs.io/en/stable/contributing.html#version-control-git-and-github). Thanks!

Thank you @tomvothecoder! -- I will take a look.
dcherian: Here's an upstream version: https://github.com/pydata/xarray/pull/9524 . I could use some help testing it out.
tomvothecoder: I want to prioritize this PR for the next release (v0.8.0), which I'm aiming by the next E3SM Unified release scheduled for 2/15/24. 

We can consider refactoring the implementation with the new feature from Xarray PR [#9524](https://github.com/pydata/xarray/pull/9524) if it makes sense to (e.g., simplifies API and/or underlying implementation, same expected behavior, performance gains, etc.). 


tomvothecoder: @lee1043 Have you tried testing this PR out recently?
lee1043: @tomvothecoder I will give another round of test for this PR in this week. Sorry, it slipped from my list and thank you for the reminder!
lee1043: I used 2 years of monthly time series data to test this PR. From the test results below, I believe the code is working as expected. 

## 1. First sanity check: compare xcdat default vs custom season

```python
import os
import xarray as xr
import xcdat as xc
import matplotlib.pyplot as plt

# input data is from 
# "/p/css03/esgf_publish/CMIP6/CMIP/AWI/AWI-CM-1-1-MR/historical/r1i1p1f1/Amon/psl/gn/v20181218/"

input_data = [
    "psl_Amon_AWI-CM-1-1-MR_historical_r1i1p1f1_gn_201201-201212.nc",
    "psl_Amon_AWI-CM-1-1-MR_historical_r1i1p1f1_gn_201301-201312.nc"
]

ds = xc.open_mfdataset(input_data)

# Sanity check

def compare_xcdat_default_vs_custom(custom_season, custom_drop_incomplete_seasons=True, ax=None, title=None):
    
    if custom_season == "DJF":
        custom_seasons = [['Dec', 'Jan', 'Feb']]
        index = 0
    elif custom_season == "MAM":
        custom_seasons = [['Mar', 'Apr', 'May']]
        index = 1
    elif custom_season == "JJA":
        custom_seasons = [['Jun', 'Jul', 'Aug']]
        index = 2
    elif custom_season == "SON":
        custom_seasons = [['Sep', 'Oct', 'Nov']]
        index = 3
        
    # season from default method
    ds_default_season = ds.temporal.climatology(
        "psl",
        freq="season",
        weighted=True,
        season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
    )

    da_default_season = ds_default_season['psl'][index]
    #print(custom_season, 'default', ds_default_season.sizes, da_default_season.shape)
    
    # season from custom season
    season_config = {'custom_seasons': custom_seasons, 'dec_mode': 'DJF', 'drop_incomplete_seasons': custom_drop_incomplete_seasons}
    ds_custom_season = ds.temporal.group_average("psl", "season", season_config=season_config)
    # get climatology
    if len(ds_custom_season.time) > 1:
        ds_custom_season = ds_custom_season.bounds.add_missing_bounds()
        ds_custom_season = ds_custom_season.temporal.average("psl")
    da_custom_season = ds_custom_season["psl"].squeeze()
    
    #print(custom_season, 'custom', ds_custom_season.sizes, da_custom_season.shape)
    
    # plot
    (da_default_season - da_custom_season).plot(ax=ax)
    
    if ax is not None and title is not None:
        ax.set_title(title)

fig, ax = plt.subplots(2, 4, figsize=(12, 5))

seasons = ["DJF", "MAM", "JJA", "SON"]

for i, season in enumerate(seasons):
    compare_xcdat_default_vs_custom(season, custom_drop_incomplete_seasons=True, ax=ax[0,i], title=season)
    compare_xcdat_default_vs_custom(season, custom_drop_incomplete_seasons=False, ax=ax[1,i], title=season)

fig.suptitle("xcdat default - custom \n upper: custom_drop_incomplete_seasons=True \n upper: custom_drop_incomplete_seasons=False")

fig.tight_layout()
```
![compare_xcdat_default_vs_custom_season](https://github.com/user-attachments/assets/0249fdc8-c397-4b40-8d0f-8b592d2be7fb)

The result is consistent to that is expected. DJF is one time step with `drop_incomplete_seasons=True` so no difference between default and custom approach (upper left). It became different when incomplete season is not dropped from custom, while it was dropped from default (lower left). Other seasons are not very much influenced by the parameter. so noise like differences are shown.

## 2. Compare to other tools
```python
custom_season = "JJAS"
custom_seasons = [['Jun', 'Jul', 'Aug', 'Sep']]

#
# xcdat
#
ds = xc.open_mfdataset(input_data)

season_config = {'custom_seasons': custom_seasons, 'dec_mode': 'DJF', 'drop_incomplete_seasons': True}

ds_cs_xcdat = ds.temporal.group_average("psl", "season", season_config=season_config)
# get climatology
if len(ds_cs_xcdat.time) > 1:
    ds_cs_xcdat = ds_cs_xcdat.bounds.add_missing_bounds()
    ds_cs_xcdat = ds_cs_xcdat.temporal.average("psl")
da_cs_xcdat = ds_cs_xcdat["psl"].squeeze()

#
# CDAT
# 

# "cdscan -x test.xml *.nc" conducted as pre-process

import cdms2
import cdutil

f = cdms2.open('test.xml')
d = f('psl')
custom_season_class = cdutil.times.Seasons(custom_season)
d_cs = custom_season_class.climatology(d)

# 
# PMP (xcdat and xarray method)
# 
from pcmdi_metrics.utils import custom_season_average

# xcdat method
ds_cs_pmp_xcdat_yearly = custom_season_average(ds, "psl", season=custom_season, method="xcdat")
ds_cs_pmp_xcdat_yearly = ds_cs_pmp_xcdat_yearly.bounds.add_missing_bounds()
ds_cs_pmp_xcdat = ds_cs_pmp_xcdat_yearly.temporal.average("psl")
da_cs_pmp_xcdat = ds_cs_pmp_xcdat['psl']

# xarray method
ds_cs_pmp_xarray_yearly = custom_season_average(ds, "psl", season=custom_season, method="xarray")
ds_cs_pmp_xarray = ds_cs_pmp_xarray_yearly.mean(dim=["year"])
da_cs_pmp_xarray = ds_cs_pmp_xarray['psl']

# 
# Compare results
# 
fig, ax = plt.subplots(2, 4, figsize=(12, 5))

da_cs_xcdat.plot(ax=ax[0,0])
da_cs_cdat.plot(ax=ax[0,1])
da_cs_pmp_xcdat.plot(ax=ax[0,2])
da_cs_pmp_xarray.plot(ax=ax[0,3])


(da_cs_xcdat.to_numpy() - da_cs_cdat).plot(ax=ax[1,0])
(da_cs_cdat - da_cs_pmp_xcdat.to_numpy()).plot(ax=ax[1,1])
(da_cs_xcdat - da_cs_pmp_xcdat).plot(ax=ax[1,2])
(da_cs_pmp_xcdat - da_cs_pmp_xarray).plot(ax=ax[1,3])

ax[0,0].set_title("xcdat")
ax[0,1].set_title("cdat")
ax[0,2].set_title("pmp (xc)")
ax[0,3].set_title("pmp (xr)")

ax[1,0].set_title("xcdat - cdat")
ax[1,1].set_title("cdat - pmp (xc)")
ax[1,2].set_title("xcdat - pmp (xc)")
ax[1,3].set_title("pmp (xc) - pmp (xr)")

fig.suptitle(custom_season)

fig.tight_layout()
```

![compare_xcdat_cdat_pmp_JJAS](https://github.com/user-attachments/assets/dc0c26d9-db87-4ebd-9f39-17aac47d2b32)

No difference between xcdat, cdat, and pmp (xcdat) results, while pmp (xarray) is inconsistent maybe because it does not consider propor temporal weighting. 
lee1043: @tomvothecoder if you plan to switch `drop_incomplete_djf` to `drop_incomplete_seasons` in this PR, do you plan to apply that to other functions as well? That would makes sense if there is a time series starting from July ending in June after years, and want to get JJA temporal average. But this can be a breaking change if `drop_incomplete_djf` is being used in one's code and they update xcdat version. I wonder if the following could be a reasonable transition plan: mapping `drop_incomplete_djf` to `drop_incomplete_seasons` during the transition, and print a deprecation warning sign if `drop_incomplete_djf` is being used.
tomvothecoder: > @tomvothecoder if you plan to switch `drop_incomplete_djf` to `drop_incomplete_seasons` in this PR, do you plan to apply that to other functions as well? That would makes sense if there is a time series starting from July ending in June after years, and want to get JJA temporal average. But this can be a breaking change if `drop_incomplete_djf` is being used in one's code and they update xcdat version.

The `drop_incomplete_seasons` functionality applies to all of the group averaging APIs (`group_average`, `climatology`, and `departures`).


>  I wonder if the following could be a reasonable transition plan: mapping `drop_incomplete_djf` to `drop_incomplete_seasons` during the transition, and print a deprecation warning sign if `drop_incomplete_djf` is being used.

I commented in this [thread](https://github.com/xCDAT/xcdat/pull/423#discussion_r1786701423) above that users will who use `drop_incomplete_djf` will receive a `DeprecationWarning` to use `drop_incomplete_seasons` instead.

In this PR right now, `drop_incomplete_djf` is mapped to `drop_incomplete_seasons`. However, I don't think we should map these configs because they behave differently. `drop_incomplete_djf` removes the leading and trailing incomplete DJF seasons, while `drop_incomplete_seasons` will drop _all_ incomplete seasons. It is better to keep the behaviors of these configs separate to not break any existing code that depends on the current `drop_incomplete_djf` behavior. Users will have time to transition over to `drop_incomplete_seasons` before `drop_incomplete_djf` is deprecated.
tomvothecoder: From our meeting on 10/23:

- Near-term solution: move forward with this PR and merge it when ready
- Long-term solution: consider refactoring the underlying logic with Xarray [#9524](https://github.com/pydata/xarray/pull/9524) in the future
  - Scope out effort in refactoring xCDAT temporal backend to move away from Pandas to Xarray (need multi variable grouping)
  - Experiment with this Xarray PR now, see if it fits our requirements for custom seasons
tomvothecoder: @xCDAT/core-developers For the final time coordinate output, the current logic maps the custom season to its middle month. For example, year 2000 and "NDJFM" will be represented by `cftime(2000, 1, 1, 0...)` after seasonal group averaging.

**How useful would it be to keep the custom seasons as auxiliary coords?** I think it would be helpful to have since it explicitly shows what the custom seasons are. In Xarray, grouping by season will result in "season" coords.  For example,
```python
  <xarray.DataArray 'season' (season: 2)> Size: 32B
  array(['DJFM', 'AMJ'], dtype='<U4')
  Coordinates:
    * season   (season) <U4 32B 'DJFM' 'AMJ'
 ```

Alternatively, the user can just reference the attributes of the data variable to see the custom seasons.
```
{'operation': 'temporal_avg',
 'mode': 'group_average',
 'freq': 'season',
 'weighted': 'True',
 'drop_incomplete_seasons': 'True',
 'custom_seasons': ['NovDecJanFebMar']}
```

Let me know your thoughts.

tomvothecoder: In commit [`0b6852f` (#423)](https://github.com/xCDAT/xcdat/pull/423/commits/0b6852f0659030c9d4d763c4673416a95f8a4b3c), I refactored the logic for shifting months to use Xarray/NumPy instead of Pandas. We'll gradually shift away from the Pandas back-end for storing and manipulating Datetime components (used to group time coordinates) towards Xarray/Numpy in #217 since `xarray >=2024.09.0` now supports grouping by multiple variables.

I will do a final walk through at the next xCDAT meeting (11/20) before merging.
tomvothecoder: Thank for testing @lee1043. I will now merge this PR!

Patch:
diff --git a/xcdat/temporal.py b/xcdat/temporal.py
index ce44bbf..3367e06 100644
--- a/xcdat/temporal.py
+++ b/xcdat/temporal.py
@@ -1,5 +1,6 @@
 """Module containing temporal functions."""
 
+import warnings
 from datetime import datetime
 from itertools import chain
 from typing import Dict, List, Literal, Optional, Tuple, TypedDict, Union, get_args
@@ -66,8 +67,10 @@ TIME_GROUPS: Dict[Mode, Dict[Frequency, Tuple[DateTimeComponent, ...]]] = {
 SeasonConfigInput = TypedDict(
     "SeasonConfigInput",
     {
-        "dec_mode": Literal["DJF", "JFD"],
+        # TODO: Deprecate incomplete_djf.
         "drop_incomplete_djf": bool,
+        "drop_incomplete_seasons": bool,
+        "dec_mode": Literal["DJF", "JFD"],
         "custom_seasons": Optional[List[List[str]]],
     },
     total=False,
@@ -76,16 +79,20 @@ SeasonConfigInput = TypedDict(
 SeasonConfigAttr = TypedDict(
     "SeasonConfigAttr",
     {
-        "dec_mode": Literal["DJF", "JFD"],
+        # TODO: Deprecate incomplete_djf.
         "drop_incomplete_djf": bool,
+        "drop_incomplete_seasons": bool,
+        "dec_mode": Literal["DJF", "JFD"],
         "custom_seasons": Optional[Dict[str, List[str]]],
     },
     total=False,
 )
 
 DEFAULT_SEASON_CONFIG: SeasonConfigInput = {
-    "dec_mode": "DJF",
+    # TODO: Deprecate incomplete_djf.
     "drop_incomplete_djf": False,
+    "drop_incomplete_seasons": False,
+    "dec_mode": "DJF",
     "custom_seasons": None,
 }
 
@@ -104,6 +111,7 @@ MONTH_INT_TO_STR: Dict[int, str] = {
     11: "Nov",
     12: "Dec",
 }
+MONTH_STR_TO_INT = {v: k for k, v in MONTH_INT_TO_STR.items()}
 
 # A dictionary mapping pre-defined seasons to their middle month. This
 # dictionary is used during the creation of datetime objects, which don't
@@ -248,6 +256,11 @@ class TemporalAccessor:
         Time bounds are used for generating weights to calculate weighted group
         averages (refer to the ``weighted`` parameter documentation below).
 
+        .. deprecated:: v0.8.0
+            The ``season_config`` dictionary argument ``"drop_incomplete_djf"``
+            is being deprecated. Please use ``"drop_incomplete_seasons"``
+            instead.
+
         Parameters
         ----------
         data_var: str
@@ -280,35 +293,48 @@ class TemporalAccessor:
         keep_weights : bool, optional
             If calculating averages using weights, keep the weights in the
             final dataset output, by default False.
-        season_config: SeasonConfigInput, optional
+        season_config : SeasonConfigInput, optional
             A dictionary for "season" frequency configurations. If configs for
             predefined seasons are passed, configs for custom seasons are
             ignored and vice versa.
 
-            Configs for predefined seasons:
-
-            * "dec_mode" (Literal["DJF", "JFD"], by default "DJF")
-                The mode for the season that includes December.
+            * "drop_incomplete_seasons" (bool, by default False)
+                Seasons are considered incomplete if they do not have all of
+                the required months to form the season. This argument supersedes
+                "drop_incomplete_djf". For example, if we have
+                the time coordinates ["2000-11-16", "2000-12-16", "2001-01-16",
+                "2001-02-16"] and we want to group seasons by "ND" ("Nov",
+                "Dec") and "JFM" ("Jan", "Feb", "Mar").
 
-                * "DJF": season includes the previous year December.
-                * "JFD": season includes the same year December.
-                    Xarray labels the season with December as "DJF", but it is
-                    actually "JFD".
+                * ["2000-11-16", "2000-12-16"] is considered a complete "ND"
+                    season since both "Nov" and "Dec" are present.
+                * ["2001-01-16", "2001-02-16"] is considered an incomplete "JFM"
+                    season because it only has "Jan" and "Feb". Therefore, these
+                    time coordinates are dropped.
 
             * "drop_incomplete_djf" (bool, by default False)
                 If the "dec_mode" is "DJF", this flag drops (True) or keeps
                 (False) time coordinates that fall under incomplete DJF seasons
                 Incomplete DJF seasons include the start year Jan/Feb and the
-                end year Dec.
+                end year Dec. This argument is superceded by
+                "drop_incomplete_seasons" and will be deprecated in a future
+                release.
 
-            Configs for custom seasons:
+            * "dec_mode" (Literal["DJF", "JFD"], by default "DJF")
+                The mode for the season that includes December in the list of
+                list of pre-defined seasons ("DJF"/"JFD", "MAM", "JJA", "SON").
+                This config is ignored if the ``custom_seasons`` config is set.
+
+                * "DJF": season includes the previous year December.
+                * "JFD": season includes the same year December.
+                    Xarray labels the season with December as "DJF", but it is
+                    actually "JFD".
 
             * "custom_seasons" ([List[List[str]]], by default None)
                 List of sublists containing month strings, with each sublist
                 representing a custom season.
 
                 * Month strings must be in the three letter format (e.g., 'Jan')
-                * Each month must be included once in a custom season
                 * Order of the months in each custom season does not matter
                 * Custom seasons can vary in length
 
@@ -349,7 +375,7 @@ class TemporalAccessor:
         >>>     "season",
         >>>     season_config={
         >>>         "dec_mode": "DJF",
-        >>>         "drop_incomplete_season": True
+        >>>         "drop_incomplete_seasons": True
         >>>     }
         >>> )
         >>> ds_season.ts
@@ -385,7 +411,7 @@ class TemporalAccessor:
             'freq': 'season',
             'weighted': 'True',
             'dec_mode': 'DJF',
-            'drop_incomplete_djf': 'False'
+            'drop_incomplete_seasons': 'False'
         }
         """
         self._set_data_var_attrs(data_var)
@@ -414,6 +440,11 @@ class TemporalAccessor:
         Time bounds are used for generating weights to calculate weighted
         climatology (refer to the ``weighted`` parameter documentation below).
 
+        .. deprecated:: v0.8.0
+            The ``season_config`` dictionary argument ``"drop_incomplete_djf"``
+            is being deprecated. Please use ``"drop_incomplete_seasons"``
+            instead.
+
         Parameters
         ----------
         data_var: str
@@ -455,35 +486,48 @@ class TemporalAccessor:
             'yyyy-mm-dd'. For example, ``('1850-01-01', '1899-12-31')``. If no
             value is provided, the climatological reference period will be the
             full period covered by the dataset.
-        season_config: SeasonConfigInput, optional
+        season_config : SeasonConfigInput, optional
             A dictionary for "season" frequency configurations. If configs for
             predefined seasons are passed, configs for custom seasons are
             ignored and vice versa.
 
-            Configs for predefined seasons:
+            * "drop_incomplete_seasons" (bool, by default False)
+                Seasons are considered incomplete if they do not have all of
+                the required months to form the season. This argument supersedes
+                "drop_incomplete_djf". For example, if we have
+                the time coordinates ["2000-11-16", "2000-12-16", "2001-01-16",
+                "2001-02-16"] and we want to group seasons by "ND" ("Nov",
+                "Dec") and "JFM" ("Jan", "Feb", "Mar").
 
-            * "dec_mode" (Literal["DJF", "JFD"], by default "DJF")
-               The mode for the season that includes December.
-
-                * "DJF": season includes the previous year December.
-                * "JFD": season includes the same year December.
-                    Xarray labels the season with December as "DJF", but it is
-                    actually "JFD".
+                * ["2000-11-16", "2000-12-16"] is considered a complete "ND"
+                    season since both "Nov" and "Dec" are present.
+                * ["2001-01-16", "2001-02-16"] is considered an incomplete "JFM"
+                    season because it only has "Jan" and "Feb". Therefore, these
+                    time coordinates are dropped.
 
             * "drop_incomplete_djf" (bool, by default False)
                 If the "dec_mode" is "DJF", this flag drops (True) or keeps
                 (False) time coordinates that fall under incomplete DJF seasons
                 Incomplete DJF seasons include the start year Jan/Feb and the
-                end year Dec.
+                end year Dec. This argument is superceded by
+                "drop_incomplete_seasons" and will be deprecated in a future
+                release.
 
-            Configs for custom seasons:
+            * "dec_mode" (Literal["DJF", "JFD"], by default "DJF")
+                The mode for the season that includes December in the list of
+                list of pre-defined seasons ("DJF"/"JFD", "MAM", "JJA", "SON").
+                This config is ignored if the ``custom_seasons`` config is set.
+
+                * "DJF": season includes the previous year December.
+                * "JFD": season includes the same year December.
+                    Xarray labels the season with December as "DJF", but it is
+                    actually "JFD".
 
             * "custom_seasons" ([List[List[str]]], by default None)
                 List of sublists containing month strings, with each sublist
                 representing a custom season.
 
                 * Month strings must be in the three letter format (e.g., 'Jan')
-                * Each month must be included once in a custom season
                 * Order of the months in each custom season does not matter
                 * Custom seasons can vary in length
 
@@ -528,7 +572,7 @@ class TemporalAccessor:
         >>>     "season",
         >>>     season_config={
         >>>         "dec_mode": "DJF",
-        >>>         "drop_incomplete_season": True
+        >>>         "drop_incomplete_seasons": True
         >>>     }
         >>> )
         >>> ds_season.ts
@@ -564,7 +608,7 @@ class TemporalAccessor:
             'freq': 'season',
             'weighted': 'True',
             'dec_mode': 'DJF',
-            'drop_incomplete_djf': 'False'
+            'drop_incomplete_seasons': 'False'
         }
         """
         self._set_data_var_attrs(data_var)
@@ -600,6 +644,11 @@ class TemporalAccessor:
         Time bounds are used for generating weights to calculate weighted
         climatology (refer to the ``weighted`` parameter documentation below).
 
+        .. deprecated:: v0.8.0
+            The ``season_config`` dictionary argument ``"drop_incomplete_djf"``
+            is being deprecated. Please use ``"drop_incomplete_seasons"``
+            instead.
+
         Parameters
         ----------
         data_var: str
@@ -642,11 +691,35 @@ class TemporalAccessor:
             ``('1850-01-01', '1899-12-31')``. If no value is provided, the
             climatological reference period will be the full period covered by
             the dataset.
-        season_config: SeasonConfigInput, optional
+        season_config : SeasonConfigInput, optional
             A dictionary for "season" frequency configurations. If configs for
             predefined seasons are passed, configs for custom seasons are
             ignored and vice versa.
 
+            General configs:
+
+            * "drop_incomplete_seasons" (bool, by default False)
+                Seasons are considered incomplete if they do not have all of
+                the required months to form the season. This argument supersedes
+                "drop_incomplete_djf". For example, if we have
+                the time coordinates ["2000-11-16", "2000-12-16", "2001-01-16",
+                "2001-02-16"] and we want to group seasons by "ND" ("Nov",
+                "Dec") and "JFM" ("Jan", "Feb", "Mar").
+
+                * ["2000-11-16", "2000-12-16"] is considered a complete "ND"
+                    season since both "Nov" and "Dec" are present.
+                * ["2001-01-16", "2001-02-16"] is considered an incomplete "JFM"
+                    season because it only has "Jan" and "Feb". Therefore, these
+                    time coordinates are dropped.
+
+            * "drop_incomplete_djf" (bool, by default False)
+                If the "dec_mode" is "DJF", this flag drops (True) or keeps
+                (False) time coordinates that fall under incomplete DJF seasons
+                Incomplete DJF seasons include the start year Jan/Feb and the
+                end year Dec. This argument is superceded by
+                "drop_incomplete_seasons" and will be deprecated in a future
+                release.
+
             Configs for predefined seasons:
 
             * "dec_mode" (Literal["DJF", "JFD"], by default "DJF")
@@ -657,12 +730,6 @@ class TemporalAccessor:
                     Xarray labels the season with December as "DJF", but it is
                     actually "JFD".
 
-            * "drop_incomplete_djf" (bool, by default False)
-                If the "dec_mode" is "DJF", this flag drops (True) or keeps
-                (False) time coordinates that fall under incomplete DJF seasons
-                Incomplete DJF seasons include the start year Jan/Feb and the
-                end year Dec.
-
             Configs for custom seasons:
 
             * "custom_seasons" ([List[List[str]]], by default None)
@@ -670,7 +737,6 @@ class TemporalAccessor:
                 representing a custom season.
 
                 * Month strings must be in the three letter format (e.g., 'Jan')
-                * Each month must be included once in a custom season
                 * Order of the months in each custom season does not matter
                 * Custom seasons can vary in length
 
@@ -730,7 +796,7 @@ class TemporalAccessor:
             'frequency': 'season',
             'weighted': 'True',
             'dec_mode': 'DJF',
-            'drop_incomplete_djf': 'False'
+            'drop_incomplete_seasons': 'False'
         }
         """
         # 1. Set the attributes for this instance of `TemporalAccessor`.
@@ -931,30 +997,46 @@ class TemporalAccessor:
             self._is_valid_reference_period(reference_period)
             self._reference_period = reference_period
 
-        # "season" frequency specific configuration attributes.
+        self._set_season_config_attr(season_config)
+
+    def _set_season_config_attr(self, season_config: SeasonConfigInput):
         for key in season_config.keys():
-            if key not in DEFAULT_SEASON_CONFIG.keys():
+            if key not in DEFAULT_SEASON_CONFIG:
                 raise KeyError(
                     f"'{key}' is not a supported season config. Supported "
                     f"configs include: {DEFAULT_SEASON_CONFIG.keys()}."
                 )
-        custom_seasons = season_config.get("custom_seasons", None)
-        dec_mode = season_config.get("dec_mode", "DJF")
-        drop_incomplete_djf = season_config.get("drop_incomplete_djf", False)
 
         self._season_config: SeasonConfigAttr = {}
-        if custom_seasons is None:
+        self._season_config["drop_incomplete_seasons"] = season_config.get(
+            "drop_incomplete_seasons", False
+        )
+
+        custom_seasons = season_config.get("custom_seasons", None)
+        if custom_seasons is not None:
+            self._season_config["custom_seasons"] = self._form_seasons(custom_seasons)
+        else:
+            dec_mode = season_config.get("dec_mode", "DJF")
             if dec_mode not in ("DJF", "JFD"):
                 raise ValueError(
                     "Incorrect 'dec_mode' key value for `season_config`. "
                     "Supported modes include 'DJF' or 'JFD'."
                 )
+
             self._season_config["dec_mode"] = dec_mode
 
+            # TODO: Deprecate incomplete_djf.
+            drop_incomplete_djf = season_config.get("drop_incomplete_djf", False)
             if dec_mode == "DJF":
+                if drop_incomplete_djf is not False:
+                    warnings.warn(
+                        "The `season_config` argument 'drop_incomplete_djf' is being "
+                        "deprecated. Please use 'drop_incomplete_seasons' instead.",
+                        DeprecationWarning,
+                        stacklevel=2,
+                    )
+
                 self._season_config["drop_incomplete_djf"] = drop_incomplete_djf
-        else:
-            self._season_config["custom_seasons"] = self._form_seasons(custom_seasons)
 
     def _is_valid_reference_period(self, reference_period: Tuple[str, str]):
         try:
@@ -997,10 +1079,6 @@ class TemporalAccessor:
         predefined_months = list(MONTH_INT_TO_STR.values())
         input_months = list(chain.from_iterable(custom_seasons))
 
-        if len(input_months) != len(predefined_months):
-            raise ValueError(
-                "Exactly 12 months were not passed in the list of custom seasons."
-            )
         if len(input_months) != len(set(input_months)):
             raise ValueError(
                 "Duplicate month(s) were found in the list of custom seasons."
@@ -1020,9 +1098,13 @@ class TemporalAccessor:
     def _preprocess_dataset(self, ds: xr.Dataset) -> xr.Dataset:
         """Preprocess the dataset based on averaging settings.
 
-        Preprocessing operations include:
-          - Drop incomplete DJF seasons (leading/trailing)
-          - Drop leap days
+        Operations include:
+            1. Drop leap days for daily climatologies.
+            2. Subset the dataset based on the reference period.
+            3. Shift years for custom seasons spanning the calendar year.
+            4. Shift Decembers for "DJF" mode and drop incomplete "DJF" seasons,
+               if specified.
+            5. Drop incomplete seasons if specified.
 
         Parameters
         ----------
@@ -1033,13 +1115,6 @@ class TemporalAccessor:
         -------
         xr.Dataset
         """
-        if (
-            self._freq == "season"
-            and self._season_config.get("dec_mode") == "DJF"
-            and self._season_config.get("drop_incomplete_djf") is True
-        ):
-            ds = self._drop_incomplete_djf(ds)
-
         if (
             self._freq == "day"
             and self._mode in ["climatology", "departures"]
@@ -1052,8 +1127,185 @@ class TemporalAccessor:
                 {self.dim: slice(self._reference_period[0], self._reference_period[1])}
             )
 
+        if (
+            self._freq == "season"
+            and self._season_config.get("custom_seasons") is not None
+        ):
+            # Get a flat list of all of the months included in the custom
+            # seasons to determine if the dataset needs to be subsetted
+            # on just those months. For example, if we define a custom season
+            # "NDJFM", we should subset the dataset for time coordinates
+            # belonging to those months.
+            months = self._season_config["custom_seasons"].values()  # type: ignore
+            months = list(chain.from_iterable(months))
+
+            if len(months) != 12:
+                ds = self._subset_coords_for_custom_seasons(ds, months)
+
+            # The years for time coordinates needs to be shifted by 1 for months
+            # that span the calendar because Xarray groups seasons by months
+            # in the same year, rather than the previous year.
+            ds = self._shift_custom_season_years(ds)
+
+        if self._freq == "season" and self._season_config.get("dec_mode") == "DJF":
+            ds = self._shift_djf_decembers(ds)
+
+            # TODO: Deprecate incomplete_djf.
+            if (
+                self._season_config.get("drop_incomplete_djf") is True
+                and self._season_config.get("drop_incomplete_seasons") is False
+            ):
+                ds = self._drop_incomplete_djf(ds)
+
+        if (
+            self._freq == "season"
+            and self._season_config["drop_incomplete_seasons"] is True
+        ):
+            ds = self._drop_incomplete_seasons(ds)
+
         return ds
 
+    def _subset_coords_for_custom_seasons(
+        self, ds: xr.Dataset, months: List[str]
+    ) -> xr.Dataset:
+        """Subsets time coordinates to the months included in custom seasons.
+
+        Parameters
+        ----------
+        ds : xr.Dataset
+            The dataset.
+        months : List[str]
+            A list of months included in custom seasons.
+            Example: ["Nov", "Dec", "Jan"]
+
+        Returns
+        -------
+        xr.Dataset
+            The dataset with time coordinate subsetted to months used in
+            custom seasons.
+        """
+        month_ints = [MONTH_STR_TO_INT[month] for month in months]
+        ds_new = ds.sel({self.dim: ds[self.dim].dt.month.isin(month_ints)})
+
+        return ds_new
+
+    def _shift_custom_season_years(self, ds: xr.Dataset) -> xr.Dataset:
+        """Shifts the year for custom seasons spanning the calendar year.
+
+        A season spans the calendar year if it includes "Jan" and "Jan" is not
+        the first month. For example, for
+        ``custom_seasons = ["Nov", "Dec", "Jan", "Feb", "Mar"]``:
+          - ["Nov", "Dec"] are from the previous year.
+          - ["Jan", "Feb", "Mar"] are from the current year.
+
+        Therefore, ["Nov", "Dec"] need to be shifted a year forward for correct
+        grouping.
+
+        Parameters
+        ----------
+        ds : xr.Dataset
+            The Dataset with time coordinates.
+
+        Returns
+        -------
+        xr.Dataset
+            The Dataset with shifted time coordinates.
+
+        Examples
+        --------
+
+        Before and after shifting months for "NDJFM" seasons:
+
+        >>> # Before shifting months
+        >>> [(2000, "NDJFM", 11), (2000, "NDJFM", 12), (2001, "NDJFM", 1),
+        >>>  (2001, "NDJFM", 2), (2001, "NDJFM", 3)]
+
+        >>> # After shifting months
+        >>> [(2001, "NDJFM", 11), (2001, "NDJFM", 12), (2001, "NDJFM", 1),
+        >>>  (2001, "NDJFM", 2), (2001, "NDJFM", 3)]
+        """
+        ds_new = ds.copy()
+        custom_seasons = self._season_config["custom_seasons"]
+
+        # Identify months that span across years in custom seasons by getting
+        # the months before "Jan" if "Jan" is not the first month of the season.
+        # Note: Only one custom season can span the calendar year.
+        span_months: List[int] = []
+        for months in custom_seasons.values():  # type: ignore
+            month_ints = [MONTH_STR_TO_INT[month] for month in months]
+
+            if 1 in month_ints and month_ints.index(1) != 0:
+                span_months.extend(month_ints[: month_ints.index(1)])
+            break
+
+        if span_months:
+            time_coords = ds_new[self.dim].copy()
+            indexes = time_coords.dt.month.isin(span_months)
+
+            if isinstance(time_coords.values[0], cftime.datetime):
+                time_coords.values[indexes] = [
+                    time.replace(year=time.year + 1)
+                    for time in time_coords.values[indexes]
+                ]
+            else:
+                time_coords.values[indexes] = [
+                    pd.Timestamp(time) + pd.DateOffset(years=1)
+                    for time in time_coords.values[indexes]
+                ]
+
+            ds_new = ds_new.assign_coords({self.dim: time_coords})
+
+        return ds_new
+
+    def _shift_djf_decembers(self, ds: xr.Dataset) -> xr.Dataset:
+        """Shifts Decembers to the next year for "DJF" seasons.
+
+        This ensures correct grouping for "DJF" seasons by shifting Decembers
+        to the next year. Without this, grouping defaults to "JFD", which
+        is the native Xarray behavior.
+
+        Parameters
+        ----------
+        ds : xr.Dataset
+            The Dataset with time coordinates.
+
+        Returns
+        -------
+        xr.Dataset
+            The Dataset with shifted time coordinates.
+
+        Examples
+        --------
+
+        Comparison of "JFD" and "DJF" seasons:
+
+        >>> # "JFD" (native xarray behavior)
+        >>> [(2000, "DJF", 1), (2000, "DJF", 2), (2000, "DJF", 12),
+        >>>  (2001, "DJF", 1), (2001, "DJF", 2)]
+
+        >>> # "DJF" (shifted Decembers)
+        >>> [(2000, "DJF", 1), (2000, "DJF", 2), (2001, "DJF", 12),
+        >>>  (2001, "DJF", 1), (2001, "DJF", 2)]
+        """
+        ds_new = ds.copy()
+        time_coords = ds_new[self.dim].copy()
+        dec_indexes = time_coords.dt.month == 12
+
+        if isinstance(time_coords.values[0], cftime.datetime):
+            time_coords.values[dec_indexes] = [
+                time.replace(year=time.year + 1)
+                for time in time_coords.values[dec_indexes]
+            ]
+        else:
+            time_coords.values[dec_indexes] = [
+                pd.Timestamp(time) + pd.DateOffset(years=1)
+                for time in time_coords.values[dec_indexes]
+            ]
+
+        ds_new = ds_new.assign_coords({self.dim: time_coords})
+
+        return ds_new
+
     def _drop_incomplete_djf(self, dataset: xr.Dataset) -> xr.Dataset:
         """Drops incomplete DJF seasons within a continuous time series.
 
@@ -1067,7 +1319,6 @@ class TemporalAccessor:
         ----------
         dataset : xr.Dataset
             The dataset with some possibly incomplete DJF seasons.
-
         Returns
         -------
         xr.Dataset
@@ -1102,6 +1353,75 @@ class TemporalAccessor:
 
         return ds_final
 
+    def _drop_incomplete_seasons(self, ds: xr.Dataset) -> xr.Dataset:
+        """Drops incomplete seasons within a continuous time series.
+
+        Seasons are considered incomplete if they do not have all of the
+        required months to form the season. For example, if we have the time
+        coordinates ["2000-11-16", "2000-12-16", "2001-01-16", "2001-02-16"]
+        and we want to group seasons by "ND" ("Nov", "Dec") and "JFM" ("Jan",
+        "Feb", "Mar").
+          - ["2000-11-16", "2000-12-16"] is considered a complete "ND" season
+            since both "Nov" and "Dec" are present.
+          - ["2001-01-16", "2001-02-16"] is considered an incomplete "JFM"
+            season because it only has "Jan" and "Feb". Therefore, these
+            time coordinates are dropped.
+
+        Parameters
+        ----------
+        df : pd.DataFrame
+            A DataFrame of seasonal datetime components with potentially
+            incomplete seasons.
+
+        Returns
+        -------
+        pd.DataFrame
+            A DataFrame of seasonal datetime components with only complete
+            seasons.
+
+        Notes
+        -----
+        TODO: Refactor this method to use pure Xarray/NumPy operations, rather
+        than Pandas.
+        """
+        # Transform the time coords into a DataFrame of seasonal datetime
+        # components based on the grouping mode.
+        time_coords = ds[self.dim].copy()
+        df = self._get_df_dt_components(time_coords, drop_obsolete_cols=False)
+
+        # Get the expected and actual number of months for each season group.
+        df["expected_months"] = df["season"].str.split(r"(?<=.)(?=[A-Z])").str.len()
+        df["actual_months"] = df.groupby(["year", "season"])["year"].transform("count")
+
+        # Get the incomplete seasons and drop the time coordinates that are in
+        # those incomplete seasons.
+        indexes_to_drop = df[df["expected_months"] != df["actual_months"]].index
+
+        if len(indexes_to_drop) == len(time_coords):
+            raise RuntimeError(
+                "No time coordinates remain with `drop_incomplete_seasons=True`. "
+                "Check the dataset has at least one complete season and/or "
+                "specify `drop_incomplete_seasons=False` instead."
+            )
+        elif len(indexes_to_drop) > 0:
+            # The dataset needs to be split into a dataset with and a dataset
+            # without the time dimension because the xarray `.where()` method
+            # adds the time dimension to non-time dimension data vars when
+            # broadcasting, which is a behavior we do not desire.
+            # https://github.com/pydata/xarray/issues/1234
+            # https://github.com/pydata/xarray/issues/8796#issuecomment-1974878267
+            ds_no_time = ds.get([v for v in ds.data_vars if self.dim not in ds[v].dims])  # type: ignore
+            ds_time = ds.get([v for v in ds.data_vars if self.dim in ds[v].dims])  # type: ignore
+
+            coords_to_drop = time_coords.values[indexes_to_drop]
+            ds_time = ds_time.where(~time_coords.isin(coords_to_drop), drop=True)
+
+            ds_new = xr.merge([ds_time, ds_no_time])
+
+            return ds_new
+
+        return ds
+
     def _drop_leap_days(self, ds: xr.Dataset):
         """Drop leap days from time coordinates.
 
@@ -1120,7 +1440,9 @@ class TemporalAccessor:
         -------
         xr.Dataset
         """
-        ds = ds.sel(**{self.dim: ~((ds.time.dt.month == 2) & (ds.time.dt.day == 29))})
+        ds = ds.sel(
+            **{self.dim: ~((ds[self.dim].dt.month == 2) & (ds[self.dim].dt.day == 29))}
+        )
         return ds
 
     def _average(self, ds: xr.Dataset, data_var: str) -> xr.DataArray:
@@ -1294,9 +1616,9 @@ class TemporalAccessor:
         This methods labels time coordinates for grouping by first extracting
         specific xarray datetime components from time coordinates and storing
         them in a pandas DataFrame. After processing (if necessary) is performed
-        on the DataFrame, it is converted to a numpy array of datetime
-        objects. This numpy serves as the data source for the final
-        DataArray of labeled time coordinates.
+        on the DataFrame, it is converted to a numpy array of datetime objects.
+        This numpy array serves as the data source for the final DataArray of
+        labeled time coordinates.
 
         Parameters
         ----------
@@ -1332,7 +1654,9 @@ class TemporalAccessor:
         >>> Coordinates:
         >>> * time     (time) datetime64[ns] 2000-01-01T00:00:00 ... 2000-04-01T00:00:00
         """
-        df_dt_components: pd.DataFrame = self._get_df_dt_components(time_coords)
+        df_dt_components: pd.DataFrame = self._get_df_dt_components(
+            time_coords, drop_obsolete_cols=True
+        )
         dt_objects = self._convert_df_to_dt(df_dt_components)
 
         time_grouped = xr.DataArray(
@@ -1346,7 +1670,9 @@ class TemporalAccessor:
 
         return time_grouped
 
-    def _get_df_dt_components(self, time_coords: xr.DataArray) -> pd.DataFrame:
+    def _get_df_dt_components(
+        self, time_coords: xr.DataArray, drop_obsolete_cols: bool
+    ) -> pd.DataFrame:
         """Returns a DataFrame of xarray datetime components.
 
         This method extracts the applicable xarray datetime components from each
@@ -1367,6 +1693,12 @@ class TemporalAccessor:
         ----------
         time_coords : xr.DataArray
             The time coordinates.
+        drop_obsolete_cols : bool
+            Drop obsolete columns after processing seasonal DataFrame when
+            ``self._freq="season"``. Set to False to keep datetime columns
+            needed for preprocessing the dataset (e.g,. removing incomplete
+            seasons), and set to True to remove obsolete columns when needing
+            to group time coordinates.
 
         Returns
         -------
@@ -1397,41 +1729,18 @@ class TemporalAccessor:
             if self._mode in ["climatology", "departures"]:
                 df["year"] = time_coords[f"{self.dim}.year"].values
                 df["month"] = time_coords[f"{self.dim}.month"].values
-
-            if self._mode == "group_average":
+            elif self._mode == "group_average":
                 df["month"] = time_coords[f"{self.dim}.month"].values
 
-            df = self._process_season_df(df)
-
-        return df
-
-    def _process_season_df(self, df: pd.DataFrame) -> pd.DataFrame:
-        """
-        Processes a DataFrame of datetime components for the season frequency.
-
-        Parameters
-        ----------
-        df : pd.DataFrame
-            A DataFrame of xarray datetime components.
-
-        Returns
-        -------
-        pd.DataFrame
-            A DataFrame of processed xarray datetime components.
-        """
-        df_new = df.copy()
-        custom_seasons = self._season_config.get("custom_seasons")
-        dec_mode = self._season_config.get("dec_mode")
+            custom_seasons = self._season_config.get("custom_seasons")
+            if custom_seasons is not None:
+                df = self._map_months_to_custom_seasons(df)
 
-        if custom_seasons is not None:
-            df_new = self._map_months_to_custom_seasons(df_new)
-        else:
-            if dec_mode == "DJF":
-                df_new = self._shift_decembers(df_new)
+            if drop_obsolete_cols:
+                df = self._drop_obsolete_columns(df)
+                df = self._map_seasons_to_mid_months(df)
 
-        df_new = self._drop_obsolete_columns(df_new)
-        df_new = self._map_seasons_to_mid_months(df_new)
-        return df_new
+        return df
 
     def _map_months_to_custom_seasons(self, df: pd.DataFrame) -> pd.DataFrame:
         """Maps the month column in the DataFrame to a custom season.
@@ -1467,45 +1776,6 @@ class TemporalAccessor:
 
         return df_new
 
-    def _shift_decembers(self, df_season: pd.DataFrame) -> pd.DataFrame:
-        """Shifts Decembers over to the next year for "DJF" seasons in-place.
-
-        For "DJF" seasons, Decembers must be shifted over to the next year in
-        order for the xarray groupby operation to correctly label and group the
-        corresponding time coordinates. If the aren't shifted over, grouping is
-        incorrectly performed with the native xarray "DJF" season (which is
-        actually "JFD").
-
-        Parameters
-        ----------
-        df_season : pd.DataFrame
-            The DataFrame of xarray datetime components produced using the
-            "season" frequency.
-
-        Returns
-        -------
-        pd.DataFrame
-            The DataFrame of xarray datetime components with Decembers shifted
-            over to the next year.
-
-        Examples
-        --------
-
-        Comparison of "JFD" and "DJF" seasons:
-
-        >>> # "JFD" (native xarray behavior)
-        >>> [(2000, "DJF", 1), (2000, "DJF", 2), (2000, "DJF", 12),
-        >>>  (2001, "DJF", 1), (2001, "DJF", 2)]
-
-        >>> # "DJF" (shifted Decembers)
-        >>> [(2000, "DJF", 1), (2000, "DJF", 2), (2001, "DJF", 12),
-        >>>  (2001, "DJF", 1), (2001, "DJF", 2)]
-
-        """
-        df_season.loc[df_season["month"] == 12, "year"] = df_season["year"] + 1
-
-        return df_season
-
     def _map_seasons_to_mid_months(self, df: pd.DataFrame) -> pd.DataFrame:
         """Maps the season column values to the integer of its middle month.
 
@@ -1682,17 +1952,22 @@ class TemporalAccessor:
         )
 
         if self._freq == "season":
-            custom_seasons = self._season_config.get("custom_seasons")
-
-            if custom_seasons is None:
-                dec_mode = self._season_config.get("dec_mode")
-                drop_incomplete_djf = self._season_config.get("drop_incomplete_djf")
+            drop_incomplete_seasons = self._season_config["drop_incomplete_seasons"]
+            drop_incomplete_djf = self._season_config.get("drop_incomplete_djf", False)
 
-                data_var.attrs["dec_mode"] = dec_mode
-                if dec_mode == "DJF":
-                    data_var.attrs["drop_incomplete_djf"] = str(drop_incomplete_djf)
+            # TODO: Deprecate drop_incomplete_djf. This attr is only set if the
+            # user does not set drop_incomplete_seasons.
+            if drop_incomplete_seasons is False and drop_incomplete_djf is not False:
+                data_var.attrs["drop_incomplete_djf"] = str(drop_incomplete_djf)
             else:
+                data_var.attrs["drop_incomplete_seasons"] = str(drop_incomplete_seasons)
+
+            custom_seasons = self._season_config.get("custom_seasons")
+            if custom_seasons is not None:
                 data_var.attrs["custom_seasons"] = list(custom_seasons.keys())
+            else:
+                dec_mode = self._season_config.get("dec_mode")
+                data_var.attrs["dec_mode"] = dec_mode
 
         return data_var
 


Test patch:
diff --git a/tests/test_temporal.py b/tests/test_temporal.py
index e5489b1..7951d3f 100644
--- a/tests/test_temporal.py
+++ b/tests/test_temporal.py
@@ -1,4 +1,5 @@
 import logging
+import warnings
 
 import cftime
 import numpy as np
@@ -209,6 +210,7 @@ class TestAverage:
         )
 
         xr.testing.assert_allclose(result, expected)
+        assert result.ts.attrs == expected.ts.attrs
 
         # Test unweighted averages
         result = ds.temporal.average("ts", weighted=False)
@@ -226,6 +228,7 @@ class TestAverage:
             },
         )
         xr.testing.assert_allclose(result, expected)
+        assert result.ts.attrs == expected.ts.attrs
 
     def test_averages_for_daily_time_series(self):
         ds = xr.Dataset(
@@ -571,27 +574,28 @@ class TestGroupAverage:
         assert result.ts.attrs == expected.ts.attrs
         assert result.time.attrs == expected.time.attrs
 
-    def test_weighted_seasonal_averages_with_DJF_and_drop_incomplete_seasons(self):
+    def test_weighted_seasonal_averages_with_DJF_without_dropping_incomplete_seasons(
+        self,
+    ):
         ds = self.ds.copy()
 
         result = ds.temporal.group_average(
             "ts",
             "season",
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": False},
         )
         expected = ds.copy()
-        # Drop the incomplete DJF seasons
-        expected = expected.isel(time=slice(2, -1))
         expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
             name="ts",
-            data=np.array([[[1]], [[1]], [[1]], [[2.0]]]),
+            data=np.array([[[2.0]], [[1.0]], [[1.0]], [[1.0]], [[2.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
                 "time": xr.DataArray(
                     data=np.array(
                         [
+                            cftime.DatetimeGregorian(2000, 1, 1),
                             cftime.DatetimeGregorian(2000, 4, 1),
                             cftime.DatetimeGregorian(2000, 7, 1),
                             cftime.DatetimeGregorian(2000, 10, 1),
@@ -614,35 +618,82 @@ class TestGroupAverage:
                 "mode": "group_average",
                 "freq": "season",
                 "weighted": "True",
+                "drop_incomplete_seasons": "False",
                 "dec_mode": "DJF",
-                "drop_incomplete_djf": "True",
             },
         )
 
         xr.testing.assert_identical(result, expected)
 
-    def test_weighted_seasonal_averages_with_DJF_without_dropping_incomplete_seasons(
-        self,
-    ):
+    def test_weighted_seasonal_averages_with_DJF_and_drop_incomplete_seasons(self):
+        ds = generate_dataset(decode_times=True, cf_compliant=True, has_bounds=True)
+
+        result = ds.temporal.group_average(
+            "ts",
+            "season",
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
+        )
+
+        expected = ds.copy()
+        expected = expected.drop_dims("time")
+        expected["ts"] = xr.DataArray(
+            name="ts",
+            data=np.ones((4, 4, 4)),
+            coords={
+                "lat": expected.lat,
+                "lon": expected.lon,
+                "time": xr.DataArray(
+                    data=np.array(
+                        [
+                            cftime.DatetimeGregorian(2000, 4, 1),
+                            cftime.DatetimeGregorian(2000, 7, 1),
+                            cftime.DatetimeGregorian(2000, 10, 1),
+                            cftime.DatetimeGregorian(2001, 1, 1),
+                        ],
+                    ),
+                    dims=["time"],
+                    attrs={
+                        "axis": "T",
+                        "long_name": "time",
+                        "standard_name": "time",
+                        "bounds": "time_bnds",
+                    },
+                ),
+            },
+            dims=["time", "lat", "lon"],
+            attrs={
+                "operation": "temporal_avg",
+                "mode": "group_average",
+                "freq": "season",
+                "weighted": "True",
+                "drop_incomplete_seasons": "True",
+                "dec_mode": "DJF",
+            },
+        )
+
+        xr.testing.assert_identical(result, expected)
+
+    def test_weighted_seasonal_averages_with_DJF_and_drop_incomplete_djf(self):
         ds = self.ds.copy()
 
         result = ds.temporal.group_average(
             "ts",
             "season",
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": False},
+            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
         )
         expected = ds.copy()
+        # Drop the incomplete DJF seasons
+        expected = expected.isel(time=slice(2, -1))
         expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
             name="ts",
-            data=np.array([[[2.0]], [[1.0]], [[1.0]], [[1.0]], [[2.0]]]),
+            data=np.array([[[1]], [[1]], [[1]], [[2.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
                 "time": xr.DataArray(
                     data=np.array(
                         [
-                            cftime.DatetimeGregorian(2000, 1, 1),
                             cftime.DatetimeGregorian(2000, 4, 1),
                             cftime.DatetimeGregorian(2000, 7, 1),
                             cftime.DatetimeGregorian(2000, 10, 1),
@@ -666,7 +717,7 @@ class TestGroupAverage:
                 "freq": "season",
                 "weighted": "True",
                 "dec_mode": "DJF",
-                "drop_incomplete_djf": "False",
+                "drop_incomplete_djf": "True",
             },
         )
 
@@ -725,12 +776,57 @@ class TestGroupAverage:
                 "mode": "group_average",
                 "freq": "season",
                 "weighted": "True",
+                "drop_incomplete_seasons": "False",
                 "dec_mode": "JFD",
             },
         )
 
         xr.testing.assert_identical(result, expected)
 
+    def test_raises_error_with_incorrect_custom_seasons_argument(self):
+        # Test raises error with non-3 letter strings
+        with pytest.raises(ValueError):
+            custom_seasons = [
+                ["J", "Feb", "Mar"],
+                ["Apr", "May", "Jun"],
+                ["Jul", "Aug", "Sep"],
+                ["Oct", "Nov", "Dec"],
+            ]
+            self.ds.temporal.group_average(
+                "ts",
+                "season",
+                season_config={"custom_seasons": custom_seasons},
+            )
+
+        # Test raises error if duplicate month(s) were found
+        with pytest.raises(ValueError):
+            custom_seasons = [
+                ["Jan", "Jan", "Mar"],
+                ["Apr", "May", "Jun"],
+                ["Jul", "Aug", "Sep"],
+                ["Oct", "Nov", "Dec"],
+            ]
+            self.ds.temporal.group_average(
+                "ts",
+                "season",
+                season_config={"custom_seasons": custom_seasons},
+            )
+
+    def test_raises_error_with_dataset_that_has_no_complete_seasons(self):
+        ds = self.ds.copy()
+        ds = ds.isel(time=slice(0, 1))
+        custom_seasons = [["Dec", "Jan"]]
+
+        with pytest.raises(RuntimeError):
+            ds.temporal.group_average(
+                "ts",
+                "season",
+                season_config={
+                    "custom_seasons": custom_seasons,
+                    "drop_incomplete_seasons": True,
+                },
+            )
+
     def test_weighted_custom_seasonal_averages(self):
         ds = self.ds.copy()
 
@@ -777,60 +873,194 @@ class TestGroupAverage:
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "season",
+                "weighted": "True",
+                "drop_incomplete_seasons": "False",
                 "custom_seasons": [
                     "JanFebMar",
                     "AprMayJun",
                     "JulAugSep",
                     "OctNovDec",
                 ],
+            },
+        )
+
+        xr.testing.assert_identical(result, expected)
+
+    def test_weighted_seasonal_averages_with_custom_seasons_and_all_complete_seasons(
+        self,
+    ):
+        ds = self.ds.copy()
+        ds["time"].values[:] = np.array(
+            [
+                "2000-01-16T12:00:00.000000000",
+                "2000-02-15T12:00:00.000000000",
+                "2000-03-16T12:00:00.000000000",
+                "2000-06-16T00:00:00.000000000",
+                "2000-09-16T00:00:00.000000000",
+            ],
+            dtype="datetime64[ns]",
+        )
+
+        result = ds.temporal.group_average(
+            "ts",
+            "season",
+            season_config={
+                "custom_seasons": [["Jan", "Mar", "Jun"], ["Feb", "Sep"]],
+                "drop_incomplete_seasons": True,
+            },
+        )
+        expected = ds.copy()
+        expected = expected.drop_dims("time")
+        expected["ts"] = xr.DataArray(
+            name="ts",
+            data=np.array([[[1.34065934]], [[1.47457627]]]),
+            coords={
+                "lat": expected.lat,
+                "lon": expected.lon,
+                "time": xr.DataArray(
+                    data=np.array(
+                        [
+                            cftime.DatetimeGregorian(2000, 3, 1),
+                            cftime.DatetimeGregorian(2000, 9, 1),
+                        ],
+                    ),
+                    dims=["time"],
+                    attrs={
+                        "axis": "T",
+                        "long_name": "time",
+                        "standard_name": "time",
+                        "bounds": "time_bnds",
+                    },
+                ),
+            },
+            dims=["time", "lat", "lon"],
+            attrs={
+                "test_attr": "test",
+                "operation": "temporal_avg",
+                "mode": "group_average",
+                "freq": "season",
                 "weighted": "True",
+                "drop_incomplete_seasons": "True",
+                "custom_seasons": ["JanMarJun", "FebSep"],
+            },
+        )
+
+        xr.testing.assert_allclose(result, expected)
+        assert result.ts.attrs == expected.ts.attrs
+
+    def test_weighted_custom_seasonal_averages_drops_incomplete_seasons(self):
+        ds = self.ds.copy()
+        ds["time"].values[:] = np.array(
+            [
+                "2000-11-16T12:00:00.000000000",
+                "2000-12-16T12:00:00.000000000",
+                "2001-01-16T00:00:00.000000000",
+                "2001-02-16T00:00:00.000000000",
+                "2001-03-16T00:00:00.000000000",
+            ],
+            dtype="datetime64[ns]",
+        )
+
+        custom_seasons = [["Nov", "Dec"], ["Feb", "Mar", "Apr"]]
+
+        result = ds.temporal.group_average(
+            "ts",
+            "season",
+            season_config={
+                "drop_incomplete_seasons": True,
+                "custom_seasons": custom_seasons,
+            },
+        )
+        expected = ds.copy()
+        expected = expected.drop_dims("time")
+        expected["ts"] = xr.DataArray(
+            name="ts",
+            data=np.array([[[1.5]]]),
+            coords={
+                "lat": expected.lat,
+                "lon": expected.lon,
+                "time": xr.DataArray(
+                    data=np.array([cftime.datetime(2000, 12, 1)], dtype=object),
+                    dims=["time"],
+                    attrs={
+                        "axis": "T",
+                        "long_name": "time",
+                        "standard_name": "time",
+                        "bounds": "time_bnds",
+                    },
+                ),
+            },
+            dims=["time", "lat", "lon"],
+            attrs={
+                "test_attr": "test",
+                "operation": "temporal_avg",
+                "mode": "group_average",
+                "freq": "season",
+                "weighted": "True",
+                "drop_incomplete_seasons": "True",
+                "custom_seasons": ["NovDec", "FebMarApr"],
             },
         )
 
         xr.testing.assert_identical(result, expected)
 
-    def test_raises_error_with_incorrect_custom_seasons_argument(self):
-        # Test raises error with non-3 letter strings
-        with pytest.raises(ValueError):
-            custom_seasons = [
-                ["J", "Feb", "Mar"],
-                ["Apr", "May", "Jun"],
-                ["Jul", "Aug", "Sep"],
-                ["Oct", "Nov", "Dec"],
-            ]
-            self.ds.temporal.group_average(
-                "ts",
-                "season",
-                season_config={"custom_seasons": custom_seasons},
-            )
+    def test_weighted_custom_seasonal_averages_with_seasons_spanning_calendar_years(
+        self,
+    ):
+        ds = self.ds.copy()
+        ds["time"].values[:] = np.array(
+            [
+                "2000-11-16T12:00:00.000000000",
+                "2000-12-16T12:00:00.000000000",
+                "2001-01-16T00:00:00.000000000",
+                "2001-02-16T00:00:00.000000000",
+                "2001-03-16T00:00:00.000000000",
+            ],
+            dtype="datetime64[ns]",
+        )
 
-        # Test raises error with missing month(s)
-        with pytest.raises(ValueError):
-            custom_seasons = [
-                ["Feb", "Mar"],
-                ["Apr", "May", "Jun"],
-                ["Jul", "Aug", "Sep"],
-                ["Oct", "Nov", "Dec"],
-            ]
-            self.ds.temporal.group_average(
-                "ts",
-                "season",
-                season_config={"custom_seasons": custom_seasons},
-            )
+        custom_seasons = [
+            ["Nov", "Dec", "Jan", "Feb", "Mar"],
+        ]
 
-        # Test raises error if duplicate month(s) were found
-        with pytest.raises(ValueError):
-            custom_seasons = [
-                ["Jan", "Jan", "Mar"],
-                ["Apr", "May", "Jun"],
-                ["Jul", "Aug", "Sep"],
-                ["Oct", "Nov", "Dec"],
-            ]
-            self.ds.temporal.group_average(
-                "ts",
-                "season",
-                season_config={"custom_seasons": custom_seasons},
-            )
+        result = ds.temporal.group_average(
+            "ts",
+            "season",
+            season_config={"custom_seasons": custom_seasons},
+        )
+        expected = ds.copy()
+        expected = expected.drop_dims("time")
+        expected["ts"] = xr.DataArray(
+            name="ts",
+            data=np.array([[[1.3933333]]]),
+            coords={
+                "lat": expected.lat,
+                "lon": expected.lon,
+                "time": xr.DataArray(
+                    data=np.array([cftime.datetime(2001, 1, 1)], dtype=object),
+                    dims=["time"],
+                    attrs={
+                        "axis": "T",
+                        "long_name": "time",
+                        "standard_name": "time",
+                        "bounds": "time_bnds",
+                    },
+                ),
+            },
+            dims=["time", "lat", "lon"],
+            attrs={
+                "test_attr": "test",
+                "operation": "temporal_avg",
+                "mode": "group_average",
+                "freq": "season",
+                "weighted": "True",
+                "drop_incomplete_seasons": "False",
+                "custom_seasons": ["NovDecJanFebMar"],
+            },
+        )
+
+        xr.testing.assert_allclose(result, expected)
+        assert result.ts.attrs == expected.ts.attrs
 
     def test_weighted_monthly_averages(self):
         ds = self.ds.copy()
@@ -1051,7 +1281,7 @@ class TestClimatology:
             ds.temporal.climatology(
                 "ts",
                 "season",
-                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+                season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
                 reference_period=("01-01-2000", "01-01-2000"),
             )
 
@@ -1059,7 +1289,7 @@ class TestClimatology:
             ds.temporal.climatology(
                 "ts",
                 "season",
-                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+                season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
                 reference_period=("01-01-2000"),
             )
 
@@ -1069,7 +1299,7 @@ class TestClimatology:
         result = ds.temporal.climatology(
             "ts",
             "season",
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
             reference_period=("2000-01-01", "2000-06-01"),
         )
 
@@ -1101,7 +1331,7 @@ class TestClimatology:
                 "freq": "season",
                 "weighted": "True",
                 "dec_mode": "DJF",
-                "drop_incomplete_djf": "True",
+                "drop_incomplete_seasons": "True",
             },
         )
 
@@ -1113,7 +1343,7 @@ class TestClimatology:
         result = ds.temporal.climatology(
             "ts",
             "season",
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
         )
 
         expected = ds.copy()
@@ -1154,8 +1384,72 @@ class TestClimatology:
                 "mode": "climatology",
                 "freq": "season",
                 "weighted": "True",
+                "drop_incomplete_seasons": "True",
                 "dec_mode": "DJF",
+            },
+        )
+
+        xr.testing.assert_identical(result, expected)
+
+    def test_raises_deprecation_warning_with_drop_incomplete_djf_season_config(self):
+        # NOTE: This will test will also cover the other public APIs that
+        # have drop_incomplete_djf as a season_config arg.
+        ds = self.ds.copy()
+
+        with warnings.catch_warnings(record=True) as w:
+            result = ds.temporal.climatology(
+                "ts",
+                "season",
+                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            )
+
+            assert len(w) == 1
+            assert issubclass(w[0].category, DeprecationWarning)
+            assert str(w[0].message) == (
+                "The `season_config` argument 'drop_incomplete_djf' is being deprecated. "
+                "Please use 'drop_incomplete_seasons' instead."
+            )
+
+        expected = ds.copy()
+        expected = expected.drop_dims("time")
+        expected_time = xr.DataArray(
+            data=np.array(
+                [
+                    cftime.DatetimeGregorian(1, 1, 1),
+                    cftime.DatetimeGregorian(1, 4, 1),
+                    cftime.DatetimeGregorian(1, 7, 1),
+                    cftime.DatetimeGregorian(1, 10, 1),
+                ],
+            ),
+            coords={
+                "time": np.array(
+                    [
+                        cftime.DatetimeGregorian(1, 1, 1),
+                        cftime.DatetimeGregorian(1, 4, 1),
+                        cftime.DatetimeGregorian(1, 7, 1),
+                        cftime.DatetimeGregorian(1, 10, 1),
+                    ],
+                ),
+            },
+            attrs={
+                "axis": "T",
+                "long_name": "time",
+                "standard_name": "time",
+                "bounds": "time_bnds",
+            },
+        )
+        expected["ts"] = xr.DataArray(
+            name="ts",
+            data=np.ones((4, 4, 4)),
+            coords={"lat": expected.lat, "lon": expected.lon, "time": expected_time},
+            dims=["time", "lat", "lon"],
+            attrs={
+                "operation": "temporal_avg",
+                "mode": "climatology",
+                "freq": "season",
+                "weighted": "True",
                 "drop_incomplete_djf": "True",
+                "dec_mode": "DJF",
             },
         )
 
@@ -1168,7 +1462,7 @@ class TestClimatology:
         result = ds.temporal.climatology(
             "ts",
             "season",
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
         )
 
         expected = ds.copy()
@@ -1210,7 +1504,7 @@ class TestClimatology:
                 "freq": "season",
                 "weighted": "True",
                 "dec_mode": "DJF",
-                "drop_incomplete_djf": "True",
+                "drop_incomplete_seasons": "True",
             },
         )
 
@@ -1261,6 +1555,7 @@ class TestClimatology:
                 "mode": "climatology",
                 "freq": "season",
                 "weighted": "True",
+                "drop_incomplete_seasons": "False",
                 "dec_mode": "JFD",
             },
         )
@@ -1319,6 +1614,7 @@ class TestClimatology:
                 "mode": "climatology",
                 "freq": "season",
                 "weighted": "True",
+                "drop_incomplete_seasons": "False",
                 "custom_seasons": [
                     "JanFebMar",
                     "AprMayJun",
@@ -1330,6 +1626,57 @@ class TestClimatology:
 
         xr.testing.assert_identical(result, expected)
 
+    def test_weighted_custom_seasonal_climatology_with_seasons_spanning_calendar_years(
+        self,
+    ):
+        ds = self.ds.copy()
+
+        custom_seasons = [["Nov", "Dec", "Jan", "Feb", "Mar"]]
+        result = ds.temporal.climatology(
+            "ts",
+            "season",
+            season_config={
+                "drop_incomplete_seasons": False,
+                "custom_seasons": custom_seasons,
+            },
+        )
+
+        expected = ds.copy()
+        expected = expected.drop_dims("time")
+        expected_time = xr.DataArray(
+            data=np.array(
+                [cftime.DatetimeGregorian(1, 1, 1)],
+            ),
+            coords={
+                "time": np.array(
+                    [cftime.DatetimeGregorian(1, 1, 1)],
+                ),
+            },
+            attrs={
+                "axis": "T",
+                "long_name": "time",
+                "standard_name": "time",
+                "bounds": "time_bnds",
+            },
+        )
+
+        expected["ts"] = xr.DataArray(
+            name="ts",
+            data=np.ones((1, 4, 4)),
+            coords={"lat": expected.lat, "lon": expected.lon, "time": expected_time},
+            dims=["time", "lat", "lon"],
+            attrs={
+                "operation": "temporal_avg",
+                "mode": "climatology",
+                "freq": "season",
+                "weighted": "True",
+                "drop_incomplete_seasons": "False",
+                "custom_seasons": ["NovDecJanFebMar"],
+            },
+        )
+
+        xr.testing.assert_identical(result, expected)
+
     def test_weighted_monthly_climatology(self):
         result = self.ds.temporal.climatology("ts", "month")
 
@@ -1749,7 +2096,7 @@ class TestDepartures:
             ds.temporal.departures(
                 "ts",
                 "season",
-                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+                season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
                 reference_period=("01-01-2000", "01-01-2000"),
             )
 
@@ -1757,7 +2104,7 @@ class TestDepartures:
             ds.temporal.departures(
                 "ts",
                 "season",
-                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+                season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
                 reference_period=("01-01-2000"),
             )
 
@@ -1768,7 +2115,7 @@ class TestDepartures:
             "ts",
             "season",
             weighted=True,
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": False},
             reference_period=("2000-01-01", "2000-06-01"),
         )
 
@@ -1776,13 +2123,14 @@ class TestDepartures:
         expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
             name="ts",
-            data=np.array([[[0.0]], [[np.nan]], [[np.nan]], [[np.nan]]]),
+            data=np.array([[[0.0]], [[0.0]], [[np.nan]], [[np.nan]], [[0.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
                 "time": xr.DataArray(
                     data=np.array(
                         [
+                            cftime.DatetimeGregorian(2000, 1, 1),
                             cftime.DatetimeGregorian(2000, 4, 1),
                             cftime.DatetimeGregorian(2000, 7, 1),
                             cftime.DatetimeGregorian(2000, 10, 1),
@@ -1806,7 +2154,7 @@ class TestDepartures:
                 "freq": "season",
                 "weighted": "True",
                 "dec_mode": "DJF",
-                "drop_incomplete_djf": "True",
+                "drop_incomplete_seasons": "False",
             },
         )
 
@@ -1821,7 +2169,7 @@ class TestDepartures:
             "ts",
             "month",
             weighted=True,
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": True},
             reference_period=("2000-01-01", "2000-06-01"),
         )
 
@@ -1904,20 +2252,21 @@ class TestDepartures:
             "ts",
             "season",
             weighted=True,
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": False},
         )
 
         expected = ds.copy()
         expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
             name="ts",
-            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
+            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
                 "time": xr.DataArray(
                     data=np.array(
                         [
+                            cftime.DatetimeGregorian(2000, 1, 1),
                             cftime.DatetimeGregorian(2000, 4, 1),
                             cftime.DatetimeGregorian(2000, 7, 1),
                             cftime.DatetimeGregorian(2000, 10, 1),
@@ -1941,7 +2290,7 @@ class TestDepartures:
                 "freq": "season",
                 "weighted": "True",
                 "dec_mode": "DJF",
-                "drop_incomplete_djf": "True",
+                "drop_incomplete_seasons": "False",
             },
         )
 
@@ -1955,20 +2304,21 @@ class TestDepartures:
             "season",
             weighted=True,
             keep_weights=True,
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": False},
         )
 
         expected = ds.copy()
         expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
             name="ts",
-            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
+            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
                 "time": xr.DataArray(
                     data=np.array(
                         [
+                            cftime.DatetimeGregorian(2000, 1, 1),
                             cftime.DatetimeGregorian(2000, 4, 1),
                             cftime.DatetimeGregorian(2000, 7, 1),
                             cftime.DatetimeGregorian(2000, 10, 1),
@@ -1992,16 +2342,17 @@ class TestDepartures:
                 "freq": "season",
                 "weighted": "True",
                 "dec_mode": "DJF",
-                "drop_incomplete_djf": "True",
+                "drop_incomplete_seasons": "False",
             },
         )
         expected["time_wts"] = xr.DataArray(
             name="ts",
-            data=np.array([1.0, 1.0, 1.0, 1.0]),
+            data=np.array([0.52542373, 1.0, 1.0, 1.0, 0.47457627]),
             coords={
                 "time_original": xr.DataArray(
                     data=np.array(
                         [
+                            "2000-01-16T12:00:00.000000000",
                             "2000-03-16T12:00:00.000000000",
                             "2000-06-16T00:00:00.000000000",
                             "2000-09-16T00:00:00.000000000",
@@ -2021,7 +2372,8 @@ class TestDepartures:
             dims=["time_original"],
         )
 
-        xr.testing.assert_identical(result, expected)
+        xr.testing.assert_allclose(result, expected)
+        assert result.ts.attrs == expected.ts.attrs
 
     def test_unweighted_seasonal_departures_with_DJF(self):
         ds = self.ds.copy()
@@ -2030,20 +2382,21 @@ class TestDepartures:
             "ts",
             "season",
             weighted=False,
-            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            season_config={"dec_mode": "DJF", "drop_incomplete_seasons": False},
         )
 
         expected = ds.copy()
         expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
             name="ts",
-            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
+            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
                 "time": xr.DataArray(
                     data=np.array(
                         [
+                            cftime.DatetimeGregorian(2000, 1, 1),
                             cftime.DatetimeGregorian(2000, 4, 1),
                             cftime.DatetimeGregorian(2000, 7, 1),
                             cftime.DatetimeGregorian(2000, 10, 1),
@@ -2066,8 +2419,8 @@ class TestDepartures:
                 "mode": "departures",
                 "freq": "season",
                 "weighted": "False",
+                "drop_incomplete_seasons": "False",
                 "dec_mode": "DJF",
-                "drop_incomplete_djf": "True",
             },
         )
 
@@ -2117,6 +2470,7 @@ class TestDepartures:
                 "mode": "departures",
                 "freq": "season",
                 "weighted": "False",
+                "drop_incomplete_seasons": "False",
                 "dec_mode": "JFD",
             },
         )
@@ -3290,7 +3644,7 @@ class Test_Averager:
                 weighted=True,
                 season_config={
                     "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
+                    "drop_incomplete_seasons": False,
                     "custom_seasons": None,
                 },
             )
@@ -3306,7 +3660,7 @@ class Test_Averager:
                 weighted=True,
                 season_config={
                     "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
+                    "drop_incomplete_seasons": False,
                     "custom_seasons": None,
                 },
             )
@@ -3318,7 +3672,7 @@ class Test_Averager:
                 weighted=True,
                 season_config={
                     "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
+                    "drop_incomplete_seasons": False,
                     "custom_seasons": None,
                 },
             )
@@ -3330,7 +3684,7 @@ class Test_Averager:
                 weighted=True,
                 season_config={
                     "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
+                    "drop_incomplete_seasons": False,
                     "custom_seasons": None,
                 },
             )
@@ -3356,7 +3710,7 @@ class Test_Averager:
                 weighted=True,
                 season_config={
                     "dec_mode": "unsupported",
-                    "drop_incomplete_djf": False,
+                    "drop_incomplete_seasons": False,
                     "custom_seasons": None,
                 },
             )