new method to calculate seasonal diurnal anomalies

mlair/helpers/filter.py

@@ -85,7 +85,10 @@ class ClimateFIRFilter:
         sampling = {1: "1d", 24: "1H"}.get(int(fs))
         logging.info(f"{plot_name}: create diurnal_anomalies")
         if apriori_diurnal is True and sampling == "1H":
-            diurnal_anomalies = self.create_hourly_mean(data, sel_opts=sel_opts, sampling=sampling, time_dim=time_dim,
+            # diurnal_anomalies = self.create_hourly_mean(data, sel_opts=sel_opts, sampling=sampling, time_dim=time_dim,
+            #                                             as_anomaly=True)
+            diurnal_anomalies = self.create_seasonal_hourly_mean(data, sel_opts=sel_opts, sampling=sampling,
+                                                                 time_dim=time_dim,
                                                                  as_anomaly=True)
         else:
             diurnal_anomalies = 0
@@ -140,7 +143,10 @@ 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),
+                    # 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)
                 else:
@@ -242,6 +248,43 @@ class ClimateFIRFilter:
             hourly.loc[{f"{time_dim}": loc}] = hourly_mean.sel(hour=hour)
         return hourly
 
+    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"
+
+        # apply selection if given (only use subset for seasonal hourly means)
+        if sel_opts is not None:
+            data = data.sel(**sel_opts)
+
+        # create unity xarray in monthly resolution with sampling point in mid of each month
+        monthly = self.create_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
+
+        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()
+
+        return hourly
+
     @staticmethod
     def extend_apriori(data, apriori, time_dim, sampling="1d"):
         """
@@ -721,7 +764,7 @@ class ClimateFIRFilter:
         rc_params = {'axes.labelsize': 'large',
                      'xtick.labelsize': 'large',
                      'ytick.labelsize': 'large',
-                     'legend.fontsize': 'large',
+                     'legend.fontsize': 'medium',
                      'axes.titlesize': 'large',
                      }
         plt.rcParams.update(rc_params)
@@ -831,7 +874,7 @@ class ClimateFIRFilter:
                             label="clim filter residuum", linewidth=2)
                     ax.set_xlim((ax_start, ax_end))
                     plt.title(f"Residuum of ClimFilter ({str(var)})")
-                    plt.legend()
+                    plt.legend(loc="upper left")
                     fig.autofmt_xdate()
                     plt.tight_layout()
                     plot_name = os.path.join(plot_folder,