many new tests for filters

mlair/helpers/filter.py

 import gc
 import warnings
-from typing import Union, Callable, Tuple
+from typing import Union, Callable, Tuple, Dict, Any
 import logging
 import os
 import time
@@ -71,15 +71,13 @@ class ClimateFIRFilter:
         sampling = {1: "1d", 24: "1H"}.get(int(fs))
         logging.debug(f"{plot_name}: create diurnal_anomalies")
         if apriori_diurnal is True and sampling == "1H":
-            diurnal_anomalies = self.create_seasonal_hourly_mean(data, sel_opts=sel_opts, sampling=sampling,
-                                                                 time_dim=time_dim,
+            diurnal_anomalies = self.create_seasonal_hourly_mean(data, time_dim, sel_opts=sel_opts, sampling=sampling,
                                                                  as_anomaly=True)
         else:
             diurnal_anomalies = 0
         logging.debug(f"{plot_name}: create monthly apriori")
         if apriori is None:
-            apriori = self.create_monthly_mean(data, sel_opts=sel_opts, sampling=sampling,
-                                               time_dim=time_dim) + diurnal_anomalies
+            apriori = self.create_monthly_mean(data, time_dim, sel_opts=sel_opts, sampling=sampling) + diurnal_anomalies
             logging.debug(f"{plot_name}: apriori shape = {apriori.shape}")
         apriori_list = to_list(apriori)
         input_data = data.__deepcopy__()
@@ -124,12 +122,9 @@ class ClimateFIRFilter:
             if len(apriori_list) <= i + 1:
                 logging.info(f"{plot_name}: create diurnal_anomalies")
                 if apriori_diurnal is True and sampling == "1H":
-                    # diurnal_anomalies = self.create_hourly_mean(input_data.sel({new_dim: 0}, drop=True),
-                    #                                             sel_opts=sel_opts, sampling=sampling,
-                    #                                             time_dim=time_dim, as_anomaly=True)
                     diurnal_anomalies = self.create_seasonal_hourly_mean(input_data.sel({new_dim: 0}, drop=True),
-                                                                         sel_opts=sel_opts, sampling=sampling,
-                                                                         time_dim=time_dim, as_anomaly=True)
+                                                                         time_dim, sel_opts=sel_opts, sampling=sampling,
+                                                                         as_anomaly=True)
                 else:
                     diurnal_anomalies = 0
                 logging.info(f"{plot_name}: create monthly apriori")
@@ -137,9 +132,8 @@ class ClimateFIRFilter:
                     apriori_list.append(xr.zeros_like(apriori_list[i]) + diurnal_anomalies)
                 elif apriori_type == "residuum_stats":  # calculate monthly statistic on residuum
                     apriori_list.append(
-                        -self.create_monthly_mean(input_data.sel({new_dim: 0}, drop=True), sel_opts=sel_opts,
-                                                  sampling=sampling,
-                                                  time_dim=time_dim) + diurnal_anomalies)
+                        -self.create_monthly_mean(input_data.sel({new_dim: 0}, drop=True), time_dim, sel_opts=sel_opts,
+                                                  sampling=sampling) + diurnal_anomalies)
                 else:
                     raise ValueError(f"Cannot handle unkown apriori type: {apriori_type}. Please choose from None, "
                                      f"`zeros` or `residuum_stats`.")
@@ -179,7 +173,7 @@ class ClimateFIRFilter:
 
         :param data: data to create monthly unity array from, must contain dimension time_dim
         :param time_dim: name of temporal dimension
