Climatological op

HISTORY.rst

@@ -4,14 +4,15 @@ History
 
 v0.8.0 (unreleased)
 -------------------
-Contributors to this version: Gabriel Rondeau-Genesse (:user:`RondeauG`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), Sarah-Claude Bourdeau-Goulet (:user:`sarahclaude`).
+Contributors to this version: Gabriel Rondeau-Genesse (:user:`RondeauG`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), Sarah-Claude Bourdeau-Goulet (:user:`sarahclaude`), Marco Braun (:user:`vindelico`).
 
 Announcements
 ^^^^^^^^^^^^^
 * N/A
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* Replaced aggregation function ``climatological_mean()`` with ``climatological_op()`` offering more types of operations to aggregate over climate periods. (:pull:`290`)
 * Added the ability to search for simulations that reach a given warming level. (:pull:`251`).
 * ``xs.spatial_mean`` now accepts the ``region="global"`` keyword to perform a global average (:issue:`94`, :pull:`260`).
 * ``xs.spatial_mean`` with ``method='xESMF'`` will also automatically segmentize polygons (down to a 1° resolution) to ensure a correct average (:pull:`260`).

tests/test_aggregate.py

@@ -10,6 +10,12 @@
 
 
 class TestClimatologicalMean:
+    def _format(self, s):
+        import xclim
+        op_format = (dict.fromkeys(("mean", "std", "var", "sum"), "adj") |
+                     dict.fromkeys(("max", "min"), "noun"))
+        return xclim.core.formatting.default_formatter.format_field(s, op_format[s])
+
     def test_daily(self):
         ds = timeseries(
             np.tile(np.arange(1, 13), 3),
@@ -42,10 +48,12 @@ def test_all_default(self, xrfreq):
         # Test metadata
         assert (
             out.tas.attrs["description"]
-            == f"30-year mean of {ds.tas.attrs['description']}"
+            # == f"30-year mean of {ds.tas.attrs['description']}" # Changed to reflect output from climatological_op
+            == f"Climatological 30-year average of {ds.tas.attrs['description']}"
         )
         assert (
-            "30-year rolling average (non-centered) with a minimum of 30 years of data"
+            # "30-year rolling average (non-centered) with a minimum of 30 years of data"
+            "Climatological 30-year average (non-centered) with a minimum of 30 years of data"
             in out.tas.attrs["history"]
         )
         assert out.attrs["cat:processing_level"] == "climatology"
@@ -76,10 +84,12 @@ def test_options(self, xrfreq):
         # Test metadata
         assert (
             out.tas.attrs["description"]
-            == f"15-year mean of {ds.tas.attrs['description']}"
+            # == f"15-year mean of {ds.tas.attrs['description']}" # Changed to reflect output from climatological_op
+            == f"Climatological 15-year average of {ds.tas.attrs['description']}"
         )
         assert (
-            "15-year rolling average (non-centered) with a minimum of 15 years of data"
+            # "15-year rolling average (non-centered) with a minimum of 15 years of data"
+            "Climatological 15-year average (non-centered) with a minimum of 15 years of data"
             in out.tas.attrs["history"]
         )
         assert out.attrs["cat:processing_level"] == "for_testing"
@@ -356,12 +366,15 @@ def test_options(self, periods, to_level):
         assert out.attrs["cat:xrfreq"] == "fx"
         assert all(v in out for v in ["tg_min", "growing_degree_days"])
         assert (
-            f"{30 if periods is None else int(periods[0][1]) - int(periods[0][0]) + 1}-year mean of"
+            # f"{30 if periods is None else int(periods[0][1]) - int(periods[0][0]) + 1}-year mean of"
+            f"Climatological {30 if periods is None else int(periods[0][1]) - int(periods[0][0]) + 1}-year average of"
             in out.tg_min.attrs["description"]
         )
         assert (
             out.tg_min.attrs["description"].split(
-                f"{30 if periods is None else int(periods[0][1]) - int(periods[0][0]) + 1}-year mean of "
+                # f"{30 if periods is None else int(periods[0][1]) - int(periods[0][0]) + 1}-year mean of "
+                f"Climatological {30 if periods is None else int(periods[0][1]) - int(periods[0][0]) + 1}"
+                f"-year average of "
             )[1]
             != self.ds.tas.attrs["description"]
         )
@@ -548,3 +561,196 @@ def test_global(self, lonstart, method, exp):
 
         avg = xs.aggregate.spatial_mean(ds, method=method, region="global")
         np.testing.assert_allclose(avg.tas, exp)
+
+
+class TestClimatologicalOp:
+    def _format(self, s):
+        import xclim
+        op_format = (dict.fromkeys(("mean", "std", "var", "sum"), "adj") |
+                     dict.fromkeys(("max", "min"), "noun"))
+        return xclim.core.formatting.default_formatter.format_field(s, op_format[s])
+
+    def test_daily(self):
+        ds = timeseries(
+            np.tile(np.arange(1, 13), 3),
+            variable="tas",
+            start="2001-01-01",
+            freq="D",
+            as_dataset=True,
+        )
+        with pytest.raises(NotImplementedError):
+            xs.climatological_op(ds, op='mean')
+
+    @pytest.mark.parametrize("xrfreq", ["MS", "AS-JAN"])
+    @pytest.mark.parametrize("op", ["max", "mean", "median", "min", "std", "sum", "var", "linregress"])
+    def test_all_default(self, xrfreq, op):
+        o = 12 if xrfreq == "MS" else 1
+
+        ds = timeseries(
+            np.tile(np.arange(1, o + 1), 30),
+            variable="tas",
+            start="2001-01-01",
+            freq=xrfreq,
+            as_dataset=True,
+        )
+        out = xs.climatological_op(ds, op=op)
+        expected = (dict.fromkeys(("max", "mean", "median", "min"), np.arange(1, o + 1)) |
+                    dict.fromkeys(("std", "var"), np.zeros(o)) |
+                    dict({"sum": np.arange(1, o + 1) * 30}) |
+                    dict({"linregress": np.array([
+                        np.zeros(o),
+                        np.arange(1, o + 1),
+                        np.zeros(o),
+                        np.ones(o),
+                        np.zeros(o),
+                        np.zeros(o)]).T})
+                    )
+        # Test output variable name, values, length, horizon
+        assert list(out.data_vars.keys()) == [f"tas_clim_{op}"]
+        np.testing.assert_array_equal(out[f"tas_clim_{op}"], expected[op])
+        assert len(out.time) == (o * len(np.unique(out.horizon.values)))
+        np.testing.assert_array_equal(out.time[0], ds.time[0])
+        assert (out.horizon == "2001-2030").all()
+        # Test metadata
+        operation = self._format(op) if op not in ['median', 'linregress'] else op
+        assert (
+            out[f"tas_clim_{op}"].attrs["description"]
+            == f"Climatological 30-year {operation} of {ds.tas.attrs['description']}"
+        )
+        assert (
+            f"Climatological 30-year {operation} (non-centered) with a minimum of 30 years of data"
+            in out[f"tas_clim_{op}"].attrs["history"]
+        )
+        assert out.attrs["cat:processing_level"] == "climatology"
+
+    @pytest.mark.parametrize("xrfreq", ["MS", "AS-JAN"])
+    @pytest.mark.parametrize("op", ['max', 'mean', 'median', 'min', 'std', 'sum', 'var', 'linregress'])
+    def test_options(self, xrfreq, op):
+        o = 12 if xrfreq == "MS" else 1
+
+        ds = timeseries(
+            np.tile(np.arange(1, o + 1), 30),
+            variable="tas",
+            start="2001-01-01",
+            freq=xrfreq,
+            as_dataset=True,
+        )
+        out = xs.climatological_op(ds, op=op, window=15, stride=5, to_level="for_testing")
+        expected = (dict.fromkeys(("max", "mean", "median", "min"),
+                                  np.tile(np.arange(1, o + 1), len(np.unique(out.horizon.values)))) |
+                    dict.fromkeys(("std", "var"),
+                                  np.tile(np.zeros(o), len(np.unique(out.horizon.values)))) |
+                    dict({"sum": np.tile(np.arange(1, o + 1) * 15, len(np.unique(out.horizon.values)))}) |
+                    dict({"linregress": np.tile(np.array([
+                        np.zeros(o),
+                        np.arange(1, o + 1),
+                        np.zeros(o),
+                        np.ones(o),
+                        np.zeros(o),
+                        np.zeros(o)]), len(np.unique(out.horizon.values))).T})
+                    )
+        # Test output values
+        np.testing.assert_array_equal(
+            out[f"tas_clim_{op}"],
+            expected[op],
+        )
+        assert len(out.time) == (o * len(np.unique(out.horizon.values)))
+        np.testing.assert_array_equal(out.time[0], ds.time[0])
+        assert {"2001-2015", "2006-2020", "2011-2025", "2016-2030"}.issubset(
+            out.horizon.values
+        )
+        # Test metadata
+        operation = self._format(op) if op not in ['median', 'linregress'] else op
+        assert (
+            out[f"tas_clim_{op}"].attrs["description"]
+            == f"Climatological 15-year {operation} of {ds.tas.attrs['description']}"
+        )
+        assert (
+            f"Climatological 15-year {operation} (non-centered) with a minimum of 15 years of data"
+            in out[f"tas_clim_{op}"].attrs["history"]
+        )
+        assert out.attrs["cat:processing_level"] == "for_testing"
+
+    @pytest.mark.parametrize("op", ["mean", "linregress"])
+    def test_minperiods(self, op):
+        ds = timeseries(
+            np.tile(np.arange(1, 5), 30),
+            variable="tas",
+            start="2001-03-01",
+            freq="QS-DEC",
+            as_dataset=True,
+        )
+        ds = ds.where(ds["time"].dt.strftime("%Y-%m-%d") != "2030-12-01")
+
+        op = 'mean'
+        out = xs.climatological_op(ds, op=op, window=30)
+        assert all(np.isreal(out[f"tas_clim_{op}"]))
+        assert len(out.time) == 4
+        np.testing.assert_array_equal(out[f"tas_clim_{op}"], np.arange(1, 5))
+
+        # min_periods as int
+        out = xs.climatological_op(ds, op=op, window=30, min_periods=30)
+        assert np.sum(np.isnan(out[f"tas_clim_{op}"])) == 1
+
+        # min_periods as float
+        out = xs.climatological_op(ds, op=op, window=30, min_periods=.5)
+        assert "minimum of 15 years of data" in out[f"tas_clim_{op}"].attrs["history"]
+        assert np.sum(np.isnan(out[f"tas_clim_{op}"])) == 0
+
+        with pytest.raises(ValueError):
+            xs.climatological_op(ds, op=op, window=5, min_periods=6)
+
+    @pytest.mark.parametrize("op", ["mean", "linregress"])
+    def test_periods(self, op):
+        ds1 = timeseries(
+            np.tile(np.arange(1, 2), 10),
+            variable="tas",
+            start="2001-01-01",
+            freq="AS-JAN",
+            as_dataset=True,
+        )
+        ds2 = timeseries(
+            np.tile(np.arange(1, 2), 10),
+            variable="tas",
+            start="2021-01-01",
+            freq="AS-JAN",
+            as_dataset=True,
+        )
+
+        ds = xr.concat([ds1, ds2], dim="time")
+        with pytest.raises(ValueError):
+            xs.climatological_op(ds, op=op)
+
+        out = xs.climatological_op(ds, op='mean', periods=[["2001", "2010"], ["2021", "2030"]])
+        assert len(out.time) == 2
+        assert {"2001-2010", "2021-2030"}.issubset(out.horizon.values)
+
+    @pytest.mark.parametrize("cal", ["proleptic_gregorian", "noleap", "360_day"])
+    def test_calendars(self, cal):
+        ds = timeseries(
+            np.tile(np.arange(1, 2), 30),
+            variable="tas",
+            start="2001-01-01",
+            freq="AS-JAN",
+            as_dataset=True,
+        )
+
+        out = xs.climatological_op(ds.convert_calendar(cal, align_on="date"), op='mean')
+        assert out.time.dt.calendar == cal
+
+    def test_periods_as_dim(self):
+        ds = timeseries(
+            np.tile(np.arange(1, 13), 30),
+            variable="tas",
+            start="2001-01-01",
+            freq="MS",
+            as_dataset=True,
+        )
+        out = xs.climatological_op(ds, op="mean", window=10, stride=5, periods_as_dim=True)
+        assert (out.tas_clim_mean.values == np.tile(np.arange(1, 13), (5, 1))).all()
+        assert out.dims == {'period': 5, 'month': 12}
+        assert out.time.dims == ('period', 'month')
+        assert (out.period.values == ['2001-2010', '2006-2015', '2011-2020', '2016-2025', '2021-2030']).all()
+        assert (out.month.values ==
+                ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
+                 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']).all()