diff --git a/requirements.txt b/requirements.txt
index 270c084865fbff00e6346b5f267c8d939a1d9902..9ccd09ef9234bafff4559aa9fc325bef7d8bf3ea 100644
--- a/requirements.txt
+++ b/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
diff --git a/src/plotting/postprocessing_plotting.py b/src/plotting/postprocessing_plotting.py
index eb3f7f8c058fae47c2703e6e53ee22bdc013f7a2..09402f0425748a09f637734455c26fdfd25ce071 100644
--- a/src/plotting/postprocessing_plotting.py
+++ b/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
+
diff --git a/src/run_modules/post_processing.py b/src/run_modules/post_processing.py
index d773c9789619c802f05f3d922d554fb1b246ffe0..ce6d719f9d51bad00b2aea9f75b4f4dcdd722fa7 100644
--- a/src/run_modules/post_processing.py
+++ b/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
+