Repository: metno/pyaerocom (language: python)

Issue:
MDA daily binning misaligned due to mislabeled data in 8hour mean calculation
**Describe the bug**
Please provide a clear and concise description of what the bug is.
 - Pyaerocom version: 0.22.dev0
 - Computing platform: NA
 - Configuration file (if applicable): NA
 - Error message (if applicable): NA

**To Reproduce**
Steps to reproduce the behavior:
1. Add the following parametrization to `test_mda8.py:test_calc_mda8`
```
        pytest.param(
            xr.date_range(start="2024-01-01 06:00:00", periods=30, freq="1h"),
            np.arange(30),
            [np.nan, 25.5]
        ),
```
2. The test fails because the last 8 hour mean included in the first day is for the period 0000-0800 the next day, instead of 1600-2400. Therefore not enough data is included in day two and it results in nan.

**Expected behavior**
MDA8 for 2024-01-02 should be a non-nan numeric value (25.5).

**Screenshots**
If applicable, add screenshots to help explain your problem.

**Additional context**
Caused by wrong alignment. xarray's DataArray.rolling labels data to the left, but labeling it to the right is expected by the current implementation. DataArray.rolling does not support a label-argument so another fix needs to be found. 

Related: https://github.com/metno/pyaerocom/pull/1322


Discussion / hints:

Fix patch:
diff --git a/pyaerocom/stats/mda8/mda8.py b/pyaerocom/stats/mda8/mda8.py
index 54b031c34..99ce4bbe4 100644
--- a/pyaerocom/stats/mda8/mda8.py
+++ b/pyaerocom/stats/mda8/mda8.py
@@ -100,10 +100,11 @@ def _calc_mda8(data: xr.DataArray) -> xr.DataArray:
 
 
 def _rolling_average_8hr(arr: xr.DataArray) -> xr.DataArray:
-    return arr.rolling(time=8, min_periods=6).mean()
+    # Xarray labels the data right based on the last data point in the period for rolling.
+    return arr.rolling(time=8, center=False, min_periods=6).mean()
 
 
 def _daily_max(arr: xr.DataArray) -> xr.DataArray:
-    return arr.resample(time="24H", offset="1h").reduce(
+    return arr.resample(time="24H", origin="start_day", label="left", offset="1h").reduce(
         lambda x, axis: np.apply_along_axis(min_periods_max, 1, x, min_periods=18)
     )


Test patch:
diff --git a/tests/stats/mda8/test_mda8.py b/tests/stats/mda8/test_mda8.py
index 476c5650d..583ac2896 100644
--- a/tests/stats/mda8/test_mda8.py
+++ b/tests/stats/mda8/test_mda8.py
@@ -49,6 +49,13 @@ def test_data(time, values) -> xr.DataArray:
             [np.nan] * 3,
             id="with-nans",
         ),
+        # https://github.com/metno/pyaerocom/issues/1323
+        pytest.param(
+            xr.date_range(start="2024-01-01 06:00:00", periods=30, freq="1h"),
+            np.arange(30),
+            [np.nan, np.nan],
+            id="#1323",
+        ),
     ),
 )
 def test_calc_mda8(test_data, exp_mda8):
@@ -90,16 +97,7 @@ def test_coldata_to_mda8(coldata):
     assert mda8.shape == (2, 8, 1)
 
     assert mda8.data.values[0, :, 0] == pytest.approx(
-        [
-            np.nan,
-            np.nan,
-            1.18741556,
-            1.18777241,
-            1.18869106,
-            1.18879322,
-            1.18807846,
-            1.18700801,
-        ],
+        [np.nan, np.nan, 1.18741556, 1.18777241, 1.18869106, 1.18879322, 1.18807846, 1.18700801],
         abs=10**-5,
         nan_ok=True,
     )
@@ -154,6 +152,26 @@ def test_rollingaverage(test_data, exp_ravg):
     assert all(ravg[0, :, 0] == pytest.approx(exp_ravg, abs=0, nan_ok=True))
 
 
+def test_rollingaverage_label():
+    """
+    Checks that the labels of rolling average is correct (we want it to be labeled based
+    on the "right-most" interval, ie. latest measurement in the interval). This is currently
+    the case in xarray but not greatly documented, so this test checks for that.
+
+    https://github.com/metno/pyaerocom/issues/1323
+    """
+    data = xr.DataArray(
+        [[[x] for x in range(24)]],
+        dims=["data_source", "time", "station_name"],
+        coords={"time": xr.date_range(start="2024-01-01 00:00", periods=24, freq="1h")},
+    )
+
+    ravg = _rolling_average_8hr(data)
+
+    assert ravg["time"].isel(time=0) == np.datetime64("2024-01-01 00:00")
+    assert ravg["time"].isel(time=23) == np.datetime64("2024-01-01 23:00")
+
+
 @pytest.mark.parametrize(
     "time,values,exp_daily_max",
     (