-        :param extend_range: number of days to extend data
+        :param extend_range: number of days to extend data (default 366)
         :returns: xarray in monthly resolution (centered at 16th day of month) with all values equal to 1
         """
         coords = data.coords
@@ -196,8 +190,8 @@ class ClimateFIRFilter:
         # loffset is required because resampling uses last day in month as resampling timestamp
         return new_array.resample({time_dim: "1m"}, loffset=datetime.timedelta(days=-15)).max()
 
-    def create_monthly_mean(self, data: xr.DataArray, sel_opts: dict = None, sampling: str = "1d",
-                            time_dim: str = "datetime") -> xr.DataArray:
+    def create_monthly_mean(self, data: xr.DataArray, time_dim: str, sel_opts: dict = None, sampling: str = "1d") \
+            -> xr.DataArray:
         """
         Calculate monthly means (12 values) and return a data array with same resolution as given data containing these
         monthly mean values. Sampling points are the 16th of each month (this value is equal to the true monthly mean)
@@ -206,11 +200,11 @@ class ClimateFIRFilter:
         calculate the monthly statistic.
 
         :param data: data to apply statistical calculation on
+        :param time_dim: name of temporal axis
         :param sel_opts: selection options as dict to select a subset of data (default None). A given sel_opts with
             `sel_opts={<time_dim>: "2006"}` forces the method e.g. to derive the monthly means only from data of the
             year 2006.
         :param sampling: sampling of the returned data (default 1d)
-        :param time_dim: name of temporal axis
         :returns: array in desired resolution containing interpolated monthly values. Months with no valid data are
              returned as np.nan which also effects data in the neighbouring months (before / after sampling points which
              are the 16th of each month).
@@ -233,28 +227,67 @@ class ClimateFIRFilter:
         return monthly.resample({time_dim: sampling}).interpolate()
 
     @staticmethod
-    def create_hourly_mean(data, sel_opts=None, sampling="1H", time_dim="datetime", as_anomaly=True):
-        """Calculate hourly statistics. Either the absolute value or the anomaly (as_anomaly=True)."""
-        # can only be used for hourly sampling rate
-        assert sampling == "1H"
-
-        # create unity xarray in hourly resolution
-        hourly = xr.ones_like(data)
-
-        # apply selection if given (only use subset for hourly means)
-        if sel_opts is not None:
-            data = data.sel(**sel_opts)
+    def _compute_hourly_mean_per_month(data: xr.DataArray, time_dim: str, as_anomaly: bool) -> Dict[int, xr.DataArray]:
+        """
+        Calculate for each hour in each month a separate mean value (12 x 24 values in total). Average is either the
+        anomaly of a monthly mean state or the raw mean value.
 
-        # create mean for each hour and replace entries in unity array, calculate anomaly if enabled
-        hourly_mean = data.groupby(f"{time_dim}.hour").mean()
+        :param data: data to calculate averages on
+        :param time_dim: name of temporal dimension
+        :param as_anomaly: indicates whether to calculate means as anomaly of a monthly mean or as raw mean values.
+        :returns: dictionary containing 12 months each with a 24-valued array (1 entry for each hour)
+        """
+        seasonal_hourly_means = {}
+        for month in data.groupby(f"{time_dim}.month").groups.keys():
+            single_month_data = data.sel({time_dim: (data[f"{time_dim}.month"] == month)})
+            hourly_mean = single_month_data.groupby(f"{time_dim}.hour").mean()
             if as_anomaly is True:
                 hourly_mean = hourly_mean - hourly_mean.mean("hour")
