Repository: xCDAT/xcdat

Issue:
[Feature]: Update `_get_weights()` method for `SpatialAccessor` and `TemporalAccessor`
### Is your feature request related to a problem?

These methods are useful to have publicly so the user can reuse weights for other operations, or to see what the generated values are for validation purposes.

Related to this comment:
>I think this issue can probably be closed since the weights can be retrieved via .get_weights() here. If a more complex implementation of this is needed in the future, that could probably be added as a feature request.
&mdash; @pochedls https://github.com/xCDAT/xcdat/issues/114#issuecomment-1146730677

### Describe the solution you'd like

- Rename `_get_weights()` to `get_weights()`
- Update references to the new name

### Describe alternatives you've considered

_No response_

### Additional context

_No response_

Discussion:
tomvothecoder: > I have no objection to making `get_weights()` public, though I didn't fully understand the the changes to the unit tests. Should I view the changes in `tests/` as a reorganization or something that I should review more carefully?

I addressed this comment in one of the PR review comments.

> 
> I think the error for `test_weights_for_region_in_lon_domain_with_both_spanning_p_meridian` is correct and we should represent the problematic line as `ds.lon_bnds.data[:] = np.array([[-1, 1], [1, 90], [90, 180], [180, 359]])`.

I ended up removing tests where the domain bounds has a larger left hand value than right hand value because **they are supposed to throw an error, rather than pass**.  There's already a test to check that an error is thrown for this case.

