Repository: xCDAT/xcdat

Issue:
[Bug]: incorrect group averaging with missing data
### What happened?

I took the annual average of some data. In one example grid cell, the values were:

> [  nan -0.02 -0.02   nan   nan   nan   nan   nan   nan   nan   nan   nan] 

Since two months have the same value and data is missing for all other months, the weighted average should be `-0.02`. Instead, I get `-0.0032786885`; `group_average` seems to be incorrectly excluding the missing data.

### What did you expect to happen?

I expected to get `-0.02`.

### Minimal Complete Verifiable Example

Load some data for an example grid cell:

```python
import xcdat
# open dataset
fn = '/p/user_pub/climate_work/pochedley1/surface/gistemp1200_GHCNv4_ERSSTv5.nc' 
ds = xcdat.open_dataset(fn)
# fix missing calendar
time_encoding = ds.time.encoding
time_encoding['calendar'] = 'standard'
ds.time.encoding = time_encoding
# select grid cell in subsetted dataset
dss = ds.isel(lat=[11], lon=[0])
```

Note that first year of data has two months (of twelve) with data values

```
print(dss.tempanomaly[0:12, 0, 0].values)
```

> [  nan -0.02 -0.02   nan   nan   nan   nan   nan   nan   nan   nan   nan] 

Taking the annual average I get an unexpected result:

```

ds_yearly = dss.temporal.group_average('tempanomaly', freq='year', weighted=True)

print(ds_yearly.tempanomaly[0, 0, 0].values)
```

> -0.0032786885

This is not the right value. Since there are two values in the first year (and they are the same), then the average should simply be `-0.02`. 

I think I see what is happening. Each month is assigned some weight in the first year (proportional to the number of days in each month): 

```
w = dss.temporal._get_weights()
```

If we consider just the first year, then:

```
w1 = w[0:12]
t1 = dss.tempanomaly[0:12, 0, 0]
print(w1.values, t1.values)
```

> [0.08469945 0.07923497 0.08469945 0.08196721 0.08469945 0.08196721                                                                             
>  0.08469945 0.08469945 0.08196721 0.08469945 0.08196721 0.08469945]
> 
> [  nan -0.02 -0.02   nan   nan   nan   nan   nan   nan   nan   nan   nan]

A simple weighted average is `WA = sum(T*W)/sum(W)`

```
print((np.sum(w1*t1)/np.sum(w1)).values)
```

> -0.0032786884513057645 

The problem is that there should be no weight assigned if there is no data for a given month/index. The weights should be corrected to reflect this: 

