first part of skill score calculation (CASE I and II implemented), added copy...

first part of skill score calculation (CASE I and II implemented), added copy of skill_score_plot from felix repository (not checked yet)

requirements.txt

@@ -22,4 +22,5 @@ pyproj
 shapely
 cartopy==0.16.0
 matplotlib
-pillow
\ No newline at end of file
+pillow
+scipy
\ No newline at end of file

src/plotting/postprocessing_plotting.py

@@ -221,3 +221,27 @@ def plot_conditional_quantiles(stations: list, plot_folder: str = ".", rolling_w
     pdf_pages.close()
     plt.close('all')
     logging.info(f"plot_conditional_quantiles() finished after {time}")
+
+
+def plot_climatological_skill_score(d: xr.DataArray, plot_folder=".", score_only=True, extra_nametag="", modelsetup=""):
+    labels = [str(i) + 'd' for i in d.coords['ahead'].values]
+    fig, ax = plt.subplots()
+    if score_only:
+        d = d.loc[:, ['CASE I', 'CASE II', 'CASE III', 'CASE IV'], :]
+        lab_add = ''
+    else:
+        fig.set_size_inches(11.7, 8.27)
+        lab_add = 'terms and '
+    d = d.to_dataframe('data').reset_index(level=[0, 1, 2])
+    sns.boxplot(x='terms', y='data', hue='ahead', data=d, ax=ax, whis=1., palette="Blues_d", showmeans=True,
+                meanprops={'markersize': 1,' markeredgecolor': 'k'}, flierprops={'marker': '.'})
+    ax.axhline(y=0, color='grey', linewidth=.5)
+    ax.set(ylabel=lab_add+'skill score', xlabel='', title='summary of all stations')
+    handles, _ = ax.get_legend_handles_labels()
+    ax.legend(handles, labels)
+    plt.tight_layout()
+    plt.savefig(plot_folder+'SS_Clim_summary_' + extra_nametag + modelsetup + '.pdf', dpi=500)
+    plt.close('all')
+
+    return d
+

src/run_modules/post_processing.py

@@ -10,6 +10,7 @@ import pandas as pd
 import xarray as xr
 import statsmodels.api as sm
 import keras
+from scipy import stats
 
 from src.run_modules.run_environment import RunEnvironment
 from src.data_handling.data_distributor import Distributor
@@ -33,6 +34,7 @@ class PostProcessing(RunEnvironment):
         self.test_data_distributed = Distributor(self.test_data, self.model, self.batch_size)
         self.train_data: DataGenerator = self.data_store.get("generator", "general.train")
         self.plot_path: str = self.data_store.get("plot_path", "general")
+        self.calculate_skill_scores()
         self._run()
 
     def _run(self):
@@ -124,7 +126,7 @@ class PostProcessing(RunEnvironment):
         return nn_prediction_all_stations
 
     @staticmethod
-    def _create_orig_forecast(data, placeholder, mean, std, transformation_method):
+    def _create_orig_forecast(data, _, mean, std, transformation_method):
         return statistics.apply_inverse_transformation(data.label, mean, std, transformation_method)
 
     def _create_ols_forecast(self, input_data, ols_prediction, mean, std, transformation_method):
@@ -195,4 +197,119 @@ class PostProcessing(RunEnvironment):
                 res.loc[match_index, :, k] = v.sel({'datetime': match_index}).squeeze('Stations').transpose()
         return res
 
+    def calculate_skill_scores(self, threshold=60):
+        path = self.data_store.get("forecast_path", "general")
+        for station in self.test_data.stations:  # TODO: replace this by a more general approach to also calculate on train/val
+            file = os.path.join(path, f"forecasts_{station}_test.nc")
+            data = xr.open_dataarray(file)
+            ss = SkillScores()
+            ss.skill_scores(data, station, 3)
+
+
+            # get scaling parameters
+            # mean, std, transformation_method = self.test_data.get_transformation_information(variable='o3')
+        # tmp_nn = statistics.apply_inverse_transformation(tmp_nn, mean, std, transformation_method)
+
+            # self.test_data.get_data_generator(station).restandardise(
+            #     self.get_data_generator(station).data.sel(variables=self.target_var).squeeze('Stations'),
+            #     variables=self.target_var)
+
+
+class SkillScores(RunEnvironment):
+
+    def __init__(self):
+        super().__init__()
+
+    def skill_scores(self, data, station_name, window_lead_time):
+        ahead_names = list(range(1, window_lead_time + 1))
+
+        all_terms_for_clim_deco = ['AI', 'AII', 'AIII', 'AIV', 'BI', 'BII', 'BIV', 'CI', 'CIV', 'CASE I', 'CASE II',
+                                   'CASE III', 'CASE IV']
+        ss_test_clim = xr.DataArray(np.full((len(all_terms_for_clim_deco), len(ahead_names)), np.nan),
+                                    coords=[all_terms_for_clim_deco, ahead_names], dims=['terms', 'ahead'])
+
+        for iahead in ahead_names:
+            ss_test_clim.loc[["CASE I", "AI", "BI", "CI"], iahead] = np.stack(self.skill_score_on_mean_squared_error(
+                data.sel(ahead=iahead), mu_type=1, forecast_name="CNN").values.flatten())
+
+            ss_test_clim.loc[["CASE II", "AII", "BII"], iahead] = np.stack(self.skill_score_on_mean_squared_error(
+                data.sel(ahead=iahead), mu_type=2, forecast_name="CNN").values.flatten())
+
+            # external_climatology = data.sel(variables="orig")
+            #
+            # ss_test_clim.loc[["CASE III", "AIII"], iahead] = np.stack(self.skill_score_on_mean_squared_error(
+            #     data.loc[: iahead, :], mu_type=3, forecast_name="CNN", external_data=external_climatology.mean()
+            # ).values.flatten())
+            #
+            # ss_test_clim.loc[["CASE IV", "AIV", "BIV", "CIV"], iahead] = np.stack(self.skill_score_on_mean_squared_error(
+            #     data.loc[: iahead, :], mu_type=4, forecast_name="CNN", external_data=external_climatology.rename({'datetime': 'index',
+            #                                                                                      'variables': 'type',
+            #                                                                                      'Stations': 'ahead'})).values.flatten())
+
+    def skill_score_on_mean_squared_error(self, data, mu_type=1, observation_name="orig", forecast_name="CNN", external_data=None):
+        kwargs = {"external_data": external_data} if external_data is not None else {}
+        return self.__getattribute__(f"skill_score_mu_case_{mu_type}")(data, observation_name, forecast_name, **kwargs)
+
+    @staticmethod
+    def skill_score_pre_calculations(data, observation_name, forecast_name):
+
+        data = data.loc[..., [observation_name, forecast_name]].drop("ahead")
+        data = data.dropna("index")
+
+        mean = data.mean("index")
+        var = data.var("index")
+        r, p = stats.spearmanr(data.loc[..., [forecast_name, observation_name]])
+
+        AI = np.array(r ** 2)
+        BI = ((r - var.loc[..., forecast_name] / var.loc[..., observation_name]) ** 2).values
+        CI = (((mean.loc[..., forecast_name] - mean.loc[..., observation_name]) / var.loc[
+            ..., observation_name]) ** 2).values
+
+        return AI, BI, CI, data
+
+    def skill_score_mu_case_1(self, data, observation_name="orig", forecast_name="CNN"):
+        AI, BI, CI, data = self.skill_score_pre_calculations(data, observation_name, forecast_name)
+        skill_score = np.array(AI - BI - CI)
+        return pd.DataFrame({"skill_score": [skill_score], "AI": [AI], "BI": [BI], "CI": [CI]}).to_xarray().to_array()
+
+    def skill_score_mu_case_2(self, data, observation_name="orig", forecast_name="CNN"):
+        AI, BI, CI, data = self.skill_score_pre_calculations(data, observation_name, forecast_name)
+        monthly_mean = self.create_monthly_mean_from_daily_data(data)
+        data = xr.concat([data, monthly_mean], dim="type")
+        var = data.var("index")
+        r, p = stats.spearmanr(data.loc[..., [observation_name, observation_name + "X"]])
+        AII = np.array(r ** 2)
+        BII = ((r - var.loc[observation_name + 'X'] / var.loc[observation_name]) ** 2).values
+        skill_score = np.array((AI - BI - CI - AII + BII) / (1 - AII + BII))
+        return pd.DataFrame({"skill_score": [skill_score], "AII": [AII], "BII": [BII]}).to_xarray().to_array()
+
+    def skill_score_mu_case_3(self, data, observation_name="orig", forecast_name="CNN"):
+        AI, BI, CI, data = self.skill_score_pre_calculations(data, observation_name, forecast_name)
+        pass
+
+    def skill_score_mu_case_4(self, data, observation_name="orig", forecast_name="CNN"):
+        AI, BI, CI, data = self.skill_score_pre_calculations(data, observation_name, forecast_name)
+        pass
+
+    @staticmethod
+    def create_monthly_mean_from_daily_data(data, external_data=None, internal_mean=True):
+
+        coordinates = [data.index, [v + "X" for v in list(data.type.values)]]  # TODO
+        empty_data = np.full((len(data.index), len(data.type)), np.nan)
+        monthly_mean = xr.DataArray(empty_data, coords=coordinates, dims=["index", "type"])
+        if internal_mean:
+            mu = data.groupby("index.month").mean()
+        elif not internal_mean and isinstance(external_data, xr.DataArray):
+            mu = external_data.groupby("index.month").mean().drop("type").drop("ahead")
+        else:
+            raise AttributeError(f"Either choose internal_mean=True to calculate the internal mean or use internal_mean"
+                                 f"=False and isinstance(external_data, xarray.DataArray) to get the external mean "
+                                 f"depending on given external data. Given was internal_mean={internal_mean} and "
+                                 f"type(external_data)={type(external_data)} .")
+
+        for month in mu.month:
+            monthly_mean[monthly_mean.index.dt.month == month, :] = mu[mu.month == month].values
+
+        return monthly_mean
+