-        for hour in hourly_mean.hour.values:
-            loc = (hourly[f"{time_dim}.hour"] == hour)
-            hourly.loc[{f"{time_dim}": loc}] = hourly_mean.sel(hour=hour)
-        return hourly
+            seasonal_hourly_means[month] = hourly_mean
+        return seasonal_hourly_means
+
+    @staticmethod
+    def _create_seasonal_cycle_of_single_hour_mean(result_arr: xr.DataArray, means: Dict[int, xr.DataArray], hour: int,
+                                                   time_dim: str, sampling: str) -> xr.DataArray:
+        """
+        Use monthly means of a given hour to create an array with interpolated values at the indicated hour for each day
+        of the full time span indicated by given result_arr.
+
+        :param result_arr: template array indicating the full time range and additional dimensions to keep
+        :param means: dictionary containing 24 hourly averages for each month (12 x 24 values in total)
+        :param hour: integer of hour of interest
+        :param time_dim: name of temporal dimension
+        :param sampling: sampling rate to interpolate
+        :returns: array with interpolated averages in sampling resolution containing only values for hour of interest
+        """
+        h_coll = xr.ones_like(result_arr) * np.nan
+        for month in means.keys():
+            hourly_mean_single_month = means[month].sel(hour=hour, drop=True)
+            h_coll = xr.where((h_coll[f"{time_dim}.month"] == month), hourly_mean_single_month, h_coll)
+        h_coll = h_coll.resample({time_dim: sampling}).interpolate()
+        h_coll = h_coll.sel({time_dim: (h_coll[f"{time_dim}.hour"] == hour)})
+        return h_coll
+
+    def create_seasonal_hourly_mean(self, data: xr.DataArray, time_dim: str, sel_opts: Dict[str, Any] = None,
+                                    sampling: str = "1H", as_anomaly: bool = True) -> xr.DataArray:
+        """
+        Compute climatological statistics on hourly base either as raw data or anomalies. For each month, an overall
+        mean value (only used if requiring anomalies) and the mean of each hour are calculated. The climatological
+        diurnal cycle is positioned on the 16th of each month and interpolated in between by using a distinct
+        interpolation for each hour of day. The returned array therefore contains data with a yearly cycle (if anomaly
+        is not calculated) or data without a yearly cycle (if using anomalies). In both cases, the data have an
+        amplitude that varies over the year.
+
+        :param data: data to apply this method to
+        :param time_dim: name of temporal axis
+        :param sel_opts: specific selection options that are applied before calculation of climatological statistics
+            (default None)
+        :param sampling: temporal resolution of data (default "1H")
+        :param as_anomaly: specify whether to use anomalies or raw data including a seasonal cycle of the mean value
+            (default: True)
+        :returns: climatological statistics for given data interpolated with given sampling rate
+        """
 
-    def create_seasonal_hourly_mean(self, data, sel_opts=None, sampling="1H", time_dim="datetime", as_anomaly=True):
         """Calculate hourly statistics. Either the absolute value or the anomaly (as_anomaly=True)."""
         # can only be used for hourly sampling rate
         assert sampling == "1H"
@@ -266,29 +299,18 @@ class ClimateFIRFilter:
         # create unity xarray in monthly resolution with sampling point in mid of each month
         monthly = self.create_monthly_unity_array(data, time_dim) * np.nan
 
-        seasonal_hourly_means = {}
-
-        for month in data.groupby(f"{time_dim}.month").groups.keys():
-            # select each month
-            single_month_data = data.sel({time_dim: (data[f"{time_dim}.month"] == month)})
-            hourly_mean = single_month_data.groupby(f"{time_dim}.hour").mean()
-            if as_anomaly is True:
-                hourly_mean = hourly_mean - hourly_mean.mean("hour")
-            seasonal_hourly_means[month] = hourly_mean
+        # calculate for each hour in each month a separate mean value
+        seasonal_hourly_means = self._compute_hourly_mean_per_month(data, time_dim, as_anomaly)
 
+        # create seasonal cycles of these hourly averages
         seasonal_coll = []
         for hour in data.groupby(f"{time_dim}.hour").groups.keys():
-            h_coll = monthly.__deepcopy__()
-            for month in seasonal_hourly_means.keys():
-                hourly_mean_single_month = seasonal_hourly_means[month].sel(hour=hour, drop=True)
-                h_coll = xr.where((h_coll[f"{time_dim}.month"] == month),
-                                  hourly_mean_single_month,
-                                  h_coll)
-            h_coll = h_coll.resample({time_dim: sampling}).interpolate()
-            h_coll = h_coll.sel({time_dim: (h_coll[f"{time_dim}.hour"] == hour)})
-            seasonal_coll.append(h_coll)
-        hourly = xr.concat(seasonal_coll, time_dim).sortby(time_dim).resample({time_dim: sampling}).interpolate()
+            mean_single_hour = self._create_seasonal_cycle_of_single_hour_mean(monthly, seasonal_hourly_means, hour,
+                                                                               time_dim, sampling)
+            seasonal_coll.append(mean_single_hour)
 
