Repository: xCDAT/xcdat

Issue:
[Feature]: Allow user selection of reference period in departures
### Is your feature request related to a problem?

The current departures functionality determines the anomalies relative to the average value over the full time series. 

Climate scientists often look at the departures relative to a specific period (rather than the full period). For example, see [IPCC AR6 WG1 Fig. 1.12](https://www.ipcc.ch/report/ar6/wg1/figures/chapter-1/figure-1-12/) which shows the global mean surface temperature departures relative to the average temperature over 1850 – 1900. This shows zero **average** anomaly between 1850 and 1900 and emphasizes the rapid warming that follows. 

We should make it possible for xCDAT users to do this kind of calculation easily. 

### Describe the solution you'd like

xCDAT should allow for flexible selection of the climatological reference period with which to define departures. I think there are several non-mutually exclusive ways to accomplish this:

* Allow the user to provide a a custom climatology when calculating departures, e.g., `ds.temporal.departures(data_var, freq, climatology=ds_climo)`
* Allow the user to specify a climatological base period when calculating departures, e.g., `ds.temporal.departures(data_var, freq, reference_period=("1850-01-01", "1899-12-31"))`
* Or both (possibly by augmenting climatology functionality to accept a `reference_period`, i.e., `ds.temporal.climatology(data_var, freq, reference_period=("1850-01-01", "1899-12-31"))`

### Describe alternatives you've considered

An alternative manual method would entail:

* Subset the dataset of interest to the climatological reference period (`ds_climo = ds.sel("1850-01-01", "1899-12-31")`)
* Take the climatology of this subsetted data (`ds_climo = ds.temporal.group_average(data_var, freq="month")`)
* Somehow subtract this climatology from the original timeseries (not sure how to do this or if this would work): `ds_departure[data_var] = ds[data_var] - ds_climo[data_var]`

### Additional context

- Related discussion: https://github.com/xCDAT/xcdat/discussions/404
- Jiwoo's xarray+xcdat implementation: https://github.com/xCDAT/xcdat_test/blob/main/cookbook/seasonal_departure_from_baseline_clim.ipynb

Discussion:
codecov-commenter: # [Codecov](https://codecov.io/gh/xCDAT/xcdat/pull/417?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) Report
Base: **100.00**% // Head: **100.00**% // No change to project coverage :thumbsup:
> Coverage data is based on head [(`7e56fb2`)](https://codecov.io/gh/xCDAT/xcdat/pull/417?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) compared to base [(`1efb807`)](https://codecov.io/gh/xCDAT/xcdat/commit/1efb8074e2cbec99c2b410b001c979926d8a44df?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).
> Patch coverage: 100.00% of modified lines in pull request are covered.

:mega: This organization is not using Codecov’s [GitHub App Integration](https://github.com/apps/codecov). We recommend you install it so Codecov can continue to function properly for your repositories. [Learn more](https://about.codecov.io/blog/codecov-is-updating-its-github-integration/?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      #417   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           14        14           
  Lines         1266      1278   +12     
=========================================
+ Hits          1266      1278   +12     
```


| [Impacted Files](https://codecov.io/gh/xCDAT/xcdat/pull/417?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) | Coverage Δ | |
|---|---|---|
| [xcdat/temporal.py](https://codecov.io/gh/xCDAT/xcdat/pull/417?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT#diff-eGNkYXQvdGVtcG9yYWwucHk=) | `100.00% <100.00%> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT)

</details>

[:umbrella: View full report at Codecov](https://codecov.io/gh/xCDAT/xcdat/pull/417?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).

tomvothecoder: >  * One thing I noticed is that calling `ds.temporal.departures` drops `time_bnds` – which I don't think needs to happen (I assume xarray doesn't know how to deal with this during the subtraction of the climatology). I think the time bounds can be added back in after subtraction. I noticed this in verifying that this PR was accurate (I get accuracy of order 10-4 – 10-7, I think because I added time bounds back in using `add_missing_bounds`, which is problematic).

Before this PR:
- `departures()` used to maintain `time_bnds` because the output frequency was the same as the input frequency (no `group_average()` call).

In this PR:
- Since we now call `group_average()`, `time_bnds` are dropped because the `time` coordinates are grouped which changes the `time` frequency and shape.

Solutions/Workarounds:
1. As you mentioned, get new `time_bnds` for the new grouped `time` coordinates using `add_missing_bounds()`. The accuracy _should_ improve with #412.
2. Keep a separate `time_original` dimension/coordinates to keep `time_bnds` in the dataset (otherwise we get `time_bnds` with `NaT` values since `time` frequency and shape changes).

---

>  * Steps zero and one of `departures` call `_set_arg_attrs()`, `.copy()`, `_preprocess_dataset()`, and `group_average()`. This overlaps strongly with `_averager()` – maybe things could be streamlined here (though I see that some of the order-of-operations is different, which may be why they are different.

I saw your point and played around with different implements using `_average()` to tried to find ways to streamline the implementation. I wasn't able to find a simple way to call `self._average()` without making significant modifications to how the `TemporalAccessor` class attributes are set (e.g., `self._mode`). **We can probably try refactoring this in the future if we wanted to.**

I updated the order-of-operations and made clearer comments in the `departures()` method so it is hopefully easier to understand. Please check out this [commit](https://github.com/xCDAT/xcdat/pull/417/commits/25832ddf5357bf043ecbd300aa03e71308840987) diff and let me know if those comments make sense.

---

>   * As noted, departures calls `group_average()` on each call – I think this may not be necesary if the input/output frequency are the same?

I also don't think we need to perform group averaging if the input/output frequencies are the same. **It still works either way** (and converts time coordinates to `cftime` if they aren't already), but it might improve performance by removing unnecessary operations.

We can try determining if `group_average()` is necessary by using the existing `_infer_freq()` method and comparing the result to the input `freq` argument. 



Patch:
diff --git a/docs/examples/climatology-and-departures.ipynb b/docs/examples/climatology-and-departures.ipynb
index 57e8420..001a936 100644
--- a/docs/examples/climatology-and-departures.ipynb
+++ b/docs/examples/climatology-and-departures.ipynb
@@ -12,7 +12,7 @@
     "\n",
     "Date: 05/27/22\n",
     "\n",
-    "Last Updated: 09/26/22 (v0.3.3)\n",
+    "Last Updated: 2/27/23\n",
     "\n",
     "Related APIs:\n",
     "\n",
@@ -1341,6 +1341,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -1353,7 +1354,11 @@
     "Helpful knowledge:\n",
     "\n",
     "* Masked (missing) data is automatically handled.\n",
-    "  * The weight of masked (missing) data is excluded when averages are calculated. This is the same as giving them a weight of 0."
+    "  * The weight of masked (missing) data is excluded when averages are calculated. This is the same as giving them a weight of 0.\n",
+    "* If desired, use the `reference_period` argument to calculate a climatology based on a\n",
+    "  climatological reference period (a subset of the entire time series). If no value is\n",
+    "  provided, the climatological reference period will be the full period covered by the\n",
+    "  dataset."
    ]
   },
   {
@@ -4238,6 +4243,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -4256,6 +4262,10 @@
     "* How is the climatology calculated?\n",
     "  * In the departures API, the reference climatology is calculated internally so there is no need  to pass one to this method.\n",
     "  * You can still calculate the reference climatology using the climatology API.\n",
+    "  * If desired, use the `reference_period` argument to calculate anomalies relative to a\n",
+    "  climatological reference period (a subset of the entire time series). If no value is\n",
+    "  provided, the climatological reference period will be the full period covered by the\n",
+    "  dataset.\n",
     "* Masked (missing) data is automatically handled.\n",
     "  * The weight of masked (missing) data is excluded when averages are calculated. This is the same as giving them a weight of 0.\n"
    ]
@@ -6608,7 +6618,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.10.8"
   },
   "toc": {
    "base_numbering": 1,
diff --git a/xcdat/temporal.py b/xcdat/temporal.py
index f4f0cc7..3120d48 100644
--- a/xcdat/temporal.py
+++ b/xcdat/temporal.py
@@ -1,4 +1,5 @@
 """Module containing temporal functions."""
+from datetime import datetime
 from itertools import chain
 from typing import Dict, List, Literal, Optional, Tuple, TypedDict, get_args
 
@@ -152,11 +153,13 @@ class TemporalAccessor:
         Returns a Dataset with the average of a data variable and the time
         dimension removed.
 
-        This method is particularly useful for calculating the weighted averages
-        of monthly or yearly time series data because the number of days per
-        month/year can vary based on the calendar type, which can affect
-        weighting. For other frequencies, the distribution of weights will be
-        equal so ``weighted=True`` is the same as ``weighted=False``.
+        This method infers the time grouping frequency by checking the distance
+        between a set of upper and lower time bounds. This method is
+        particularly useful for calculating the weighted averages of monthly or
+        yearly time series data because the number of days per month/year can
+        vary based on the calendar type, which can affect weighting. For other
+        frequencies, the distribution of weights will be equal so
+        ``weighted=True`` is the same as ``weighted=False``.
 
         Parameters
         ----------
@@ -192,11 +195,14 @@ class TemporalAccessor:
         >>> ds_month = ds.temporal.average("ts")
         >>> ds_month.ts
         """
+        # Set the data variable related attributes (e.g., dim name, calendar)
         self._set_data_var_attrs(data_var)
 
         freq = self._infer_freq()
 
-        return self._averager(data_var, "average", freq, weighted, keep_weights)
+        return self._averager(
+            data_var, "average", freq, weighted=weighted, keep_weights=keep_weights
+        )
 
     def group_average(
         self,
@@ -332,7 +338,12 @@ class TemporalAccessor:
         self._set_data_var_attrs(data_var)
 
         return self._averager(
-            data_var, "group_average", freq, weighted, keep_weights, season_config
+            data_var,
+            "group_average",
+            freq,
+            weighted=weighted,
+            keep_weights=keep_weights,
+            season_config=season_config,
         )
 
     def climatology(
@@ -341,6 +352,7 @@ class TemporalAccessor:
         freq: Frequency,
         weighted: bool = True,
         keep_weights: bool = False,
+        reference_period: Optional[Tuple[str, str]] = None,
         season_config: SeasonConfigInput = DEFAULT_SEASON_CONFIG,
     ):
         """Returns a Dataset with the climatology of a data variable.
@@ -375,6 +387,12 @@ class TemporalAccessor:
         keep_weights : bool, optional
             If calculating averages using weights, keep the weights in the
             final dataset output, by default False.
+        reference_period : Optional[Tuple[str, str]], optional
+            The climatological reference period, which is a subset of the entire
+            time series. This parameter accepts a tuple of strings in the format
+            '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
             A dictionary for "season" frequency configurations. If configs for
             predefined seasons are passed, configs for custom seasons are
@@ -476,7 +494,13 @@ class TemporalAccessor:
         self._set_data_var_attrs(data_var)
 
         return self._averager(
-            data_var, "climatology", freq, weighted, keep_weights, season_config
+            data_var,
+            "climatology",
+            freq,
+            weighted,
+            keep_weights,
+            reference_period,
+            season_config,
         )
 
     def departures(
@@ -485,6 +509,7 @@ class TemporalAccessor:
         freq: Frequency,
         weighted: bool = True,
         keep_weights: bool = False,
+        reference_period: Optional[Tuple[str, str]] = None,
         season_config: SeasonConfigInput = DEFAULT_SEASON_CONFIG,
     ) -> xr.Dataset:
         """
@@ -496,17 +521,6 @@ class TemporalAccessor:
         temperature) and the long-term average value for that time interval
         (e.g., the average surface temperature over the last 30 Januaries).
 
-        This method uses xarray's grouped arithmetic as a shortcut for mapping
-        over all unique labels. Grouped arithmetic works by assigning a grouping
-        label to each time coordinate of the observation data based on the
-        averaging mode and frequency. Afterwards, the corresponding climatology
-        is removed from the observation data at each time coordinate based on
-        the matching labels.
-
-        xarray's grouped arithmetic operates over each value of the DataArray
-        corresponding to each grouping label without changing the size of the
-        DataArray. For example,the original monthly time coordinates are
-        maintained when calculating seasonal departures on monthly data.
 
         Parameters
         ----------
@@ -539,6 +553,13 @@ class TemporalAccessor:
         keep_weights : bool, optional
             If calculating averages using weights, keep the weights in the
             final dataset output, by default False.
+        reference_period : Optional[Tuple[str, str]], optional
+            The climatological reference period, which is a subset of the entire
+            time series and used for calculating departures. This parameter
+            accepts a tuple of strings in the format '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
             A dictionary for "season" frequency configurations. If configs for
             predefined seasons are passed, configs for custom seasons are
@@ -586,6 +607,13 @@ class TemporalAccessor:
 
         Notes
         -----
+        This method uses xarray's grouped arithmetic as a shortcut for mapping
+        over all unique labels. Grouped arithmetic works by assigning a grouping
+        label to each time coordinate of the observation data based on the
+        averaging mode and frequency. Afterwards, the corresponding climatology
+        is removed from the observation data at each time coordinate based on
+        the matching labels.
+
         Refer to [3]_ to learn more about how xarray's grouped arithmetic works.
 
         References
@@ -611,48 +639,88 @@ class TemporalAccessor:
             'drop_incomplete_djf': 'False'
         }
         """
-        ds = self._dataset.copy()
+        # 1. Set the attributes for this instance of `TemporalAccessor`.
+        # ----------------------------------------------------------------------
+        self._set_arg_attrs(
+            "departures", freq, weighted, reference_period, season_config
+        )
         self._set_data_var_attrs(data_var)
-        self._set_arg_attrs("departures", freq, weighted, season_config)
 
-        # Preprocess the dataset based on method argument values.
+        # 2. Copy the original dataset and preprocess (if needed) for reuse.
+        # ----------------------------------------------------------------------
+        ds = self._dataset.copy()
+
+        # Preprocess ahead of time to avoid needing to preprocess again when
+        # calling the group average or climatology APIs in step #3.
         ds = self._preprocess_dataset(ds)
 
-        # Group the observation data variable.
-        dv_obs = _get_data_var(ds, data_var)
-        self._labeled_time = self._label_time_coords(dv_obs[self.dim])
-        dv_obs_grouped = self._group_data(dv_obs)
+        # 3. Calculate the grouped average and climatology of the data variable.
+        # ----------------------------------------------------------------------
+        # The climatology and grouped average APIs are called on the copied
+        # dataset to create separate instances of the `TemporalAccessor` class.
+        # This is done to avoid overriding the attributes of the current
+        # instance of `TemporalAccessor` (set in step #1 above).
+
+        # Group averaging is only required if the dataset's frequency (input)
+        # differs from the `freq` arg (output).
+        ds_obs = ds.copy()
+        inferred_freq = self._infer_freq()
+        if inferred_freq != freq:
+            ds_obs = ds_obs.temporal.group_average(
+                data_var,
+                freq,
+                weighted,
+                keep_weights,
+                season_config,
+            )
 
-        # Calculate the climatology of the data variable.
         ds_climo = ds.temporal.climatology(
-            data_var, freq, weighted, keep_weights, season_config
+            data_var,
+            freq,
+            weighted,
+            keep_weights,
+            reference_period,
+            season_config,
         )
-        dv_climo = ds_climo[data_var]
 
-        # Rename the time dimension using the name from the grouped observation
-        # data. The dimension names must align for xarray's grouped arithmetic
-        # to work. Otherwise, the error below is thrown: `ValueError:
-        # incompatible dimensions for a grouped binary operation: the group
-        # variable '<CHOSEN_FREQ>' is not a dimension on the other argument`
+        # 4. Group the averaged data variable values by the time `freq`.
+        # ----------------------------------------------------------------------
+        # This step allows us to perform xarray's grouped arithmetic to
+        # calculate departures.
+        dv_obs = ds_obs[data_var].copy()
+        self._labeled_time = self._label_time_coords(dv_obs[self.dim])
+        dv_obs_grouped = self._group_data(dv_obs)
+
+        # 5. Align time dimension names using the labeled time dimension name.
+        # ----------------------------------------------------------------------
+        # The climatology's time dimension is renamed to the labeled time
+        # dimension in step #4 above (e.g., "time" -> "season"). xarray requires
+        # dimension names to be aligned to perform grouped arithmetic, which we
+        # use for calculating departures in step #5. Otherwise, this error is
+        # raised: "`ValueError: incompatible dimensions for a grouped binary
+        # operation: the group variable '<FREQ ARG>' is not a dimension on the
+        # other argument`".
+        dv_climo = ds_climo[data_var]
         dv_climo = dv_climo.rename({self.dim: self._labeled_time.name})
 
-        # Calculate the departures for the data variable, which uses the formula
-        # observation - climatology.
+        # 6. Calculate the departures for the data variable.
+        # ----------------------------------------------------------------------
+        # departures = observation - climatology
         with xr.set_options(keep_attrs=True):
-            ds_departs = ds.copy()
-            ds_departs[data_var] = dv_obs_grouped - dv_climo
-            ds_departs[data_var] = self._add_operation_attrs(ds_departs[data_var])
+            dv_departs = dv_obs_grouped - dv_climo
+            dv_departs = self._add_operation_attrs(dv_departs)
+            ds_obs[data_var] = dv_departs
 
-            # The original time coordinates are restored after grouped
-            # subtraction. The grouped time coordinates are dropped since they
-            # are no longer required.
-            ds_departs = ds_departs.drop_vars(self._labeled_time.name)
+            # The original time dimension name is restored after grouped
+            # arithmetic, so the labeled time dimension name is no longer needed
+            # and therefore dropped.
+            ds_obs = ds_obs.drop_vars(self._labeled_time.name)
 
         if weighted and keep_weights:
             self._weights = ds_climo.time_wts
-            ds_departs = self._keep_weights(ds_departs)
+            ds_obs = self._keep_weights(ds_obs)
 
-        return ds_departs
+        return ds_obs
 
     def _infer_freq(self) -> Frequency:
         """Infers the time frequency from the coordinates.
@@ -689,11 +757,12 @@ class TemporalAccessor:
         freq: Frequency,
         weighted: bool = True,
         keep_weights: bool = False,
+        reference_period: Optional[Tuple[str, str]] = None,
         season_config: SeasonConfigInput = DEFAULT_SEASON_CONFIG,
     ) -> xr.Dataset:
         """Averages a data variable based on the averaging mode and frequency."""
         ds = self._dataset.copy()
-        self._set_arg_attrs(mode, freq, weighted, season_config)
+        self._set_arg_attrs(mode, freq, weighted, reference_period, season_config)
 
         # Preprocess the dataset based on method argument values.
         ds = self._preprocess_dataset(ds)
@@ -775,6 +844,7 @@ class TemporalAccessor:
         mode: Mode,
         freq: Frequency,
         weighted: bool,
+        reference_period: Optional[Tuple[str, str]] = None,
         season_config: SeasonConfigInput = DEFAULT_SEASON_CONFIG,
     ):
         """Validates method arguments and sets them as object attributes.
@@ -819,6 +889,11 @@ class TemporalAccessor:
         self._freq = freq
         self._weighted = weighted
 
+        self._reference_period = None
+        if reference_period is not None:
+            self._is_valid_reference_period(reference_period)
+            self._reference_period = reference_period
+
         # "season" frequency specific configuration attributes.
         for key in season_config.keys():
             if key not in DEFAULT_SEASON_CONFIG.keys():
@@ -844,6 +919,17 @@ class TemporalAccessor:
         else:
             self._season_config["custom_seasons"] = self._form_seasons(custom_seasons)
 
+    def _is_valid_reference_period(self, reference_period: Tuple[str, str]):
+        try:
+            datetime.strptime(reference_period[0], "%Y-%m-%d")
+            datetime.strptime(reference_period[1], "%Y-%m-%d")
+        except (IndexError, ValueError):
+            raise ValueError(
+                "Reference periods must be a tuple of strings with the format "
+                "'yyyy-mm-dd'. For example, reference_period=('1850-01-01', "
+                "'1899-12-31')."
+            )
+
     def _form_seasons(self, custom_seasons: List[List[str]]) -> Dict[str, List[str]]:
         """Forms custom seasons from a nested list of months.
 
@@ -924,6 +1010,11 @@ class TemporalAccessor:
         ):
             ds = self._drop_leap_days(ds)
 
+        if self._mode == "climatology" and self._reference_period is not None:
+            ds = ds.sel(
+                {self.dim: slice(self._reference_period[0], self._reference_period[1])}
+            )
+
         return ds
 
     def _drop_incomplete_djf(self, dataset: xr.Dataset) -> xr.Dataset:
@@ -1533,11 +1624,6 @@ class TemporalAccessor:
             # Only keep the original time coordinates, not the ones labeled
             # by group.
             self._weights = self._weights.drop_vars(self._labeled_time.name)
-        # Strip "_original" from the name of the weights` time coordinates
-        # because the final departures Dataset has the original time coordinates
-        # restored after performing grouped subtraction.
-        elif self._mode == "departures":
-            self._weights = self._weights.rename({f"{self.dim}_original": self.dim})
 
         ds[self._weights.name] = self._weights
 


Test patch:
diff --git a/tests/test_temporal.py b/tests/test_temporal.py
index cf48116..22be7fe 100644
--- a/tests/test_temporal.py
+++ b/tests/test_temporal.py
@@ -1053,6 +1053,69 @@ class TestClimatology:
             "ds['time'].encoding['calendar'] = 'noleap') and try again."
         ) in caplog.text
 
+    def test_raises_error_if_reference_period_arg_is_incorrect(self):
+        ds = self.ds.copy()
+
+        with pytest.raises(ValueError):
+            ds.temporal.climatology(
+                "ts",
+                "season",
+                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+                reference_period=("01-01-2000", "01-01-2000"),
+            )
+
+        with pytest.raises(ValueError):
+            ds.temporal.climatology(
+                "ts",
+                "season",
+                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+                reference_period=("01-01-2000"),
+            )
+
+    def test_subsets_climatology_based_on_reference_period(self):
+        ds = self.ds.copy()
+
+        result = ds.temporal.climatology(
+            "ts",
+            "season",
+            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            reference_period=("2000-01-01", "2000-06-01"),
+        )
+
+        # The first month Jan/Feb are dropped (incomplete DJF). This means
+        # only the MAM season will be present, with April being the middle month
+        # (represented by month number 4).
+        expected = ds.copy()
+        expected = expected.drop_dims("time")
+        expected_time = xr.DataArray(
+            data=np.array([cftime.DatetimeGregorian(1, 4, 1)]),
+            coords={
+                "time": np.array([cftime.DatetimeGregorian(1, 4, 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",
+                "dec_mode": "DJF",
+                "drop_incomplete_djf": "True",
+            },
+        )
+
+        assert result.identical(expected)
+
     def test_weighted_seasonal_climatology_with_DJF(self):
         ds = self.ds.copy()
 
@@ -1618,15 +1681,54 @@ class TestClimatology:
 
 
 class TestDepartures:
-    # TODO: Update TestDepartures tests to use other numbers rather than 1's for
-    # better test reliability and accuracy. This may require subsetting.
     @pytest.fixture(autouse=True)
     def setup(self):
-        self.ds: xr.Dataset = generate_dataset(
-            decode_times=True, cf_compliant=False, has_bounds=True
+        time = 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",
+                    "2001-02-15T12:00:00.000000000",
+                ],
+                dtype="datetime64[ns]",
+            ),
+            dims=["time"],
+            attrs={"axis": "T", "long_name": "time", "standard_name": "time"},
+        )
+        time.encoding = {"calendar": "standard"}
+        time_bnds = xr.DataArray(
+            name="time_bnds",
+            data=np.array(
+                [
+                    ["2000-01-01T00:00:00.000000000", "2000-02-01T00:00:00.000000000"],
+                    ["2000-03-01T00:00:00.000000000", "2000-04-01T00:00:00.000000000"],
+                    ["2000-06-01T00:00:00.000000000", "2000-07-01T00:00:00.000000000"],
+                    ["2000-09-01T00:00:00.000000000", "2000-10-01T00:00:00.000000000"],
+                    ["2001-02-01T00:00:00.000000000", "2001-03-01T00:00:00.000000000"],
+                ],
+                dtype="datetime64[ns]",
+            ),
+            coords={"time": time},
+            dims=["time", "bnds"],
+            attrs={"xcdat_bounds": "True"},
         )
 
-        self.seasons = ["JJA", "MAM", "SON", "DJF"]
+        self.ds = xr.Dataset(
+            data_vars={"time_bnds": time_bnds},
+            coords={"lat": [-90], "lon": [0], "time": time},
+        )
+        self.ds.time.attrs["bounds"] = "time_bnds"
+
+        self.ds["ts"] = xr.DataArray(
+            data=np.array(
+                [[[2.0]], [[1.0]], [[1.0]], [[1.0]], [[2.0]]], dtype="float64"
+            ),
+            coords={"time": self.ds.time, "lat": self.ds.lat, "lon": self.ds.lon},
+            dims=["time", "lat", "lon"],
+            attrs={"test_attr": "test"},
+        )
 
     def test_raises_error_if_time_coords_are_not_decoded(self):
         ds: xr.Dataset = generate_dataset(
@@ -1654,28 +1756,65 @@ class TestDepartures:
             "ds['time'].encoding['calendar'] = 'noleap') and try again."
         ) in caplog.text
 
-    def test_weighted_seasonal_departures_with_DJF(self):
+    def test_raises_error_if_reference_period_arg_is_incorrect(self):
         ds = self.ds.copy()
 
-        # Drop incomplete DJF seasons
-        ds = ds.isel(time=slice(2, -1))
+        with pytest.raises(ValueError):
+            ds.temporal.departures(
+                "ts",
+                "season",
+                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+                reference_period=("01-01-2000", "01-01-2000"),
+            )
+
+        with pytest.raises(ValueError):
+            ds.temporal.departures(
+                "ts",
+                "season",
+                season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+                reference_period=("01-01-2000"),
+            )
+
+    def test_seasonal_departures_relative_to_climatology_reference_period(self):
+        ds = self.ds.copy()
 
-        # Compare result of the method against the expected.
         result = ds.temporal.departures(
             "ts",
             "season",
+            weighted=True,
             season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            reference_period=("2000-01-01", "2000-06-01"),
         )
+
         expected = ds.copy()
+        expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
-            data=np.zeros((12, 4, 4)),
+            name="ts",
+            data=np.array([[[0.0]], [[np.nan]], [[np.nan]], [[np.nan]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
-                "time": ds.time,
+                "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={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "departures",
                 "freq": "season",
@@ -1687,13 +1826,144 @@ class TestDepartures:
 
         assert result.identical(expected)
 
-    def test_weighted_seasonal_departures_with_DJF_and_keep_weights(self):
+    def test_monthly_departures_relative_to_climatology_reference_period_with_same_output_freq(
+        self,
+    ):
         ds = self.ds.copy()
 
-        # Drop incomplete DJF seasons
-        ds = ds.isel(time=slice(2, -1))
+        result = ds.temporal.departures(
+            "ts",
+            "month",
+            weighted=True,
+            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+            reference_period=("2000-01-01", "2000-06-01"),
+        )
+
+        expected = ds.copy()
+        expected = expected.drop_dims("time")
+        expected["ts"] = xr.DataArray(
+            name="ts",
+            data=np.array([[[0.0]], [[0.0]], [[np.nan]], [[np.nan]], [[np.nan]]]),
+            coords={
+                "lat": expected.lat,
+                "lon": expected.lon,
+                "time": 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",
+                            "2001-02-15T12:00:00.000000000",
+                        ],
+                        dtype="datetime64[ns]",
+                    ),
+                    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": "departures",
+                "freq": "month",
+                "weighted": "True",
+            },
+        )
+        expected["time_bnds"] = xr.DataArray(
+            name="time_bnds",
+            data=np.array(
+                [
+                    [
+                        "2000-01-01T00:00:00.000000000",
+                        "2000-02-01T00:00:00.000000000",
+                    ],
+                    [
+                        "2000-03-01T00:00:00.000000000",
+                        "2000-04-01T00:00:00.000000000",
+                    ],
+                    [
+                        "2000-06-01T00:00:00.000000000",
+                        "2000-07-01T00:00:00.000000000",
+                    ],
+                    [
+                        "2000-09-01T00:00:00.000000000",
+                        "2000-10-01T00:00:00.000000000",
+                    ],
+                    [
+                        "2001-02-01T00:00:00.000000000",
+                        "2001-03-01T00:00:00.000000000",
+                    ],
+                ],
+                dtype="datetime64[ns]",
+            ),
+            dims=["time", "bnds"],
+            attrs={
+                "xcdat_bounds": "True",
+            },
+        )
+
+        assert result.identical(expected)
+
+    def test_weighted_seasonal_departures_with_DJF(self):
+        ds = self.ds.copy()
+
+        result = ds.temporal.departures(
+            "ts",
+            "season",
+            weighted=True,
+            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+        )
+
+        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]]]),
+            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={
+                "test_attr": "test",
+                "operation": "temporal_avg",
+                "mode": "departures",
+                "freq": "season",
+                "weighted": "True",
+                "dec_mode": "DJF",
+                "drop_incomplete_djf": "True",
+            },
+        )
+
+        assert result.identical(expected)
+
+    def test_weighted_seasonal_departures_with_DJF_and_keep_weights(self):
+        ds = self.ds.copy()
 
-        # Compare result of the method against the expected.
         result = ds.temporal.departures(
             "ts",
             "season",
@@ -1701,16 +1971,36 @@ class TestDepartures:
             keep_weights=True,
             season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
         )
+
         expected = ds.copy()
+        expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
-            data=np.zeros((12, 4, 4)),
+            name="ts",
+            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
-                "time": ds.time,
+                "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={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "departures",
                 "freq": "season",
@@ -1720,51 +2010,72 @@ class TestDepartures:
             },
         )
         expected["time_wts"] = xr.DataArray(
-            data=np.array(
-                [
-                    0.33695652,
-                    0.32608696,
-                    0.33695652,
-                    0.32608696,
-                    0.33695652,
-                    0.33695652,
-                    0.32967033,
-                    0.34065934,
-                    0.32967033,
-                    0.34444444,
-                    0.34444444,
-                    0.31111111,
-                ]
-            ),
-            coords={"time": ds.time},
-            dims=["time"],
+            name="ts",
+            data=np.array([1.0, 1.0, 1.0, 1.0]),
+            coords={
+                "time_original": xr.DataArray(
+                    data=np.array(
+                        [
+                            "2000-03-16T12:00:00.000000000",
+                            "2000-06-16T00:00:00.000000000",
+                            "2000-09-16T00:00:00.000000000",
+                            "2001-02-15T12:00:00.000000000",
+                        ],
+                        dtype="datetime64[ns]",
+                    ),
+                    dims=["time_original"],
+                    attrs={
+                        "axis": "T",
+                        "long_name": "time",
+                        "standard_name": "time",
+                        "bounds": "time_bnds",
+                    },
+                )
+            },
+            dims=["time_original"],
         )
 
-        xr.testing.assert_allclose(result, expected)
-        assert result.ts.attrs == expected.ts.attrs
+        assert result.identical(expected)
 
     def test_unweighted_seasonal_departures_with_DJF(self):
         ds = self.ds.copy()
-        # Drop incomplete DJF seasons
-        ds = ds.isel(time=slice(2, -1))
 
-        # Compare result of the method against the expected.
         result = ds.temporal.departures(
             "ts",
             "season",
             weighted=False,
             season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
         )
+
         expected = ds.copy()
+        expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
-            data=np.zeros((12, 4, 4)),
+            name="ts",
+            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
-                "time": ds.time,
+                "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={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "departures",
                 "freq": "season",
@@ -1779,23 +2090,43 @@ class TestDepartures:
     def test_unweighted_seasonal_departures_with_JFD(self):
         ds = self.ds.copy()
 
-        # Compare result of the method against the expected.
         result = ds.temporal.departures(
             "ts",
             "season",
             weighted=False,
             season_config={"dec_mode": "JFD"},
         )
+
         expected = ds.copy()
+        expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
-            data=np.zeros((15, 4, 4)),
+            name="ts",
+            data=np.array([[[0.0]], [[0.0]], [[0.0]], [[0.0]], [[0.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
-                "time": ds.time,
+                "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),
+                            cftime.DatetimeGregorian(2001, 1, 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": "departures",
                 "freq": "season",
@@ -1866,12 +2197,11 @@ class TestDepartures:
                     "time": xr.DataArray(
                         data=np.array(
                             [
-                                "2000-01-16T12:00:00.000000000",
-                                "2000-06-16T00:00:00.000000000",
-                                "2000-09-16T00:00:00.000000000",
-                                "2001-03-15T12:00:00.000000000",
+                                cftime.DatetimeGregorian(2000, 1, 16),
+                                cftime.DatetimeGregorian(2000, 6, 16),
+                                cftime.DatetimeGregorian(2000, 9, 16),
+                                cftime.DatetimeGregorian(2001, 3, 15),
                             ],
-                            dtype="datetime64[ns]",
                         ),
                         dims=["time"],
                         attrs={
@@ -1895,32 +2225,6 @@ class TestDepartures:
                             "weighted": "True",
                         },
                     ),
-                    "time_bnds": xr.DataArray(
-                        name="time_bnds",
-                        data=np.array(
-                            [
-                                [
-                                    "2000-01-01T00:00:00.000000000",
-                                    "2000-02-01T00:00:00.000000000",
-                                ],
-                                [
-                                    "2000-06-01T00:00:00.000000000",
-                                    "2000-07-01T00:00:00.000000000",
-                                ],
-                                [
-                                    "2000-09-01T00:00:00.000000000",
-                                    "2000-10-01T00:00:00.000000000",
-                                ],
-                                [
-                                    "2001-03-01T00:00:00.000000000",
-                                    "2001-04-01T00:00:00.000000000",
-                                ],
-                            ],
-                            dtype="datetime64[ns]",
-                        ),
-                        dims=["time", "bnds"],
-                        attrs={"xcdat_bounds": "True"},
-                    ),
                 },
             )