I kept `test_weights_for_region_in_lon_domain_with_region_spanning_p_meridian` since **regions can have a larger lower bound vs. upper bound**.
codecov-commenter: # [Codecov](https://codecov.io/gh/xCDAT/xcdat/pull/252?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) Report
> Merging [#252](https://codecov.io/gh/xCDAT/xcdat/pull/252?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) (d4a49f1) into [main](https://codecov.io/gh/xCDAT/xcdat/commit/1fdc8a97d685ae38345a053e4c9ed32ef9cf5f14?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) (1fdc8a9) will **not change** coverage.
> The diff coverage is `100.00%`.

```diff
@@            Coverage Diff            @@
##              main      #252   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files            9         8    -1     
  Lines          742       727   -15     
=========================================
- Hits           742       727   -15     
```


| [Impacted Files](https://codecov.io/gh/xCDAT/xcdat/pull/252?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) | Coverage Δ | |
|---|---|---|
| [xcdat/axis.py](https://codecov.io/gh/xCDAT/xcdat/pull/252/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT#diff-eGNkYXQvYXhpcy5weQ==) | `100.00% <ø> (ø)` | |
| [xcdat/utils.py](https://codecov.io/gh/xCDAT/xcdat/pull/252/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT#diff-eGNkYXQvdXRpbHMucHk=) | `100.00% <ø> (ø)` | |
| [xcdat/\_\_init\_\_.py](https://codecov.io/gh/xCDAT/xcdat/pull/252/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT#diff-eGNkYXQvX19pbml0X18ucHk=) | `100.00% <100.00%> (ø)` | |
| [xcdat/bounds.py](https://codecov.io/gh/xCDAT/xcdat/pull/252/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT#diff-eGNkYXQvYm91bmRzLnB5) | `100.00% <100.00%> (ø)` | |
| [xcdat/dataset.py](https://codecov.io/gh/xCDAT/xcdat/pull/252/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT#diff-eGNkYXQvZGF0YXNldC5weQ==) | `100.00% <100.00%> (ø)` | |
| [xcdat/spatial.py](https://codecov.io/gh/xCDAT/xcdat/pull/252/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT#diff-eGNkYXQvc3BhdGlhbC5weQ==) | `100.00% <100.00%> (ø)` | |
| [xcdat/temporal.py](https://codecov.io/gh/xCDAT/xcdat/pull/252/diff?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%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/xCDAT/xcdat/pull/252?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/xCDAT/xcdat/pull/252?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT). Last update [223e73b...d4a49f1](https://codecov.io/gh/xCDAT/xcdat/pull/252?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT).

tomvothecoder: Hi @pochedls, @lee1043, and @chengzhuzhang, instead of making `_get_weights()` public for temporal averaging so the user can get weights, I implemented a `keep_weights` keyword argument to the temporal averaging APIs. This gives users the choice with whether or not they want to couple `time_weights` with its related output.

| Option 	| Make `_get_weights()` public 	| Option to keep `time_weights` in the output `xr.Dataset` 	|
|---	|---	|---	|
| Pros 	| * User can calculate weights independently of temporal averaging, which makes it faster to get weights instead of running temporal averaging first (do scientists normally calculate weights independently? does CDAT allow have a method for this?) 	| * Makes sense to couple the weights to its related averaging output <br>* Implementation is straightforward, just add the `time_weights` DataArray to the final output Dataset                         <br>* User has the option to not include the weights in the output 	|
| Cons 	| * User needs to provide the same arguments as temporal averaging APIs (`data_var`, `mode`, `freq`, `season_config`) which makes the method less intuitive to implement and use <br>* Weights are decoupled from the output, which removes context on how they should be used for averaging<br>* We don't accept weights in our temporal averaging APIs right now (does it make sense to?) 	| * For group_average and climatology we must append `"_original"` to the name of `time_weights` time dimension/coordinates since datasets share dimensions and coordinates. Otherwise, it will use the reduced  time<br>dimension/coordinates (not correct, will result in  `np.nan` weights) 	|


Example for `ds.temporal.climatology()`:
https://github.com/xCDAT/xcdat/blob/44620c0953350fffa4df69c1ae2050088d31f058/xcdat/temporal.py#L332-L339
https://github.com/xCDAT/xcdat/blob/44620c0953350fffa4df69c1ae2050088d31f058/xcdat/temporal.py#L742-L745

Example for `ds.temporal.departures()`
https://github.com/xCDAT/xcdat/blob/44620c0953350fffa4df69c1ae2050088d31f058/xcdat/temporal.py#L624-L627
https://github.com/xCDAT/xcdat/blob/44620c0953350fffa4df69c1ae2050088d31f058/xcdat/temporal.py#L465-L472

`_keep_weights()` method:
https://github.com/xCDAT/xcdat/blob/44620c0953350fffa4df69c1ae2050088d31f058/xcdat/temporal.py#L1369-L1402


**If this looks good, we can implement `keep_weights` with the spatial averaging API as well.**

chengzhuzhang: Hi @tomvothecoder Thank you for the discussion. I don't recall that I have a use case that requires to get/keep the weights to pass down to another operation. With the exception that, they are useful to be examed during xcdat development/validation. So I don't have strong opinion on the implementation, maybe keeping it private is fine. I will defer to Steve and Jiwoo. 
pochedls: I like `keep_weights`. I didn't totally understand the downside for keep_weights: why can't we just overwrite the original time_weights? Are both the full time series and the climatology time series saved in the dataset when calculating a climatology?

I think it would be good to keep `.get_weights` in spatial averaging, though I don't object to adding a `keep_weights` in that routine. One thing that might confuse people is that the weights we generate for spatial averaging are 2D (we never explicitly generate a time resolved weighting matrix).
tomvothecoder: > I like `keep_weights`. I didn't totally understand the downside for keep_weights: why can't we just overwrite the original time_weights? 

There is only one `time_wts` generated using the `time_bounds`, which is saved in the Dataset when `keep_weights=True`. I'm not sure what you mean by overwriting the original `time_wts`. 

#### Keeping weights with Climatology

* `time_wts` uses the original time coordinates, so it CANNOT share the reduced time dimension and coordinates from the climatology data variable.
  *  `"time_original"` is the original time dimension (1919 coord pts)
  * `"time"` is the reduced time dimension (4 coord pts)

```
ds.temporal.climatology("tas", freq="season", weighted=True, keep_weights=True)

<xarray.Dataset>
Dimensions:        (lon: 320, nv: 2, lat: 160, bnds: 2, time: 4,
                    time_original: 1919)
Coordinates:
  * lon            (lon) float64 0.0 1.125 2.25 3.375 ... 356.6 357.8 358.9
  * lat            (lat) float64 -89.14 -88.03 -86.91 ... 86.91 88.03 89.14
  * time           (time) object 0001-01-01 00:00:00 ... 0001-10-01 00:00:00
  * time_original  (time_original) datetime64[ns] 1850-01-16T12:00:00 ... 200...
Dimensions without coordinates: nv, bnds
Data variables:
    lon_bounds     (lon, nv) float64 dask.array<chunksize=(320, 2), meta=np.ndarray>
    lat_bounds     (lat, nv) float64 dask.array<chunksize=(160, 2), meta=np.ndarray>
    lon_bnds       (lon, bnds) float64 -0.5625 0.5625 0.5625 ... 358.3 359.4
    lat_bnds       (lat, bnds) float64 -89.7 -88.59 -88.59 ... 88.59 88.59 89.7
    tas            (time, lat, lon) float64 dask.array<chunksize=(1, 160, 320), meta=np.ndarray>
    time_wts       (time_original) float64 0.002038 0.002038 ... 0.002083

```

#### Keeping weights with Departures

* `time_wts` uses the original "time" dimension and coordinates, which are maintained from the departures data variable so they can be shared.
* Notice there is no `"time_original"` dim or coords

```
ds.temporal.departures("tas", freq="season", weighted=True, keep_weights=True)

<xarray.Dataset>
Dimensions:     (lon: 320, nv: 2, lat: 160, time: 1919, bnds: 2)
Coordinates:
  * lon         (lon) float64 0.0 1.125 2.25 3.375 ... 355.5 356.6 357.8 358.9
  * lat         (lat) float64 -89.14 -88.03 -86.91 -85.79 ... 86.91 88.03 89.14
  * time        (time) datetime64[ns] 1850-01-16T12:00:00 ... 2009-11-16
Dimensions without coordinates: nv, bnds
Data variables:
    lon_bounds  (lon, nv) float64 dask.array<chunksize=(320, 2), meta=np.ndarray>
    lat_bounds  (lat, nv) float64 dask.array<chunksize=(160, 2), meta=np.ndarray>
    time_bnds   (time, bnds) datetime64[ns] 1850-01-01T18:00:00 ... 2009-12-0...
    tas         (time, lat, lon) float64 dask.array<chunksize=(2, 160, 320), meta=np.ndarray>
    lon_bnds    (lon, bnds) float64 -0.5625 0.5625 0.5625 ... 358.3 358.3 359.4
    lat_bnds    (lat, bnds) float64 -89.7 -88.59 -88.59 ... 88.59 88.59 89.7
    time_wts    (time) float64 0.002038 0.002038 0.002059 ... 0.002083 0.002083

```


> Are both the full time series and the climatology time series saved in the dataset when calculating a climatology?
Only the climatology time series is stored in the climatology output.

> I think it would be good to keep `.get_weights` in spatial averaging, though I don't object to adding a `keep_weights` in that routine. 

Sounds good, I added `keep_weights` as a keyword argument to `ds.spatial.average()`.

> One thing that might confuse people is that the weights we generate for spatial averaging are 2D (we never explicitly generate a time resolved weighting matrix).

Hmmm, we can address this if users open up an issue about it, unless you think it is worth addressing now.



pochedls: Got it – thanks Tom.

Patch:
diff --git a/xcdat/axis.py b/xcdat/axis.py
index 5d1a256..2d57d57 100644
--- a/xcdat/axis.py
+++ b/xcdat/axis.py
@@ -350,6 +350,6 @@ def _get_coord_var(dataset: xr.Dataset, axis: GenericAxis):
         raise KeyError(
             f"A '{axis}' coordinate variable was found in the dataset. Make sure "
             "the coordinate variable exists in the Dataset, and its 'axis' attribute "
-            "is set to '{axis}'."
+            f"is set to '{axis}'."
         )
     return coord_var
diff --git a/xcdat/spatial.py b/xcdat/spatial.py
index bcf2002..fefdd7b 100644
--- a/xcdat/spatial.py
+++ b/xcdat/spatial.py
@@ -64,6 +64,7 @@ class SpatialAccessor:
         data_var: str,
         axis: List[SpatialAxis] = ["X", "Y"],
         weights: Union[Literal["generate"], xr.DataArray] = "generate",
+        keep_weights: bool = False,
         lat_bounds: Optional[RegionAxisBounds] = None,
         lon_bounds: Optional[RegionAxisBounds] = None,
     ) -> xr.Dataset:
@@ -101,6 +102,9 @@ class SpatialAccessor:
             averaging. ``weights`` must include the same spatial axis dimensions
             and have the same dimensional sizes as the data variable, by default
             "generate".
+        keep_weights : bool, optional
+            If calculating averages using weights, keep the weights in the
+            final dataset output, by default False.
         lat_bounds : Optional[RegionAxisBounds], optional
             A tuple of floats/ints for the regional latitude lower and upper
             boundaries. This arg is used when calculating axis weights, but is
@@ -169,8 +173,8 @@ class SpatialAccessor:
         >>> ts_global = ds.spatial.average("tas", axis=["X", "Y"],
         >>>     weights=weights)["tas"]
         """
-        dataset = self._dataset.copy()
-        dv = _get_data_var(dataset, data_var)
+        ds = self._dataset.copy()
+        dv = _get_data_var(ds, data_var)
         self._validate_axis_arg(axis)
 
         if isinstance(weights, str) and weights == "generate":
@@ -178,13 +182,100 @@ class SpatialAccessor:
                 self._validate_region_bounds("Y", lat_bounds)
             if lon_bounds is not None:
                 self._validate_region_bounds("X", lon_bounds)
-            dv_weights = self._get_weights(axis, lat_bounds, lon_bounds)
+            self._weights = self.get_weights(axis, lat_bounds, lon_bounds)
         elif isinstance(weights, xr.DataArray):
-            dv_weights = weights
+            self._weights = weights
 
-        self._validate_weights(dv, axis, dv_weights)
-        dataset[dv.name] = self._averager(dv, axis, dv_weights)
-        return dataset
+        self._validate_weights(dv, axis)
+        ds[dv.name] = self._averager(dv, axis)
+
+        if keep_weights:
+            ds[self._weights.name] = self._weights
+
+        return ds
+
+    def get_weights(
+        self,
+        axis: List[SpatialAxis],
+        lat_bounds: Optional[RegionAxisBounds] = None,
+        lon_bounds: Optional[RegionAxisBounds] = None,
+    ) -> xr.DataArray:
+        """
+        Get area weights for specified axis keys and an optional target domain.
+
+        This method first determines the weights for an individual axis based on
+        the difference between the upper and lower bound. For latitude the
+        weight is determined by the difference of sine(latitude). All axis
+        weights are then combined to form a DataArray of weights that can be
+        used to perform a weighted (spatial) average.
+
+        If ``lat_bounds`` or ``lon_bounds`` are supplied, then grid cells
+        outside this selected regional domain are given zero weight. Grid cells
+        that are partially in this domain are given partial weight.
+
+        Parameters
+        ----------
+        axis : List[SpatialAxis]
+            List of axis dimensions to average over.
+        lat_bounds : Optional[RegionAxisBounds]
+            Tuple of latitude boundaries for regional selection, by default
+            None.
+        lon_bounds : Optional[RegionAxisBounds]
+            Tuple of longitude boundaries for regional selection, by default
+            None.
+
+        Returns
+        -------
+        xr.DataArray
+           A DataArray containing the region weights to use during averaging.
+           ``weights`` are 1-D and correspond to the specified axes (``axis``)
+           in the region.
+
+        Notes
+        -----
+        This method was developed for rectilinear grids only. ``get_weights()``
+        recognizes and operate on latitude and longitude, but could be extended
+        to work with other standard geophysical dimensions (e.g., time, depth,
+        and pressure).
+        """
+        Bounds = TypedDict(
+            "Bounds",
+            {"domain": xr.DataArray, "region": Optional[RegionAxisBounds]},
+        )
+        axis_bounds: Dict[SpatialAxis, Bounds] = {
+            "X": {
+                "domain": self._dataset.bounds.get_bounds("lon").copy(),
+                "region": lon_bounds,
+            },
+            "Y": {
+                "domain": self._dataset.bounds.get_bounds("lat").copy(),
+                "region": lat_bounds,
+            },
+        }
+
+        axis_weights: AxisWeights = {}
+
+        for key in axis:
+            d_bounds = axis_bounds[key]["domain"]
+            self._validate_domain_bounds(d_bounds)
+
+            r_bounds: Union[Optional[RegionAxisBounds], np.ndarray] = axis_bounds[key][
+                "region"
+            ]
+            if r_bounds is not None:
+                r_bounds = np.array(r_bounds, dtype="float")
+
+            if key == "X":
+                weights = self._get_longitude_weights(d_bounds, r_bounds)
+            elif key == "Y":
+                weights = self._get_latitude_weights(d_bounds, r_bounds)
+
+            weights.attrs = d_bounds.attrs
+            axis_weights[key] = weights
+
+        weights = self._combine_weights(axis_weights)
+
+        return weights
 
     def _validate_axis_arg(self, axis: List[SpatialAxis]):
         """
@@ -294,88 +385,6 @@ class SpatialAccessor:
                 "Pass a tuple with the format (lower, upper)."
             )
 
-    def _get_weights(
-        self,
-        axis: List[SpatialAxis],
-        lat_bounds: Optional[RegionAxisBounds],
-        lon_bounds: Optional[RegionAxisBounds],
-    ) -> xr.DataArray:
-        """
-        Get area weights for specified axis keys and an optional target domain.
-
-        This method first determines the weights for an individual axis based on
-        the difference between the upper and lower bound. For latitude the
-        weight is determined by the difference of sine(latitude). All axis
-        weights are then combined to form a DataArray of weights that can be
-        used to perform a weighted (spatial) average.
-
-        If ``lat_bounds`` or ``lon_bounds`` are supplied, then grid cells
-        outside this selected regional domain are given zero weight. Grid cells
-        that are partially in this domain are given partial weight.
-
-        Parameters
-        ----------
-        axis : List[SpatialAxis]
-            List of axis dimensions to average over.
-        lat_bounds : Optional[RegionAxisBounds]
-            Tuple of latitude boundaries for regional selection.
-        lon_bounds : Optional[RegionAxisBounds]
-            Tuple of longitude boundaries for regional selection.
-
-        Returns
-        -------
-        xr.DataArray
-           A DataArray containing the region weights to use during averaging.
-           ``weights`` are 1-D and correspond to the specified axes (``axis``)
-           in the region.
-
-        Notes
-        -----
-        This method was developed for rectilinear grids only. ``_get_weights()``
-        recognizes and operate on latitude and longitude, but could be extended
-        to work with other standard geophysical dimensions (e.g., time, depth,
-        and pressure).
-        """
-        Bounds = TypedDict(
-            "Bounds",
-            {"domain": xr.DataArray, "region": Optional[RegionAxisBounds]},
-        )
-        axis_bounds: Dict[SpatialAxis, Bounds] = {
-            "X": {
-                "domain": self._dataset.bounds.get_bounds("lon").copy(),
-                "region": lon_bounds,
-            },
-            "Y": {
-                "domain": self._dataset.bounds.get_bounds("lat").copy(),
-                "region": lat_bounds,
-            },
-        }
-
-        axis_weights: AxisWeights = {}
-
-        for key in axis:
-            d_bounds = axis_bounds[key]["domain"]
-            self._validate_domain_bounds(d_bounds)
-
-            r_bounds: Union[Optional[RegionAxisBounds], np.ndarray] = axis_bounds[key][
-                "region"
-            ]
-            if r_bounds is not None:
-                r_bounds = np.array(r_bounds, dtype="float")
-
-            if key == "X":
-                weights = self._get_longitude_weights(d_bounds, r_bounds)
-            elif key == "Y":
-                weights = self._get_latitude_weights(d_bounds, r_bounds)
-
-            weights.attrs = d_bounds.attrs
-            weights.name = key + "_wts"
-            axis_weights[key] = weights
-
-        weights = self._combine_weights(axis_weights)
-
-        return weights
-
     def _get_longitude_weights(
         self, domain_bounds: xr.DataArray, region_bounds: Optional[np.ndarray]
     ) -> xr.DataArray:
@@ -639,11 +648,11 @@ class SpatialAccessor:
             (``axis``) in the region.
         """
         region_weights = reduce((lambda x, y: x * y), axis_weights.values())
+        region_weights.name = "_".join(sorted(region_weights.coords.keys())) + "_wts"
+
         return region_weights
 
-    def _validate_weights(
-        self, data_var: xr.DataArray, axis: List[SpatialAxis], weights: xr.DataArray
-    ):
+    def _validate_weights(self, data_var: xr.DataArray, axis: List[SpatialAxis]):
         """Validates the ``weights`` arg based on a set of criteria.
 
         This methods checks for the dimensional alignment between the
@@ -676,13 +685,13 @@ class SpatialAccessor:
         x_key = data_var.cf.axes["X"][0]
         y_key = data_var.cf.axes["Y"][0]
 
-        if "X" in axis and x_key not in weights.dims:
+        if "X" in axis and x_key not in self._weights.dims:
             raise KeyError(
                 "The weights DataArray either does not include an X axis, "
                 "or the X axis coordinates does not have the 'axis' attribute "
                 "set to 'X'."
             )
-        if "Y" in axis and y_key not in weights.dims:
+        if "Y" in axis and y_key not in self._weights.dims:
             raise KeyError(
                 "The weights DataArray either does not include an Y axis, "
                 "or the Y axis coordinates does not have the 'axis' attribute "
@@ -690,17 +699,15 @@ class SpatialAccessor:
             )
 
         # Check the weight dim sizes equal data var dim sizes.
-        dim_sizes = {key: data_var.sizes[key] for key in weights.sizes.keys()}
-        for dim, size in weights.sizes.items():
+        dim_sizes = {key: data_var.sizes[key] for key in self._weights.sizes.keys()}
+        for dim, size in self._weights.sizes.items():
             if size != dim_sizes[dim]:
                 raise ValueError(
-                    f"The axis dimension sizes between supplied `weights` {dict(weights.sizes)} "
+                    f"The axis dimension sizes between supplied `weights` {dict(self._weights.sizes)} "
                     f"and the data variable {dim_sizes} are misaligned."
                 )
 
-    def _averager(
-        self, data_var: xr.DataArray, axis: List[SpatialAxis], weights: xr.DataArray
-    ):
+    def _averager(self, data_var: xr.DataArray, axis: List[SpatialAxis]):
         """Perform a weighted average of a data variable.
 
         This method assumes all specified keys in ``axis`` exists in the data
@@ -718,10 +725,6 @@ class SpatialAccessor:
             Data variable inside a Dataset.
         axis : List[SpatialAxis]
             List of axis dimensions to average over.
-        weights : xr.DataArray
-            A DataArray containing the region area weights for averaging.
-            ``weights`` must include the same spatial axis dimensions and have
-            the same sizes as the data variable.
 
         Returns
         -------
@@ -733,7 +736,7 @@ class SpatialAccessor:
         ``weights`` must be a DataArray and cannot contain missing values.
         Missing values are replaced with 0 using ``weights.fillna(0)``.
         """
-        weights = weights.fillna(0)
+        weights = self._weights.fillna(0)
 
         with xr.set_options(keep_attrs=True):
             weighted_mean = data_var.cf.weighted(weights).mean(axis)
diff --git a/xcdat/temporal.py b/xcdat/temporal.py
index 6d04c78..fb5016e 100644
--- a/xcdat/temporal.py
+++ b/xcdat/temporal.py
@@ -144,12 +144,12 @@ class TemporalAccessor:
 
     def __init__(self, dataset: xr.Dataset):
         self._dataset: xr.Dataset = dataset
+
         self._dim_name = _get_coord_var(self._dataset, "T").name
 
-        # The weights for time coordinates, which are based on a chosen frequency.
-        self._weights: Optional[xr.DataArray] = None
+        self._time_bounds = self._dataset.bounds.get_bounds("time").copy()
 
-    def average(self, data_var: str, weighted: bool = True):
+    def average(self, data_var: str, weighted: bool = True, keep_weights: bool = False):
         """
         Returns a Dataset with the average of a data variable and the time
         dimension removed.
@@ -176,6 +176,9 @@ class TemporalAccessor:
 
             The weight of masked (missing) data is excluded when averages are
             taken. This is the same as giving them a weight of 0.
+        keep_weights : bool, optional
+            If calculating averages using weights, keep the weights in the
+            final dataset output, by default False.
 
         Returns
         -------
@@ -193,41 +196,14 @@ class TemporalAccessor:
         """
         freq = self._infer_freq()
 
-        return self._averager(data_var, "average", freq, weighted)
-
-    def _infer_freq(self) -> Frequency:
-        """Infers the time frequency from the coordinates.
-
-        This method infers the time frequency from the coordinates by
-        calculating the minimum delta and comparing it against a set of
-        conditionals.
-
-        The native ``xr.infer_freq()`` method does not work for all cases
-        because the frequency can be irregular (e.g., different hour
-        measurements), which ends up returning None.
-
-        Returns
-        -------
-        Frequency
-            The time frequency.
-        """
-        time_coords = self._dataset[self._dim_name]
-        min_delta = pd.to_timedelta(np.diff(time_coords).min(), unit="ns")
-
-        if min_delta < pd.Timedelta(days=1):
-            return "hour"
-        elif min_delta >= pd.Timedelta(days=1) and min_delta < pd.Timedelta(days=28):
-            return "day"
-        elif min_delta >= pd.Timedelta(days=28) and min_delta < pd.Timedelta(days=365):
-            return "month"
-        else:
-            return "year"
+        return self._averager(data_var, "average", freq, weighted, keep_weights)
 
     def group_average(
         self,
         data_var: str,
         freq: Frequency,
         weighted: bool = True,
+        keep_weights: bool = False,
         season_config: SeasonConfigInput = DEFAULT_SEASON_CONFIG,
     ):
         """Returns a Dataset with average of a data variable by time group.
@@ -256,6 +232,9 @@ class TemporalAccessor:
 
             The weight of masked (missing) data is excluded when averages are
             calculated. This is the same as giving them a weight of 0.
+        keep_weights : bool, optional
+            If calculating averages using weights, keep the weights in the
+            final dataset output, by default False.
         season_config: SeasonConfigInput, optional
             A dictionary for "season" frequency configurations. If configs for
             predefined seasons are passed, configs for custom seasons are
@@ -350,13 +329,16 @@ class TemporalAccessor:
             'drop_incomplete_djf': 'False'
         }
         """
-        return self._averager(data_var, "group_average", freq, weighted, season_config)
+        return self._averager(
+            data_var, "group_average", freq, weighted, keep_weights, season_config
+        )
 
     def climatology(
         self,
         data_var: str,
         freq: Frequency,
         weighted: bool = True,
+        keep_weights: bool = False,
         season_config: SeasonConfigInput = DEFAULT_SEASON_CONFIG,
     ):
         """Returns a Dataset with the climatology of a data variable.
@@ -383,6 +365,9 @@ class TemporalAccessor:
 
             The weight of masked (missing) data is excluded when averages are
             taken. This is the same as giving them a weight of 0.
+        keep_weights : bool, optional
+            If calculating averages using weights, keep the weights in the
+            final dataset output, by default False.
         season_config: SeasonConfigInput, optional
             A dictionary for "season" frequency configurations. If configs for
             predefined seasons are passed, configs for custom seasons are
@@ -477,13 +462,16 @@ class TemporalAccessor:
             'drop_incomplete_djf': 'False'
         }
         """
-        return self._averager(data_var, "climatology", freq, weighted, season_config)
+        return self._averager(
+            data_var, "climatology", freq, weighted, keep_weights, season_config
+        )
 
     def departures(
         self,
         data_var: str,
         freq: Frequency,
         weighted: bool = True,
+        keep_weights: bool = False,
         season_config: SeasonConfigInput = DEFAULT_SEASON_CONFIG,
     ) -> xr.Dataset:
         """
@@ -530,6 +518,9 @@ class TemporalAccessor:
 
             The weight of masked (missing) data is excluded when averages are
             taken. This is the same as giving them a weight of 0.
+        keep_weights : bool, optional
+            If calculating averages using weights, keep the weights in the
+            final dataset output, by default False.
         season_config: SeasonConfigInput, optional
             A dictionary for "season" frequency configurations. If configs for
             predefined seasons are passed, configs for custom seasons are
@@ -601,18 +592,18 @@ class TemporalAccessor:
             'drop_incomplete_djf': 'False'
         }
         """
-        self._set_obj_attrs("departures", freq, weighted, season_config)
-
         ds = self._dataset.copy()
+        self._set_obj_attrs("departures", freq, weighted, season_config)
 
         # Group the observation data variable.
         dv_obs = _get_data_var(ds, data_var)
         dv_obs_grouped = self._group_data(dv_obs)
 
         # Calculate the climatology of the data variable.
-        dv_climo = ds.temporal.climatology(data_var, freq, weighted, season_config)[
-            data_var
-        ]
+        ds_climo = ds.temporal.climatology(
+            data_var, freq, weighted, keep_weights, 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
@@ -633,14 +624,47 @@ class TemporalAccessor:
             # are no longer required.
             ds_departs = ds_departs.drop_vars(self._labeled_time.name)
 
+        if weighted and keep_weights:
+            self._weights = ds_climo.time_wts
+            ds_departs = self._keep_weights(ds_departs)
+
         return ds_departs
 
+    def _infer_freq(self) -> Frequency:
+        """Infers the time frequency from the coordinates.
+
+        This method infers the time frequency from the coordinates by
+        calculating the minimum delta and comparing it against a set of
+        conditionals.
+
+        The native ``xr.infer_freq()`` method does not work for all cases
+        because the frequency can be irregular (e.g., different hour
+        measurements), which ends up returning None.
+
+        Returns
+        -------
+        Frequency
+            The time frequency.
+        """
+        time_coords = self._dataset[self._dim_name]
+        min_delta = pd.to_timedelta(np.diff(time_coords).min(), unit="ns")
+
+        if min_delta < pd.Timedelta(days=1):
+            return "hour"
+        elif min_delta >= pd.Timedelta(days=1) and min_delta < pd.Timedelta(days=28):
+            return "day"
+        elif min_delta >= pd.Timedelta(days=28) and min_delta < pd.Timedelta(days=365):
+            return "month"
+        else:
+            return "year"
+
     def _averager(
         self,
         data_var: str,
         mode: Mode,
         freq: Frequency,
         weighted: bool = True,
+        keep_weights: bool = False,
         season_config: SeasonConfigInput = DEFAULT_SEASON_CONFIG,
     ) -> xr.Dataset:
         """Averages a data variable based on the averaging mode and frequency."""
@@ -660,13 +684,16 @@ class TemporalAccessor:
         elif self._mode in ["group_average", "climatology", "departures"]:
             dv = self._group_average(dv)
 
-        # The original time dimension is dropped from the Dataset because
-        # it becomes obsolete after the data variable is averaged. A new time
-        # dimension will be added to the Dataset when adding the averaged
-        # data variable.
+        # The original time dimension is dropped from the dataset because
+        # it becomes obsolete after the data variable is averaged. When the
+        # averaged data variable is added to the dataset, the new time dimension
+        # and its associated coordinates are also added.
         ds = ds.drop_dims("time")
         ds[dv.name] = dv
 
+        if keep_weights:
+            ds = self._keep_weights(ds)
+
         return ds
 
     def _set_obj_attrs(
@@ -714,7 +741,6 @@ class TemporalAccessor:
                 f"include: {list(freq_keys)}."
             )
 
-        self._time_bounds = self._dataset.bounds.get_bounds("time").copy()
         self._mode = mode
         self._freq = freq
         self._weighted = weighted
@@ -878,7 +904,6 @@ class TemporalAccessor:
 
         if self._weighted:
             self._weights = self._get_weights()
-
             dv *= self._weights
             dv = self._group_data(dv).sum()  # type: ignore
         else:
@@ -903,6 +928,86 @@ class TemporalAccessor:
 
         return dv
 
+    def _get_weights(self) -> xr.DataArray:
+        """Calculates weights for a data variable using time bounds.
+
+        This method gets the length of time for each coordinate point by using
+        the difference in the upper and lower time bounds. This approach ensures
+        that the correct time lengths are calculated regardless of how time
+        coordinates are recorded (e.g., monthly, daily, hourly) and the calendar
+        type used.
+
+        The time lengths are labeled and grouped, then each time length is
+        divided by the total sum of the time lengths in its group to get its
+        corresponding weight.
+
+        The sum of the weights for each group is validated to ensure it equals
+        1.0.
+
+        Returns
+        -------
+        xr.DataArray
+            The weights based on a specified frequency.
+
+        Notes
+        -----
+        Refer to [5]_ for the supported CF convention calendar types.
+
+        References
+        ----------
+        .. [5] https://cfconventions.org/cf-conventions/cf-conventions.html#calendar
+        """
+        with xr.set_options(keep_attrs=True):
+            time_lengths: xr.DataArray = (
+                self._time_bounds[:, 1] - self._time_bounds[:, 0]
+            )
+
+        # Must be cast dtype from "timedelta64[ns]" to "float64", specifically
+        # when using Dask arrays. Otherwise, the numpy warning below is thrown:
+        # `DeprecationWarning: The `dtype` and `signature` arguments to ufuncs
+        # only select the general DType and not details such as the byte order
+        # or time unit (with rare exceptions see release notes). To avoid this
+        # warning please use the scalar types `np.float64`, or string notation.`
+        time_lengths = time_lengths.astype(np.float64)
+
+        grouped_time_lengths = self._group_data(time_lengths)
+        weights: xr.DataArray = grouped_time_lengths / grouped_time_lengths.sum()  # type: ignore
+        weights.name = f"{self._dim_name}_wts"
+
+        # Validate the sum of weights for each group is 1.0.
+        actual_sum = self._group_data(weights).sum().values  # type: ignore
+        expected_sum = np.ones(len(grouped_time_lengths.groups))
+        np.testing.assert_allclose(actual_sum, expected_sum)
+
+        return weights
+
+    def _group_data(self, data_var: xr.DataArray) -> DataArrayGroupBy:
+        """Groups a data variable.
+
+        This method groups a data variable by a single datetime component for
+        the "average" mode or labeled time coordinates for all other modes.
+
+        Parameters
+        ----------
+        data_var : xr.DataArray
+            A data variable.
+
+        Returns
+        -------
+        DataArrayGroupBy
+            A data variable grouped by label.
+        """
+        dv = data_var.copy()
+
+        if self._mode == "average":
+            dv_gb = dv.groupby(f"{self._dim_name}.{self._freq}")
+        else:
+            self._labeled_time = self._label_time_coords(dv[self._dim_name])
+            dv.coords[self._labeled_time.name] = self._labeled_time
+            dv_gb = dv.groupby(self._labeled_time.name)
+
+        return dv_gb
+
     def _label_time_coords(self, time_coords: xr.DataArray) -> xr.DataArray:
         """Labels time coordinates with a group for grouping.
 
@@ -1251,97 +1356,40 @@ class TemporalAccessor:
 
         return np.array(dates)
 
-    def _get_weights(self) -> xr.DataArray:
-        """Calculates weights for a data variable using time bounds.
-
-        This method gets the length of time for each coordinate point by using
-        the difference in the upper and lower time bounds. This approach ensures
-        that the correct time lengths are calculated regardless of how time
-        coordinates are recorded (e.g., monthly, daily, hourly) and the calendar
-        type used.
-
-        The time lengths are labeled and grouped, then each time length is
-        divided by the total sum of the time lengths in its group to get its
-        corresponding weight.
-
-        The sum of the weights for each group is validated to ensure it equals
-        1.0.
-
-        Returns
-        -------
-        xr.DataArray
-            The weights based on a specified frequency.
-
-        Notes
-        -----
-        Refer to [5]_ for the supported CF convention calendar types.
-
-        References
-        ----------
-        .. [5] https://cfconventions.org/cf-conventions/cf-conventions.html#calendar
-        """
-        with xr.set_options(keep_attrs=True):
-            time_lengths: xr.DataArray = (
-                self._time_bounds[:, 1] - self._time_bounds[:, 0]
-            )
-
-        # Must be cast dtype from "timedelta64[ns]" to "float64", specifically
-        # when using Dask arrays. Otherwise, the numpy warning below is thrown:
-        # `DeprecationWarning: The `dtype` and `signature` arguments to ufuncs
-        # only select the general DType and not details such as the byte order
-        # or time unit (with rare exceptions see release notes). To avoid this
-        # warning please use the scalar types `np.float64`, or string notation.`
-        time_lengths = time_lengths.astype(np.float64)
-
-        grouped_time_lengths = self._group_data(time_lengths)
-        weights: xr.DataArray = grouped_time_lengths / grouped_time_lengths.sum()  # type: ignore
-
-        num_groups = len(grouped_time_lengths.groups)
-        self._validate_weights(weights, num_groups)
-
-        return weights
-
-    def _validate_weights(self, weights: xr.DataArray, num_groups: int):
-        """Validates that the sum of the weights for each group equals 1.0.
+    def _keep_weights(self, ds: xr.Dataset) -> xr.Dataset:
+        """Keep the weights in the dataset.
 
         Parameters
         ----------
-        weights : xr.DataArray
-            The weights for the time coordinates.
-        num_groups : int
-            The number of groups.
-        """
-        actual_sum = self._group_data(weights).sum().values  # type: ignore
-        expected_sum = np.ones(num_groups)
-
-        np.testing.assert_allclose(actual_sum, expected_sum)
-
-    def _group_data(self, data_var: xr.DataArray) -> DataArrayGroupBy:
-        """Groups a data variable.
-
-        This method groups a data variable by a single datetime component for
-        the "average" mode or labeled time coordinates for all other modes.
-
-        Parameters
-        ----------
-        data_var : xr.DataArray
-            A data variable.
+        ds : xr.Dataset
+            The dataset.
 
         Returns
         -------
-        DataArrayGroupBy
-            A data variable grouped by label.
+        xr.Dataset
+            The dataset with the weights used for averaging.
         """
-        dv = data_var.copy()
+        # Append "_original" to the name of the weights` time coordinates to
+        # avoid conflict with the grouped time coordinates in the Dataset (can
+        # have a different shape).
+        if self._mode in ["group_average", "climatology"]:
+            self._weights = self._weights.rename(
+                {self._dim_name: f"{self._dim_name}_original"}
+            )
+            # 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_name}_original": self._dim_name}
+            )
 
-        if self._mode == "average":
-            dv_gb = dv.groupby(f"{self._dim_name}.{self._freq}")
-        else:
-            self._labeled_time = self._label_time_coords(dv[self._dim_name])
-            dv.coords[self._labeled_time.name] = self._labeled_time
-            dv_gb = dv.groupby(self._labeled_time.name)
+        ds[self._weights.name] = self._weights
 
-        return dv_gb
+        return ds
 
     def _add_operation_attrs(self, data_var: xr.DataArray) -> xr.DataArray:
         """Adds attributes to the data variable describing the operation.


Test patch:
diff --git a/tests/test_spatial.py b/tests/test_spatial.py
index 49e01b8..3209882 100644
--- a/tests/test_spatial.py
+++ b/tests/test_spatial.py
@@ -39,58 +39,115 @@ class TestAverage:
 
     def test_raises_error_if_data_var_not_in_dataset(self):
         with pytest.raises(KeyError):
-            self.ds.spatial.average(
-                "not_a_data_var",
-                axis=["Y", "incorrect_axis"],
-            )
+            self.ds.spatial.average("not_a_data_var", axis=["Y", "incorrect_axis"])
 
-    def test_spatial_average_for_lat_and_lon_region_using_custom_weights(self):
+    def test_raises_error_if_axis_list_contains_unsupported_axis(self):
+        with pytest.raises(ValueError):
+            self.ds.spatial.average("ts", axis=["Y", "incorrect_axis"])
+
+    def test_raises_error_if_lat_axis_does_not_exist(self):
         ds = self.ds.copy()
+        ds.lat.attrs["axis"] = None
+        with pytest.raises(KeyError):
+            ds.spatial.average("ts", axis=["X", "Y"])
 
+    def test_raises_error_if_lon_axis_does_not_exist(self):
+        ds = self.ds.copy()
+        ds.lon.attrs["axis"] = None
+        with pytest.raises(KeyError):
+            ds.spatial.average("ts", axis=["X", "Y"])
+
+    def test_raises_error_if_bounds_type_is_not_a_tuple(self):
+        with pytest.raises(TypeError):
+            self.ds.spatial.average("ts", axis=["Y"], lat_bounds=[1, 1])
+
+        with pytest.raises(TypeError):
+            self.ds.spatial.average("ts", axis=["Y"], lat_bounds="str")
+
+    def test_raises_error_if_there_are_0_elements_in_the_bounds(self):
+        with pytest.raises(ValueError):
+            self.ds.spatial.average("ts", axis=["Y"], lat_bounds=())
+
+    def test_raises_error_if_there_are_more_than_two_elements_in_the_bounds(self):
+        with pytest.raises(ValueError):
+            self.ds.spatial.average("ts", axis=["Y"], lat_bounds=(1, 2, 3))
+
+    def test_raises_error_if_lower_bound_is_not_a_float_or_int(self):
+        with pytest.raises(TypeError):
+            self.ds.spatial.average("ts", axis=["Y"], lat_bounds=("invalid", 1))
+
+    def test_raises_error_if_upper_bound_is_not_a_float_or_int(self):
+        with pytest.raises(TypeError):
+            self.ds.spatial.average("ts", axis=["Y"], lat_bounds=(1, "invalid"))
+
+    def test_raises_error_if_lower_lat_bound_is_larger_than_upper(self):
+        with pytest.raises(ValueError):
+            self.ds.spatial.average("ts", axis=["Y"], lat_bounds=(2, 1))
+
+    def test_does_not_raise_error_if_lower_lon_bound_is_larger_than_upper(self):
+        self.ds.spatial.average("ts", axis=["X"], lon_bounds=(2, 1))
+
+    def test_raises_error_if_lat_axis_is_specified_but_lat_is_not_in_weights_dims(
+        self,
+    ):
         weights = xr.DataArray(
-            data=np.array([[1, 2, 3, 4], [2, 4, 6, 8], [3, 6, 9, 12], [4, 8, 12, 16]]),
-            coords={"lat": ds.lat, "lon": ds.lon},
-            dims=["lat", "lon"],
+            data=np.ones(4), coords={"lon": self.ds.lon}, dims=["lon"]
         )
-        result = ds.spatial.average(
-            axis=["X", "Y"],
-            lat_bounds=(-5.0, 5),
-            lon_bounds=(-170, -120.1),
-            weights=weights,
-            data_var="ts",
+        with pytest.raises(KeyError):
+            self.ds.spatial.average("ts", axis=["X", "Y"], weights=weights)
+
+    def test_raises_error_if_lon_axis_is_specified_but_lon_is_not_in_weights_dims(
+        self,
+    ):
+        weights = xr.DataArray(
+            data=np.ones(4), coords={"lat": self.ds.lat}, dims=["lat"]
         )
+        with pytest.raises(KeyError):
+            self.ds.spatial.average("ts", axis=["X", "Y"], weights=weights)
 
-        expected = self.ds.copy()
-        expected["ts"] = xr.DataArray(
-            data=np.array([2.25, 1.0, 1.0]),
-            coords={"time": expected.time},
-            dims="time",
+    def test_raises_error_if_weights_lat_and_lon_dims_dont_align_with_data_var_dims(
+        self,
+    ):
+        # Get a slice of the dataset to reduce the size of the dimensions for
+        # simpler testing.
+        ds = self.ds.isel(lat=slice(0, 3), lon=slice(0, 3))
+        weights = xr.DataArray(
+            data=np.ones((3, 3)),
+            coords={"lat": ds.lat, "lon": ds.lon},
+            dims=["lat", "lon"],
         )
 
-        assert result.identical(expected)
+        with pytest.raises(ValueError):
+            self.ds.spatial.average("ts", axis=["X", "Y"], weights=weights)
 
-    def test_spatial_average_for_lat_and_lon_region(self):
+    def test_spatial_average_for_lat_region_and_keep_weights(self):
         ds = self.ds.copy()
+
         result = ds.spatial.average(
-            "ts", axis=["X", "Y"], lat_bounds=(-5.0, 5), lon_bounds=(-170, -120.1)
+            "ts", axis=["Y"], lat_bounds=(-5.0, 5), keep_weights=True
         )
 
         expected = self.ds.copy()
         expected["ts"] = xr.DataArray(
-            data=np.array([2.25, 1.0, 1.0]),
-            coords={"time": expected.time},
-            dims="time",
+            data=np.array(
+                [[2.25, 2.25, 2.25, 2.25], [1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0]]
+            ),
+            coords={"time": expected.time, "lon": expected.lon},
+            dims=["time", "lon"],
+        )
+        expected["lat_wts"] = xr.DataArray(
+            name="lat_wts",
+            data=np.array([0.0, 0.08715574, 0.08715574, 0.0]),
+            dims=["lat"],
         )
 
-        assert result.identical(expected)
+        xr.testing.assert_allclose(result, expected)
 
     def test_spatial_average_for_lat_region(self):
         ds = self.ds.copy()
 
         # Specifying axis as a str instead of list of str.
-        result = ds.spatial.average(
-            "ts", axis=["Y"], lat_bounds=(-5.0, 5), lon_bounds=(-170, -120.1)
-        )
+        result = ds.spatial.average("ts", axis=["Y"], lat_bounds=(-5.0, 5))
 
         expected = self.ds.copy()
         expected["ts"] = xr.DataArray(
@@ -104,12 +161,12 @@ class TestAverage:
         assert result.identical(expected)
 
     @requires_dask
-    def test_chunked_spatial_average_for_lat_region(self):
+    def test_spatial_average_for_lat_region_and_keep_weights_with_dask(self):
         ds = self.ds.copy().chunk(2)
 
         # Specifying axis as a str instead of list of str.
         result = ds.spatial.average(
-            "ts", axis=["Y"], lat_bounds=(-5.0, 5), lon_bounds=(-170, -120.1)
+            "ts", axis=["Y"], lat_bounds=(-5.0, 5), keep_weights=True
         )
 
         expected = self.ds.copy()
@@ -120,132 +177,190 @@ class TestAverage:
             coords={"time": expected.time, "lon": expected.lon},
             dims=["time", "lon"],
         )
+        expected["lat_wts"] = xr.DataArray(
+            name="lat_wts",
+            data=np.array([0.0, 0.08715574, 0.08715574, 0.0]),
+            dims=["lat"],
+        )
 
-        assert result.identical(expected)
+        xr.testing.assert_allclose(result, expected)
 
+    def test_spatial_average_for_lat_and_lon_region_and_keep_weights(self):
+        ds = self.ds.copy()
+        result = ds.spatial.average(
+            "ts",
+            axis=["X", "Y"],
+            lat_bounds=(-5.0, 5),
+            lon_bounds=(-170, -120.1),
+            keep_weights=True,
+        )
 
-class TestValidateAxisArg:
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
+        expected = self.ds.copy()
+        expected["ts"] = xr.DataArray(
+            data=np.array([2.25, 1.0, 1.0]),
+            coords={"time": expected.time},
+            dims="time",
+        )
+        expected["lat_lon_wts"] = xr.DataArray(
+            name="ts_weights",
+            data=np.array(
+                [
+                    [0.0, 0.0, 0.0, 0.0],
+                    [0.0, 0.0, 0.0, 0.0],
+                    [0.0, 4.34907156, 4.34907156, 0.0],
+                    [0.0, 0.0, 0.0, 0.0],
+                ]
+            ),
+            dims=["lon", "lat"],
+        )
 
-    def test_raises_error_if_axis_list_contains_unsupported_axis(self):
-        with pytest.raises(ValueError):
-            self.ds.spatial._validate_axis_arg(axis=["Y", "incorrect_axis"])
+        xr.testing.assert_allclose(result, expected)
 
-    def test_raises_error_if_lat_axis_does_not_exist(self):
+    def test_spatial_average_for_lat_and_lon_region_with_custom_weights(self):
         ds = self.ds.copy()
-        ds.lat.attrs["axis"] = None
-        with pytest.raises(KeyError):
-            ds.spatial._validate_axis_arg(axis=["X", "Y"])
 
-    def test_raises_error_if_lon_axis_does_not_exist(self):
-        ds = self.ds.copy()
-        ds.lon.attrs["axis"] = None
-        with pytest.raises(KeyError):
-            ds.spatial._validate_axis_arg(axis=["X", "Y"])
+        weights = xr.DataArray(
+            data=np.array([[1, 2, 3, 4], [2, 4, 6, 8], [3, 6, 9, 12], [4, 8, 12, 16]]),
+            coords={"lat": ds.lat, "lon": ds.lon},
+            dims=["lat", "lon"],
+        )
+        result = ds.spatial.average(
+            axis=["X", "Y"],
+            lat_bounds=(-5.0, 5),
+            lon_bounds=(-170, -120.1),
+            weights=weights,
+            data_var="ts",
+        )
 
+        expected = self.ds.copy()
+        expected["ts"] = xr.DataArray(
+            data=np.array([2.25, 1.0, 1.0]),
+            coords={"time": expected.time},
+            dims="time",
+        )
 
-class TestValidateRegionBounds:
+        assert result.identical(expected)
+
+
+class TestGetWeights:
     @pytest.fixture(autouse=True)
     def setup(self):
         self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
 
-    def test_raises_error_if_bounds_type_is_not_a_tuple(self):
-        with pytest.raises(TypeError):
-            self.ds.spatial._validate_region_bounds("lon", [1, 1])
-
-        with pytest.raises(TypeError):
-            self.ds.spatial._validate_region_bounds("lon", "str")
-
-    def test_raises_error_if_there_are_0_elements_in_the_bounds(self):
-        with pytest.raises(ValueError):
-            self.ds.spatial._validate_region_bounds("lon", ())
-
-    def test_raises_error_if_there_are_more_than_two_elements_in_the_bounds(self):
-        with pytest.raises(ValueError):
-            self.ds.spatial._validate_region_bounds("lon", (1, 1, 2))
-
-    def test_does_not_raise_error_if_lower_and_upper_bounds_are_floats_or_ints(self):
-        self.ds.spatial._validate_region_bounds("lon", (1, 1))
-        self.ds.spatial._validate_region_bounds("lon", (1, 1.2))
-
-    def test_raises_error_if_lower_bound_is_not_a_float_or_int(self):
-        with pytest.raises(TypeError):
-            self.ds.spatial._validate_region_bounds("Y", ("invalid", 1))
-
-    def test_raises_error_if_upper_bound_is_not_a_float_or_int(self):
-        with pytest.raises(TypeError):
-            self.ds.spatial._validate_region_bounds("X", (1, "invalid"))
+    def test_raises_error_if_domain_lower_bound_exceeds_upper_bound(self):
+        ds = self.ds.copy()
+        ds.lon_bnds.data[:] = np.array([[1, 0], [1, 2], [2, 3], [3, 4]])
 
-    def test_raises_error_if_lower_lat_bound_is_bigger_than_upper(self):
         with pytest.raises(ValueError):
-            self.ds.spatial._validate_region_bounds("Y", (2, 1))
+            ds.spatial.get_weights(axis=["X"])
 
-    def test_does_not_raise_error_if_lon_lower_bound_is_larger_than_upper(self):
-        self.ds.spatial._validate_region_bounds("X", (2, 1))
+    def test_weights_for_region_in_lat_and_lon_domains(self):
+        result = self.ds.spatial.get_weights(
+            axis=["Y", "X"], lat_bounds=(-5, 5), lon_bounds=(-170, -120)
+        )
+        expected = xr.DataArray(
+            data=np.array(
+                [
+                    [0.0, 0.0, 0.0, 0.0],
+                    [0.0, 0.0, 4.35778714, 0.0],
+                    [0.0, 0.0, 4.35778714, 0.0],
+                    [0.0, 0.0, 0.0, 0.0],
+                ]
+            ),
+            coords={"lat": self.ds.lat, "lon": self.ds.lon},
+            dims=["lat", "lon"],
+        )
 
+        xr.testing.assert_allclose(result, expected)
 
-class TestValidateWeights:
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-        self.weights = xr.DataArray(
-            data=np.ones((4, 4)),
+    def test_area_weights_for_region_in_lat_domain(self):
+        result = self.ds.spatial.get_weights(
+            axis=["Y", "X"], lat_bounds=(-5, 5), lon_bounds=None
+        )
+        expected = xr.DataArray(
+            data=np.array(
+                [
+                    [0.0, 0.0, 0.0, 0.0],
+                    [0.16341702, 15.52461668, 15.52461668, 0.16341702],
+                    [0.16341702, 15.52461668, 15.52461668, 0.16341702],
+                    [0.0, 0.0, 0.0, 0.0],
+                ]
+            ),
             coords={"lat": self.ds.lat, "lon": self.ds.lon},
             dims=["lat", "lon"],
         )
 
-    def test_no_error_is_raised_when_spatial_dim_sizes_align_between_weights_and_data_var(
-        self,
-    ):
-        weights = xr.DataArray(
-            data=np.ones((4, 4)),
+        xr.testing.assert_allclose(result, expected)
+
+    def test_weights_for_region_in_lon_domain(self):
+        expected = xr.DataArray(
+            data=np.array(
+                [
+                    [0.0, 0.0, 0.00297475, 0.0],
+                    [0.0, 0.0, 49.99702525, 0.0],
+                    [0.0, 0.0, 49.99702525, 0.0],
+                    [0.0, 0.0, 0.00297475, 0.0],
+                ]
+            ),
             coords={"lat": self.ds.lat, "lon": self.ds.lon},
             dims=["lat", "lon"],
         )
-        self.ds.spatial._validate_weights(self.ds["ts"], axis=["Y"], weights=weights)
+        result = self.ds.spatial.get_weights(
+            axis=["Y", "X"], lat_bounds=None, lon_bounds=(-170, -120)
+        )
 
-    def test_error_is_raised_when_lat_axis_is_specified_but_lat_is_not_in_weights_dims(
-        self,
-    ):
-        weights = xr.DataArray(
-            data=np.ones(4), coords={"lon": self.ds.lon}, dims=["lon"]
+        xr.testing.assert_allclose(result, expected)
+
+    def test_weights_for_region_in_lon_domain_with_region_spanning_p_meridian(self):
+        ds = self.ds.copy()
+
+        result = ds.spatial._get_longitude_weights(
+            domain_bounds=ds.lon_bnds,
+            # Region spans prime meridian.
+            region_bounds=np.array([359, 1]),
+        )
+        expected = xr.DataArray(
+            data=np.array([1.875, 0.0625, 0.0, 0.0625]),
+            coords={"lon": ds.lon},
+            dims=["lon"],
         )
-        with pytest.raises(KeyError):
-            self.ds.spatial._validate_weights(
-                self.ds["ts"], axis=["X", "Y"], weights=weights
-            )
 
-    def test_error_is_raised_when_lon_axis_is_specified_but_lon_is_not_in_weights_dims(
-        self,
-    ):
-        weights = xr.DataArray(
-            data=np.ones(4), coords={"lat": self.ds.lat}, dims=["lat"]
+        xr.testing.assert_allclose(result, expected)
+
+    def test_weights_all_longitudes_for_equal_region_bounds(self):
+        expected = xr.DataArray(
+            data=np.array(
+                [1.875, 178.125, 178.125, 1.875],
+            ),
+            coords={"lon": self.ds.lon},
+            dims=["lon"],
+        )
+        result = self.ds.spatial.get_weights(
+            axis=["X"],
+            lat_bounds=None,
+            lon_bounds=np.array([0.0, 360.0]),
         )
-        with pytest.raises(KeyError):
-            self.ds.spatial._validate_weights(
-                self.ds["ts"], axis=["X", "Y"], weights=weights
-            )
 
-    def test_error_is_raised_when_weights_lat_and_lon_dims_dont_align_with_data_var_dims(
-        self,
-    ):
-        # Get a slice of the dataset to reduce the size of the dimensions for
-        # simpler testing.
-        ds = self.ds.isel(lat=slice(0, 3), lon=slice(0, 3))
-        weights = xr.DataArray(
-            data=np.ones((3, 3)),
-            coords={"lat": ds.lat, "lon": ds.lon},
-            dims=["lat", "lon"],
+        xr.testing.assert_allclose(result, expected)
+
+    def test_weights_for_equal_region_bounds_representing_entire_lon_domain(self):
+        expected = xr.DataArray(
+            data=np.array(
+                [1.875, 178.125, 178.125, 1.875],
+            ),
+            coords={"lon": self.ds.lon},
+            dims=["lon"],
+        )
+        result = self.ds.spatial.get_weights(
+            axis=["X"], lat_bounds=None, lon_bounds=np.array([10.0, 10.0])
         )
 
-        with pytest.raises(ValueError):
-            self.ds.spatial._validate_weights(
-                self.ds["ts"], axis=["X", "Y"], weights=weights
-            )
+        xr.testing.assert_allclose(result, expected)
 
 
-class TestSwapLonAxis:
+class Test_SwapLonAxis:
+    # NOTE: This private method is tested because we might want to consider
+    # converting it to a public method in the future.
     @pytest.fixture(autouse=True)
     def setup(self):
         self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
@@ -391,233 +506,9 @@ class TestSwapLonAxis:
         assert np.array_equal(result, expected)
 
 
-class TestGetWeights:
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-
-    def test_weights_for_region_in_lat_and_lon_domains(self):
-        result = self.ds.spatial._get_weights(
-            axis=["Y", "X"], lat_bounds=(-5, 5), lon_bounds=(-170, -120)
-        )
-        expected = xr.DataArray(
-            data=np.array(
-                [
-                    [0.0, 0.0, 0.0, 0.0],
-                    [0.0, 0.0, 4.35778714, 0.0],
-                    [0.0, 0.0, 4.35778714, 0.0],
-                    [0.0, 0.0, 0.0, 0.0],
-                ]
-            ),
-            coords={"lat": self.ds.lat, "lon": self.ds.lon},
-            dims=["lat", "lon"],
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-    def test_area_weights_for_region_in_lat_domain(self):
-        result = self.ds.spatial._get_weights(
-            axis=["Y", "X"], lat_bounds=(-5, 5), lon_bounds=None
-        )
-        expected = xr.DataArray(
-            data=np.array(
-                [
-                    [0.0, 0.0, 0.0, 0.0],
-                    [0.16341702, 15.52461668, 15.52461668, 0.16341702],
-                    [0.16341702, 15.52461668, 15.52461668, 0.16341702],
-                    [0.0, 0.0, 0.0, 0.0],
-                ]
-            ),
-            coords={"lat": self.ds.lat, "lon": self.ds.lon},
-            dims=["lat", "lon"],
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-    def test_weights_for_region_in_lon_domain(self):
-        expected = xr.DataArray(
-            data=np.array(
-                [
-                    [0.0, 0.0, 0.00297475, 0.0],
-                    [0.0, 0.0, 49.99702525, 0.0],
-                    [0.0, 0.0, 49.99702525, 0.0],
-                    [0.0, 0.0, 0.00297475, 0.0],
-                ]
-            ),
-            coords={"lat": self.ds.lat, "lon": self.ds.lon},
-            dims=["lat", "lon"],
-        )
-        result = self.ds.spatial._get_weights(
-            axis=["Y", "X"], lat_bounds=None, lon_bounds=(-170, -120)
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-
-class TestGetLongitudeWeights:
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-
-    def test_weights_for_region_in_lon_domain(self):
-        # Longitude axis orientation swaps from (-180, 180) to (0, 360).
-        result = self.ds.spatial._get_longitude_weights(
-            domain_bounds=self.ds.lon_bnds.copy(),
-            region_bounds=np.array([-170.0, -120.0]),
-        )
-        expected = xr.DataArray(
-            data=np.array([0.0, 0.0, 50.0, 0.0]),
-            coords={"lon": self.ds.lon},
-            dims=["lon"],
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-    def test_weights_for_region_in_lon_domain_with_both_spanning_p_meridian(self):
-        ds = self.ds.copy()
-        # Domain spans prime meridian.
-        ds.lon_bnds.data[:] = np.array([[359, 1], [1, 90], [90, 180], [180, 359]])
-
-        result = ds.spatial._get_longitude_weights(
-            domain_bounds=ds.lon_bnds,
-            # Region spans prime meridian.
-            region_bounds=np.array([359, 1]),
-        )
-        expected = xr.DataArray(
-            data=np.array([2.0, 0.0, 0.0, 0.0]),
-            coords={"lon": ds.lon},
-            dims=["lon"],
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-    def test_weights_for_region_in_lon_domain_with_domain_spanning_p_meridian(self):
-        ds = self.ds.copy()
-        # Domain spans prime meridian.
-        ds.lon_bnds.data[:] = np.array([[359, 1], [1, 90], [90, 180], [180, 359]])
-
-        # Longitude axis orientation swaps from (-180, 180) to (0, 360).
-        result = ds.spatial._get_longitude_weights(
-            domain_bounds=ds.lon_bnds,
-            region_bounds=np.array([-170.0, -120.0]),
-        )
-        expected = xr.DataArray(
-            data=np.array([0.0, 0.0, 0.0, 50.0]),
-            coords={"lon": ds.lon},
-            dims=["lon"],
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-    def test_weights_for_region_in_lon_domain_with_region_spanning_p_meridian(self):
-        ds = self.ds.copy()
-
-        result = ds.spatial._get_longitude_weights(
-            domain_bounds=ds.lon_bnds,
-            # Region spans prime meridian.
-            region_bounds=np.array([359, 1]),
-        )
-        expected = xr.DataArray(
-            data=np.array([1.875, 0.0625, 0.0, 0.0625]),
-            coords={"lon": ds.lon},
-            dims=["lon"],
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-    def test_weights_all_longitudes_for_equal_region_bounds(self):
-        expected = xr.DataArray(
-            data=np.array(
-                [1.875, 178.125, 178.125, 1.875],
-            ),
-            coords={"lon": self.ds.lon},
-            dims=["lon"],
-        )
-        result = self.ds.spatial._get_longitude_weights(
-            domain_bounds=self.ds.lon_bnds.copy(),
-            region_bounds=np.array([0.0, 360.0]),
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-    def test_weights_for_equal_region_bounds_representing_entire_lon_domain(self):
-        expected = xr.DataArray(
-            data=np.array(
-                [1.875, 178.125, 178.125, 1.875],
-            ),
-            coords={"lon": self.ds.lon},
-            dims=["lon"],
-        )
-        result = self.ds.spatial._get_longitude_weights(
-            domain_bounds=self.ds.lon_bnds.copy(), region_bounds=np.array([10.0, 10.0])
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-
-class TestGetLatitudeWeights:
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-
-    def test_weights_for_region_in_lat_domain(self):
-        expected = xr.DataArray(
-            data=np.array([0.0, 0.087156, 0.087156, 0.0]),
-            coords={"lat": self.ds.lat},
-            dims=["lat"],
-        )
-        result = self.ds.spatial._get_latitude_weights(
-            domain_bounds=self.ds.lat_bnds, region_bounds=np.array([-5.0, 5.0])
-        )
-
-        xr.testing.assert_allclose(result, expected)
-
-
-class TestValidateDomainBounds:
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-
-    def test_raises_error_if_low_bounds_exceeds_high_bound(self):
-        domain_bounds = xr.DataArray(
-            name="lon_bnds",
-            data=np.array([[1, 0], [1, 2], [2, 3], [3, 4]]),
-            dims=["lon", "bnds"],
-        )
-        with pytest.raises(ValueError):
-            self.ds.spatial._validate_domain_bounds(domain_bounds)
-
-
-class TestCalculateWeights:
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-
-    def test_returns_weights_as_the_absolute_difference_of_upper_and_lower_bounds(self):
-        lat = xr.DataArray(
-            name="lat",
-            data=np.array([-90.0, -88.75, 88.75, 90.0]),
-            coords={"lat": np.array([-90.0, -88.75, 88.75, 90.0])},
-            dims=["lat"],
-        )
-        lat_bounds = xr.DataArray(
-            data=np.array(
-                [[-90.0, -89.375], [-89.375, 0.0], [0.0, 89.375], [89.375, 90.0]]
-            ),
-            coords={"lat": lat},
-            dims=["lat", "bnds"],
-        )
-
-        result = self.ds.spatial._calculate_weights(lat_bounds)
-        expected = xr.DataArray(
-            data=np.array([0.625, 89.375, 89.375, 0.625]),
-            coords={"lat": lat},
-            dims=["lat"],
-        )
-        assert result.identical(expected)
-
-
-class TestScaleDimToRegion:
+class Test_ScaleDimToRegion:
+    # NOTE: This private method is tested because this is an in-house algorithm
+    # that has edge cases with some complexities.
     @pytest.fixture(autouse=True)
     def setup(self):
         self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
@@ -766,115 +657,3 @@ class TestScaleDimToRegion:
         )
 
         assert result.identical(expected)
-
-
-class TestCombineWeights:
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-        self.axis_weights = {
-            "lat": xr.DataArray(
-                name="lat_wts",
-                data=np.array([1, 2, 3, 4]),
-                coords={"lat": self.ds.lat},
-                dims=["lat"],
-            ),
-            "lon": xr.DataArray(
-                name="lon_wts",
-                data=np.array([1, 2, 3, 4]),
-                coords={"lon": self.ds.lon},
-                dims=["lon"],
-            ),
-        }
-
-    def test_weights_for_single_axis_are_identical(self):
-        axis_weights = self.axis_weights
-        del axis_weights["lon"]
-
-        result = self.ds.spatial._combine_weights(axis_weights=self.axis_weights)
-        expected = self.axis_weights["lat"]
-
-        assert result.identical(expected)
-
-    def test_weights_for_multiple_axis_is_the_product_of_matrix_multiplication(self):
-        result = self.ds.spatial._combine_weights(axis_weights=self.axis_weights)
-        expected = xr.DataArray(
-            data=np.array([[1, 2, 3, 4], [2, 4, 6, 8], [3, 6, 9, 12], [4, 8, 12, 16]]),
-            coords={"lat": self.ds.lat, "lon": self.ds.lon},
-            dims=["lat", "lon"],
-        )
-
-        assert result.identical(expected)
-
-
-class TestAverager:
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-
-    @requires_dask
-    def test_chunked_weighted_avg_over_lat_and_lon_axes(self):
-        ds = self.ds.copy().chunk(2)
-
-        weights = xr.DataArray(
-            data=np.array([[1, 2, 3, 4], [2, 4, 6, 8], [3, 6, 9, 12], [4, 8, 12, 16]]),
-            coords={"lat": ds.lat, "lon": ds.lon},
-            dims=["lat", "lon"],
-        )
-
-        result = ds.spatial._averager(ds.ts, axis=["X", "Y"], weights=weights)
-        expected = xr.DataArray(
-            name="ts", data=np.ones(15), coords={"time": ds.time}, dims=["time"]
-        )
-
-        assert result.identical(expected)
-
-    def test_weighted_avg_over_lat_axis(self):
-        weights = xr.DataArray(
-            name="lat_wts",
-            data=np.array([1, 2, 3, 4]),
-            coords={"lat": self.ds.lat},
-            dims=["lat"],
-        )
-
-        result = self.ds.spatial._averager(self.ds.ts, axis=["Y"], weights=weights)
-        expected = xr.DataArray(
-            name="ts",
-            data=np.ones((15, 4)),
-            coords={"time": self.ds.time, "lon": self.ds.lon},
-            dims=["time", "lon"],
-        )
-
-        assert result.identical(expected)
-
-    def test_weighted_avg_over_lon_axis(self):
-        weights = xr.DataArray(
-            name="lon_wts",
-            data=np.array([1, 2, 3, 4]),
-            coords={"lon": self.ds.lon},
-            dims=["lon"],
-        )
-
-        result = self.ds.spatial._averager(self.ds.ts, axis=["X"], weights=weights)
-        expected = xr.DataArray(
-            name="ts",
-            data=np.ones((15, 4)),
-            coords={"time": self.ds.time, "lat": self.ds.lat},
-            dims=["time", "lat"],
-        )
-
-        assert result.identical(expected)
-
-    def test_weighted_avg_over_lat_and_lon_axis(self):
-        weights = xr.DataArray(
-            data=np.array([[1, 2, 3, 4], [2, 4, 6, 8], [3, 6, 9, 12], [4, 8, 12, 16]]),
-            coords={"lat": self.ds.lat, "lon": self.ds.lon},
-            dims=["lat", "lon"],
-        )
-
-        result = self.ds.spatial._averager(self.ds.ts, axis=["X", "Y"], weights=weights)
-        expected = xr.DataArray(
-            name="ts", data=np.ones(15), coords={"time": self.ds.time}, dims=["time"]
-        )
-
-        assert result.identical(expected)
diff --git a/tests/test_temporal.py b/tests/test_temporal.py
index 73ff885..4a853be 100644
--- a/tests/test_temporal.py
+++ b/tests/test_temporal.py
@@ -1421,8 +1421,8 @@ class TestDepartures:
         self.seasons = ["JJA", "MAM", "SON", "DJF"]
 
     def test_weighted_seasonal_departures_with_DJF(self):
-        # Create a post-climatology dataset.
         ds = self.ds.copy()
+
         # Drop incomplete DJF seasons
         ds = ds.isel(time=slice(2, -1))
 
@@ -1453,6 +1453,62 @@ class TestDepartures:
 
         assert result.identical(expected)
 
+    def test_weighted_seasonal_departures_with_DJF_and_keep_weights(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=True,
+            keep_weights=True,
+            season_config={"dec_mode": "DJF", "drop_incomplete_djf": True},
+        )
+        expected = ds.copy()
+        expected["ts"] = xr.DataArray(
+            data=np.zeros((12, 4, 4)),
+            coords={
+                "lat": expected.lat,
+                "lon": expected.lon,
+                "time": ds.time,
+            },
+            dims=["time", "lat", "lon"],
+            attrs={
+                "operation": "temporal_avg",
+                "mode": "departures",
+                "freq": "season",
+                "weighted": "True",
+                "dec_mode": "DJF",
+                "drop_incomplete_djf": "True",
+            },
+        )
+        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"],
+        )
+
+        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()
         # Drop incomplete DJF seasons
@@ -1517,133 +1573,7 @@ class TestDepartures:
         assert result.identical(expected)
 
 
-class Test_SetObjAttrs:
-    # NOTE: Testing this private method directly instead of through the public
-    # methods because it eliminates redundancy.
-    @pytest.fixture(autouse=True)
-    def setup(self):
-        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
-
-    def test_raises_error_if_operation_is_not_supported(self):
-        with pytest.raises(ValueError):
-            self.ds.temporal._set_obj_attrs(
-                "unsupported",
-                freq="season",
-                weighted=True,
-                season_config={
-                    "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
-                    "custom_seasons": None,
-                },
-            )
-
-    def test_raises_error_if_freq_arg_is_not_supported_by_operation(self):
-        ds = self.ds.copy()
-
-        with pytest.raises(ValueError):
-            ds.temporal._set_obj_attrs(
-                "group_average",
-                freq="unsupported",
-                weighted=True,
-                season_config={
-                    "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
-                    "custom_seasons": None,
-                },
-            )
-        with pytest.raises(ValueError):
-            ds.temporal._set_obj_attrs(
-                "climatology",
-                freq="unsupported",
-                weighted=True,
-                season_config={
-                    "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
-                    "custom_seasons": None,
-                },
-            )
-        with pytest.raises(ValueError):
-            ds.temporal._set_obj_attrs(
-                "departures",
-                freq="unsupported",
-                weighted=True,
-                season_config={
-                    "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
-                    "custom_seasons": None,
-                },
-            )
-
-    def test_does_not_raise_error_if_freq_arg_is_supported_by_operation(self):
-        ds = self.ds.copy()
-        climatology_freqs = ["season", "month", "day"]
-        departure_freqs = ["season", "month", "day"]
-        time_series_freqs = ["year", "season", "month", "day", "hour"]
-
-        for freq in time_series_freqs:
-            ds.temporal._set_obj_attrs(
-                "group_average",
-                freq=freq,
-                weighted=True,
-                season_config={
-                    "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
-                    "custom_seasons": None,
-                },
-            )
-
-        for freq in climatology_freqs:
-            ds.temporal._set_obj_attrs(
-                "climatology",
-                freq=freq,
-                weighted=True,
-                season_config={
-                    "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
-                    "custom_seasons": None,
-                },
-            )
-
-        for freq in departure_freqs:
-            ds.temporal._set_obj_attrs(
-                "departures",
-                freq=freq,
-                weighted=True,
-                season_config={
-                    "dec_mode": "DJF",
-                    "drop_incomplete_djf": False,
-                    "custom_seasons": None,
-                },
-            )
-
-    def test_raises_error_if_season_config_key_is_not_supported(self):
-        with pytest.raises(KeyError):
-            self.ds.temporal._set_obj_attrs(
-                "climatology",
-                freq="season",
-                weighted=True,
-                season_config={
-                    "not_supported": "invalid",
-                },
-            )
-
-    def test_raises_error_if_december_mode_is_not_supported(self):
-        with pytest.raises(ValueError):
-            self.ds.temporal._set_obj_attrs(
-                "climatology",
-                freq="season",
-                weighted=True,
-                season_config={
-                    "dec_mode": "unsupported",
-                    "drop_incomplete_djf": False,
-                    "custom_seasons": None,
-                },
-            )
-
-
 class Test_GetWeights:
-    # NOTE: Testing this private method directly instead of through the public
-    # methods because there is potential for this method to become public.
     class TestWeightsForAverageMode:
         @pytest.fixture(autouse=True)
         def setup(self):
@@ -1653,7 +1583,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "average"
             ds.temporal._freq = "year"
             ds.temporal._weighted = "True"
@@ -1716,7 +1645,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "average"
             ds.temporal._freq = "month"
             ds.temporal._weighted = "True"
@@ -1783,7 +1711,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "group_average"
             ds.temporal._freq = "year"
             ds.temporal._weighted = "True"
@@ -1846,7 +1773,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "group_average"
             ds.temporal._freq = "month"
             ds.temporal._weighted = "True"
@@ -1917,6 +1843,7 @@ class Test_GetWeights:
                     "axis": "T",
                     "long_name": "time",
                     "standard_name": "time",
+                    "bounds": "time_bnds",
                 },
             )
             ds["ts"] = xr.DataArray(
@@ -1988,7 +1915,6 @@ class Test_GetWeights:
             )
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "group_average"
             ds.temporal._freq = "season"
             ds.temporal._weighted = "True"
@@ -2050,7 +1976,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "group_average"
             ds.temporal._freq = "season"
             ds.temporal._weighted = "True"
@@ -2120,7 +2045,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "group_average"
             ds.temporal._freq = "season"
             ds.temporal._weighted = "True"
@@ -2197,7 +2121,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "group_average"
             ds.temporal._freq = "day"
             ds.temporal._weighted = "True"
@@ -2242,7 +2165,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "group_average"
             ds.temporal._freq = "hour"
             ds.temporal._weighted = "True"
@@ -2317,6 +2239,7 @@ class Test_GetWeights:
                     "axis": "T",
                     "long_name": "time",
                     "standard_name": "time",
+                    "bounds": "time_bnds",
                 },
             )
             ds["ts"] = xr.DataArray(
@@ -2388,7 +2311,6 @@ class Test_GetWeights:
             )
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "climatology"
             ds.temporal._freq = "season"
             ds.temporal._weighted = "True"
@@ -2445,7 +2367,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "climatology"
             ds.temporal._freq = "season"
             ds.temporal._weighted = "True"
@@ -2510,7 +2431,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "climatology"
             ds.temporal._freq = "month"
             ds.temporal._weighted = "True"
@@ -2573,7 +2493,6 @@ class Test_GetWeights:
             ds = self.ds.copy()
 
             # Set object attrs required to test the method.
-            ds.temporal._time_bounds = ds.time_bnds.copy()
             ds.temporal._mode = "climatology"
             ds.temporal._freq = "day"
             ds.temporal._weighted = "True"
@@ -2629,3 +2548,91 @@ class Test_GetWeights:
                 ]
             )
             assert np.allclose(result, expected)
+
+
+class Test_Averager:
+    # NOTE: This private method is tested because it is more redundant to
+    # test these cases for the public methods that call this private method.
+    @pytest.fixture(autouse=True)
+    def setup(self):
+        self.ds = generate_dataset(cf_compliant=True, has_bounds=True)
+
+    def test_raises_error_with_incorrect_mode_arg(self):
+        with pytest.raises(ValueError):
+            self.ds.temporal._averager(
+                "ts",
+                "unsupported",
+                freq="season",
+                weighted=True,
+                season_config={
+                    "dec_mode": "DJF",
+                    "drop_incomplete_djf": False,
+                    "custom_seasons": None,
+                },
+            )
+
+    def test_raises_error_if_freq_arg_is_not_supported_by_operation(self):
+        ds = self.ds.copy()
+
+        with pytest.raises(ValueError):
+            ds.temporal._averager(
+                "ts",
+                "group_average",
+                freq="unsupported",
+                weighted=True,
+                season_config={
+                    "dec_mode": "DJF",
+                    "drop_incomplete_djf": False,
+                    "custom_seasons": None,
+                },
+            )
+        with pytest.raises(ValueError):
+            ds.temporal._averager(
+                "ts",
+                "climatology",
+                freq="unsupported",
+                weighted=True,
+                season_config={
+                    "dec_mode": "DJF",
+                    "drop_incomplete_djf": False,
+                    "custom_seasons": None,
+                },
+            )
+        with pytest.raises(ValueError):
+            ds.temporal._averager(
+                "ts",
+                "departures",
+                freq="unsupported",
+                weighted=True,
+                season_config={
+                    "dec_mode": "DJF",
+                    "drop_incomplete_djf": False,
+                    "custom_seasons": None,
+                },
+            )
+
+    def test_raises_error_if_season_config_key_is_not_supported(self):
+        with pytest.raises(KeyError):
+            self.ds.temporal._averager(
+                "ts",
+                "climatology",
+                freq="season",
+                weighted=True,
+                season_config={
+                    "not_supported": "invalid",
+                },
+            )
+
+    def test_raises_error_if_december_mode_is_not_supported(self):
+        with pytest.raises(ValueError):
+            self.ds.temporal._averager(
+                "ts",
+                "climatology",
+                freq="season",
+                weighted=True,
+                season_config={
+                    "dec_mode": "unsupported",
+                    "drop_incomplete_djf": False,
+                    "custom_seasons": None,
+                },
+            )