+        # combine all cycles in a common data array
+        hourly = xr.concat(seasonal_coll, time_dim).sortby(time_dim).resample({time_dim: sampling}).interpolate()
         return hourly
 
     @staticmethod
@@ -383,7 +405,7 @@ class ClimateFIRFilter:
         return filter_input_data
 
     def create_visualization(self, filtered, data, filter_input_data, plot_dates, time_dim, new_dim, sampling, extend_length_history,
-                             extend_length_future, minimum_length, h, variable_name):
+                             extend_length_future, minimum_length, h, variable_name):  # pragma: no cover
         plot_data = []
         for viz_date in set(plot_dates).intersection(filtered.coords[time_dim].values):
             try:
@@ -456,7 +478,7 @@ class ClimateFIRFilter:
 
         # calculate apriori information from data if not given and extend its range if not sufficient long enough
         if apriori is None:
-            apriori = self.create_monthly_mean(data, sel_opts=sel_opts, sampling=sampling, time_dim=time_dim)
+            apriori = self.create_monthly_mean(data, time_dim, sel_opts=sel_opts, sampling=sampling)
         apriori = apriori.astype(data.dtype)
         apriori = self.extend_apriori(data, apriori, time_dim, sampling, station_name=station_name)
 

test/test_helpers/test_filter.py

@@ -28,14 +28,16 @@ class TestClimateFIRFilter:
     def xr_array(self, data, time_dim):
         start = np.datetime64("2010-01-01 00:00")
         time_index = [start + np.timedelta64(h, "h") for h in range(len(data))]
-        array = xr.DataArray(data, dims=time_dim, coords={time_dim: time_index})
+        array = xr.DataArray(data.reshape(len(data), 1), dims=[time_dim, "station"],
+                             coords={time_dim: time_index, "station": ["DE266X"]})
         return array
 
     @pytest.fixture
     def xr_array_long(self, data, time_dim):
         start = np.datetime64("2010-01-01 00:00")
         time_index = [start + np.timedelta64(175 * h, "h") for h in range(len(data))]
-        array = xr.DataArray(data, dims=time_dim, coords={time_dim: time_index})
+        array = xr.DataArray(data.reshape(len(data), 1), dims=[time_dim, "station"],
+                             coords={time_dim: time_index, "station": ["DE266X"]})
         return array
 
     def test_combine_observation_and_apriori_no_new_dim(self, xr_array, time_dim):
@@ -61,11 +63,12 @@ class TestClimateFIRFilter:
         assert xr.testing.assert_equal(first_entry.sel(window=range(1, 10)), apriori.sel({time_dim: date_pos})) is None
 
     def test_shift_data(self, xr_array, time_dim):
+        remaining_dims = set(xr_array.dims).difference([time_dim])
         obj = object.__new__(ClimateFIRFilter)
         index_values = range(-15, 1)
         res = obj._shift_data(xr_array, index_values, time_dim, new_dim="window")
-        assert len(res.dims) == 2
-        assert len(set(res.dims).difference([time_dim, "window"])) == 0
+        assert len(res.dims) == len(remaining_dims) + 2
+        assert len(set(res.dims).difference([time_dim, "window", *remaining_dims])) == 0
         assert np.testing.assert_array_equal(res.coords["window"].values, np.arange(-15, 1)) is None
         sel = res.sel({time_dim: res.coords[time_dim].values[15]})
         assert sel.sel(window=-15).values == xr_array.sel({time_dim: xr_array.coords[time_dim].values[0]}).values