> [0.0 0.07923497 0.08469945 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

This weighting matrix would yield `-0.02` (which is the correct answer).


### Relevant log output

_No response_

### Anything else we need to know?

I think this is probably handled correctly in [`_average`](https://github.com/xCDAT/xcdat/blob/main/xcdat/temporal.py#L897) (via xarray `.mean`), but [`_group_average`](https://github.com/xCDAT/xcdat/blob/main/xcdat/temporal.py#L923) is calling `.sum` instead of `.mean`. 

### Environment

main branch

Discussion:
codecov-commenter: # [Codecov](https://codecov.io/gh/xCDAT/xcdat/pull/320?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) Report
> Merging [#320](https://codecov.io/gh/xCDAT/xcdat/pull/320?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) (428675f) into [main](https://codecov.io/gh/xCDAT/xcdat/commit/afd7a8215b307c27f6cbbdcfb84cf55f3ae561b2?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=xCDAT) (afd7a82) will **not change** coverage.
> The diff coverage is `100.00%`.

```diff
@@            Coverage Diff            @@
##              main      #320   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           14        14           
  Lines         1191      1201   +10     
=========================================
+ Hits          1191      1201   +10     
```


| [Impacted Files](https://codecov.io/gh/xCDAT/xcdat/pull/320?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/320/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%> (ø)` | |

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)

tomvothecoder: Thanks for opening this PR. FYI, I realized I won't have time to review this until next week.

Patch:
diff --git a/xcdat/temporal.py b/xcdat/temporal.py
index 4342d9d..7142c23 100644
--- a/xcdat/temporal.py
+++ b/xcdat/temporal.py
@@ -921,8 +921,23 @@ class TemporalAccessor:
 
         if self._weighted:
             self._weights = self._get_weights()
+            # Weight the data variable.
             dv *= self._weights
-            dv = self._group_data(dv).sum()
+
+            # Ensure missing data (`np.nan`) receives no weight (zero). To
+            # achieve this, first broadcast the one-dimensional (temporal
+            # dimension) shape of the `weights` DataArray to the
+            # multi-dimensional shape of its corresponding data variable.
+            weights, _ = xr.broadcast(self._weights, dv)
+            weights = xr.where(np.isnan(dv), 0.0, weights)
+
+            # Perform weighted average using the formula
+            # WA = sum(data*weights) / sum(weights). The denominator must be
+            # included to take into account zero weight for missing data.
+            with xr.set_options(keep_attrs=True):
+                dv = self._group_data(dv).sum() / self._group_data(weights).sum()
+            # Restore the data variable's name.
+            dv.name = data_var.name
         else:
             dv = self._group_data(dv).mean()
 
@@ -931,7 +946,7 @@ class TemporalAccessor:
         # For example, grouping on "year_season" replaces the time dimension
         # with "year_season". This dimension needs to be renamed back to
         # the original time dimension name before the data variable is added
-        # back to the dataset so that the CF compliant name is maintained.
+        # back to the dataset so that the original name is preserved.
         dv = dv.rename({self._labeled_time.name: self._dim})
 
         # After grouping and aggregating, the grouped time dimension's
@@ -1021,7 +1036,12 @@ class TemporalAccessor:
         if self._mode == "average":
             dv_gb = dv.groupby(f"{self._dim}.{self._freq}")
         else:
-            self._labeled_time = self._label_time_coords(dv[self._dim])
+            # Reuse the labeled time coordinates if they are already generated.
+            # This specifically applies when generating weights, which calls
+            # this method again after it has already been called to group
+            # the data variable with the same time coordinates.
+            if not hasattr(self, "_labeled_time"):
+                self._labeled_time = self._label_time_coords(dv[self._dim])
             dv.coords[self._labeled_time.name] = self._labeled_time
             dv_gb = dv.groupby(self._labeled_time.name)
 


Test patch:
diff --git a/tests/test_temporal.py b/tests/test_temporal.py
index ffa49f6..218eef8 100644
--- a/tests/test_temporal.py
+++ b/tests/test_temporal.py
@@ -418,6 +418,7 @@ class TestGroupAverage:
             ),
             coords={"time": self.ds.time, "lat": self.ds.lat, "lon": self.ds.lon},
             dims=["time", "lat", "lon"],
+            attrs={"test_attr": "test"},
         )
 
     def test_weighted_annual_averages(self):
@@ -458,6 +459,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "year",
@@ -467,6 +469,7 @@ class TestGroupAverage:
 
         xr.testing.assert_allclose(result, expected)
         assert result.ts.attrs == expected.ts.attrs
+        assert result.time.attrs == expected.time.attrs
 
     @requires_dask
     def test_weighted_annual_averages_with_chunking(self):
@@ -507,6 +510,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "year",
@@ -516,6 +520,7 @@ class TestGroupAverage:
 
         xr.testing.assert_allclose(result, expected)
         assert result.ts.attrs == expected.ts.attrs
+        assert result.time.attrs == expected.time.attrs
 
     def test_weighted_seasonal_averages_with_DJF_and_drop_incomplete_seasons(self):
         ds = self.ds.copy()
@@ -555,6 +560,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "season",
@@ -605,6 +611,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "season",
@@ -664,6 +671,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "season",
@@ -716,6 +724,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "season",
@@ -807,6 +816,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "month",
@@ -824,6 +834,7 @@ class TestGroupAverage:
             ),
             coords={"time": self.ds.time, "lat": self.ds.lat, "lon": self.ds.lon},
             dims=["time", "lat", "lon"],
+            attrs={"test_attr": "test"},
         )
 
         result = ds.temporal.group_average("ts", "month")
@@ -831,7 +842,7 @@ class TestGroupAverage:
         expected = expected.drop_dims("time")
         expected["ts"] = xr.DataArray(
             name="ts",
-            data=np.array([[[2.0]], [[0.0]], [[1.0]], [[1.0]], [[2.0]]]),
+            data=np.array([[[2.0]], [[np.nan]], [[1.0]], [[1.0]], [[2.0]]]),
             coords={
                 "lat": expected.lat,
                 "lon": expected.lon,
@@ -856,6 +867,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "month",
@@ -898,6 +910,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "day",
@@ -941,6 +954,7 @@ class TestGroupAverage:
             },
             dims=["time", "lat", "lon"],
             attrs={
+                "test_attr": "test",
                 "operation": "temporal_avg",
                 "mode": "group_average",
                 "freq": "hour",