@@ -135,8 +138,8 @@ class TestClimateFIRFilter:
 
     def test_create_monthly_mean(self, xr_array_long, time_dim):
         obj = object.__new__(ClimateFIRFilter)
-        res = obj.create_monthly_mean(xr_array_long, time_dim=time_dim)
-        assert res.shape == (1462,)
+        res = obj.create_monthly_mean(xr_array_long, time_dim)
+        assert res.shape == (1462, 1)
         assert np.datetime64("2008-12-16") == res.coords[time_dim][0].values
         assert np.datetime64("2012-12-16") == res.coords[time_dim][-1].values
         mean_jan = xr_array_long[xr_array_long[f"{time_dim}.month"] == 1].mean()
@@ -148,10 +151,10 @@ class TestClimateFIRFilter:
 
     def test_create_monthly_mean_sampling(self, xr_array_long, time_dim):
         obj = object.__new__(ClimateFIRFilter)
-        res = obj.create_monthly_mean(xr_array_long, time_dim=time_dim, sampling="1m")
-        assert res.shape == (49,)
-        res = obj.create_monthly_mean(xr_array_long, time_dim=time_dim, sampling="1H")
-        assert res.shape == (35065,)
+        res = obj.create_monthly_mean(xr_array_long, time_dim, sampling="1m")
+        assert res.shape == (49, 1)
+        res = obj.create_monthly_mean(xr_array_long, time_dim, sampling="1H")
+        assert res.shape == (35065, 1)
         mean_jun = xr_array_long[xr_array_long[f"{time_dim}.month"] == 6].mean()
         assert res.sel({time_dim: "2010-06-15T00:00:00"}) == mean_jun
         assert res.sel({time_dim: "2011-06-15T00:00:00"}) == mean_jun
@@ -159,15 +162,57 @@ class TestClimateFIRFilter:
     def test_create_monthly_mean_sel_opts(self, xr_array_long, time_dim):
         obj = object.__new__(ClimateFIRFilter)
         sel_opts = {time_dim: slice("2010-05", "2010-08")}
-        res = obj.create_monthly_mean(xr_array_long, time_dim=time_dim, sel_opts=sel_opts)
+        res = obj.create_monthly_mean(xr_array_long, time_dim, sel_opts=sel_opts)
         assert res.dropna(time_dim)[f"{time_dim}.month"].min() == 5
         assert res.dropna(time_dim)[f"{time_dim}.month"].max() == 8
         mean_jun_2010 = xr_array_long[xr_array_long[f"{time_dim}.month"] == 6].sel({time_dim: "2010"}).mean()
         assert res.sel({time_dim: "2010-06-15T00:00:00"}) == mean_jun_2010
 
-    def test_create_seasonal_hourly_mean(self):
-        #Todo: stopped here
-        pass
+    def test_compute_hourly_mean_per_month(self, xr_array_long, time_dim):
+        obj = object.__new__(ClimateFIRFilter)
+        xr_array_long = xr_array_long.resample({time_dim: "1H"}).interpolate()
+        res = obj._compute_hourly_mean_per_month(xr_array_long, time_dim, True)
+        assert len(res.keys()) == 12
+        assert 6 in res.keys()
+        assert np.testing.assert_almost_equal(res[12].mean(), 0) is None
+        assert np.testing.assert_almost_equal(res[3].mean(), 0) is None
+        assert res[8].shape == (24, 1)
+
+    def test_compute_hourly_mean_per_month_no_anomaly(self, xr_array_long, time_dim):
+        obj = object.__new__(ClimateFIRFilter)
+        xr_array_long = xr_array_long.resample({time_dim: "1H"}).interpolate()
+        res = obj._compute_hourly_mean_per_month(xr_array_long, time_dim, False)
+        assert len(res.keys()) == 12
+        assert 9 in res.keys()
+        assert np.testing.assert_array_less(res[2], res[1]) is None
+
+    def test_create_seasonal_cycle_of_hourly_mean(self, xr_array_long, time_dim):
+        obj = object.__new__(ClimateFIRFilter)
+        xr_array_long = xr_array_long.resample({time_dim: "1H"}).interpolate()
+        monthly = obj.create_monthly_unity_array(xr_array_long, time_dim) * np.nan
+        seasonal_hourly_means = obj._compute_hourly_mean_per_month(xr_array_long, time_dim, True)
+        res = obj._create_seasonal_cycle_of_single_hour_mean(monthly, seasonal_hourly_means, 0, time_dim, "1h")
+        assert res[f"{time_dim}.hour"].sum() == 0
+        assert np.testing.assert_almost_equal(res.sel({time_dim: "2010-12-01"}), res.sel({time_dim: "2011-12-01"})) is None
+        res = obj._create_seasonal_cycle_of_single_hour_mean(monthly, seasonal_hourly_means, 13, time_dim, "1h")
+        assert res[f"{time_dim}.hour"].mean() == 13
+        assert np.testing.assert_almost_equal(res.sel({time_dim: "2010-12-01"}), res.sel({time_dim: "2011-12-01"})) is None
+
+    def test_create_seasonal_hourly_mean(self, xr_array_long, time_dim):
+        obj = object.__new__(ClimateFIRFilter)
+        xr_array_long = xr_array_long.resample({time_dim: "1H"}).interpolate()
+        res = obj.create_seasonal_hourly_mean(xr_array_long, time_dim)
+        assert len(set(res.dims).difference(xr_array_long.dims)) == 0
+        assert res.coords[time_dim][0] < xr_array_long.coords[time_dim][0]
+        assert res.coords[time_dim][-1] > xr_array_long.coords[time_dim][-1]
+
+    def test_create_seasonal_hourly_mean_sel_opts(self, xr_array_long, time_dim):
+        obj = object.__new__(ClimateFIRFilter)
+        xr_array_long = xr_array_long.resample({time_dim: "1H"}).interpolate()
+        sel_opts = {time_dim: slice("2010-05", "2010-08")}
+        res = obj.create_seasonal_hourly_mean(xr_array_long, time_dim, sel_opts=sel_opts)
+        assert res.dropna(time_dim)[f"{time_dim}.month"].min() == 5
+        assert res.dropna(time_dim)[f"{time_dim}.month"].max() == 8
 
     def test_create_unity_array(self, xr_array, time_dim):
         obj = object.__new__(ClimateFIRFilter)
@@ -176,12 +221,12 @@ class TestClimateFIRFilter:
         assert np.datetime64("2011-01-16") == res.coords[time_dim][-1].values
         assert res.max() == res.min()
         assert res.max() == 1
-        assert res.shape == (26,)
+        assert res.shape == (26, 1)
         res = obj.create_monthly_unity_array(xr_array, time_dim, extend_range=0)
-        assert res.shape == (1,)
+        assert res.shape == (1, 1)
         assert np.datetime64("2010-01-16") == res.coords[time_dim][0].values
         res = obj.create_monthly_unity_array(xr_array, time_dim, extend_range=28)
-        assert res.shape == (3,)
+        assert res.shape == (3, 1)
 
     def test_extend_apriori_at_end(self, xr_array_long, time_dim):
         obj = object.__new__(ClimateFIRFilter)
@@ -221,6 +266,15 @@ class TestClimateFIRFilter:
         assert res.stop == np.datetime64("1993-01-01T01:00:00")
         assert res.step is None
 
+    def test_properties(self):
+        obj = object.__new__(ClimateFIRFilter)
+        obj._h = [1, 2, 3]
+        obj._filtered = [4, 5, 63]
+        obj._apriori = [10, 11, 12, 13]
+        assert obj.filter_coefficients == [1, 2, 3]
+        assert obj.filtered_data == [4, 5, 63]
+        assert obj.apriori_data == [10, 11, 12, 13]
+        assert obj.initial_apriori_data == 10
 
 
 class TestFirwinKzf: