From 729e4d399653ed074faa8ecabba4e460fc6009fc Mon Sep 17 00:00:00 2001
From: leufen1 <l.leufen@fz-juelich.de>
Date: Mon, 15 Feb 2021 19:24:11 +0100
Subject: [PATCH] create now forecasts for internal (train, val) and external
(test) data, updated skill score calculation, case IV is still under
construction
---
mlair/configuration/defaults.py | 1 +
mlair/helpers/statistics.py | 45 ++++++++++++-----------
mlair/run_modules/experiment_setup.py | 6 ++-
mlair/run_modules/post_processing.py | 53 +++++++++++++++++----------
4 files changed, 62 insertions(+), 43 deletions(-)
diff --git a/mlair/configuration/defaults.py b/mlair/configuration/defaults.py
index 55a18a1a..04e441fe 100644
--- a/mlair/configuration/defaults.py
+++ b/mlair/configuration/defaults.py
@@ -26,6 +26,7 @@ DEFAULT_EPOCHS = 20
DEFAULT_TARGET_VAR = "o3"
DEFAULT_TARGET_DIM = "variables"
DEFAULT_WINDOW_LEAD_TIME = 3
+DEFAULT_WINDOW_DIM = "window"
DEFAULT_TIME_DIM = "datetime"
DEFAULT_ITER_DIM = "Stations"
DEFAULT_DIMENSIONS = {"new_index": [DEFAULT_TIME_DIM, DEFAULT_ITER_DIM]}
diff --git a/mlair/helpers/statistics.py b/mlair/helpers/statistics.py
index e1887e62..1b1b9284 100644
--- a/mlair/helpers/statistics.py
+++ b/mlair/helpers/statistics.py
@@ -186,21 +186,21 @@ class SkillScores:
Skill score according to Murphy:
Follow climatological skill score definition of Murphy (1988). External data is data from another time period
than the internal data set on initialisation. In other terms, this should be the train and validation data
- whereas the internal data is the test data. This sounds perhaps counter-intuitive, but if a skill score is
- evaluated to a model to another, this must be performend test data set. Therefore, for this case the foreign
- data is train and val data.
+ whereas the external data is the test data. This sounds perhaps counter-intuitive, but if a skill score is
+ evaluated to a model to another, this must be performed on test data set. Therefore, for this case the foreign
+ data is test.
.. code-block:: python
- skill_scores_class = SkillScores(internal_data) # must contain columns obs and CNN.
- skill_scores_clim = skill_scores_class.climatological_skill_scores(external_data, window_lead_time=3)
+ skill_scores_class = SkillScores(external_data) # must contain columns obs and CNN.
+ skill_scores_clim = skill_scores_class.climatological_skill_scores(internal_data, window_lead_time=3)
"""
models_default = ["cnn", "persi", "ols"]
- def __init__(self, internal_data: Data, models=None, observation_name="obs"):
+ def __init__(self, external_data: Data, models=None, observation_name="obs"):
"""Set internal data."""
- self.internal_data = internal_data
+ self.external_data = external_data
self.models = self.set_model_names(models)
self.observation_name = observation_name
@@ -236,7 +236,7 @@ class SkillScores:
combinations, combination_strings = self.get_model_name_combinations()
skill_score = pd.DataFrame(index=combination_strings)
for iahead in ahead_names:
- data = self.internal_data.sel(ahead=iahead)
+ data = self.external_data.sel(ahead=iahead)
skill_score[iahead] = [self.general_skill_score(data,
forecast_name=first,
reference_name=second,
@@ -244,7 +244,7 @@ class SkillScores:
for (first, second) in combinations]
return skill_score
- def climatological_skill_scores(self, external_data: Data, window_lead_time: int,
+ def climatological_skill_scores(self, internal_data: Data, window_lead_time: int,
forecast_name: str) -> xr.DataArray:
"""
Calculate climatological skill scores according to Murphy (1988).
@@ -252,7 +252,7 @@ class SkillScores:
Calculate all CASES I - IV and terms [ABC][I-IV]. Internal data has to be set by initialisation, external data
is part of parameters.
- :param external_data: external data
+ :param internal_data: internal data
:param window_lead_time: interested time step of forecast horizon to select data
:param forecast_name: name of the forecast to use for this calculation (must be available in `data`)
@@ -266,8 +266,8 @@ class SkillScores:
dims=['terms', 'ahead'])
for iahead in ahead_names:
-
- data = self.internal_data.sel(ahead=iahead)
+ data = internal_data.sel(ahead=iahead)
+ external_data = self.external_data.sel(ahead=iahead, type=self.observation_name)
skill_score.loc[["CASE I", "AI", "BI", "CI"], iahead] = np.stack(self._climatological_skill_score(
data, mu_type=1, forecast_name=forecast_name, observation_name=self.observation_name).values.flatten())
@@ -275,14 +275,14 @@ class SkillScores:
skill_score.loc[["CASE II", "AII", "BII"], iahead] = np.stack(self._climatological_skill_score(
data, mu_type=2, forecast_name=forecast_name, observation_name=self.observation_name).values.flatten())
- if external_data is not None:
- skill_score.loc[["CASE III", "AIII"], iahead] = np.stack(self._climatological_skill_score(
- data, mu_type=3, forecast_name=forecast_name, observation_name=self.observation_name,
- external_data=external_data).values.flatten())
+ # if external_data is not None:
+ skill_score.loc[["CASE III", "AIII"], iahead] = np.stack(self._climatological_skill_score(
+ data, mu_type=3, forecast_name=forecast_name, observation_name=self.observation_name,
+ external_data=external_data).values.flatten())
- skill_score.loc[["CASE IV", "AIV", "BIV", "CIV"], iahead] = np.stack(self._climatological_skill_score(
- data, mu_type=4, forecast_name=forecast_name, observation_name=self.observation_name,
- external_data=external_data).values.flatten())
+ skill_score.loc[["CASE IV", "AIV", "BIV", "CIV"], iahead] = np.stack(self._climatological_skill_score(
+ data, mu_type=4, forecast_name=forecast_name, observation_name=self.observation_name,
+ external_data=external_data).values.flatten())
return skill_score
@@ -371,14 +371,15 @@ class SkillScores:
mean, sigma = suffix["mean"], suffix["sigma"]
AIII = (((external_data.mean().values - mean.sel(type=observation_name, drop=True)) / sigma.sel(
type=observation_name, drop=True)) ** 2).values
- skill_score = np.array((AI - BI - CI + AIII) / 1 + AIII)
+ skill_score = np.array((AI - BI - CI + AIII) / (1 + AIII))
return pd.DataFrame({"skill_score": [skill_score], "AIII": [AIII]}).to_xarray().to_array()
def skill_score_mu_case_4(self, data, observation_name, forecast_name, external_data=None):
"""Calculate CASE IV."""
AI, BI, CI, data, suffix = self.skill_score_pre_calculations(data, observation_name, forecast_name)
- monthly_mean_external = self.create_monthly_mean_from_daily_data(external_data, columns=data.type.values,
- index=data.index)
+ monthly_mean_external = self.create_monthly_mean_from_daily_data(external_data, columns=data.type.values)
+ # monthly_mean_external = self.create_monthly_mean_from_daily_data(external_data, columns=data.type.values,
+ # index=data.index)
data = xr.concat([data, monthly_mean_external], dim="type")
mean, sigma = suffix["mean"], suffix["sigma"]
monthly_mean_external = self.create_monthly_mean_from_daily_data(external_data, columns=data.type.values)
diff --git a/mlair/run_modules/experiment_setup.py b/mlair/run_modules/experiment_setup.py
index 18ef98f8..af540fc2 100644
--- a/mlair/run_modules/experiment_setup.py
+++ b/mlair/run_modules/experiment_setup.py
@@ -13,7 +13,7 @@ from mlair.configuration.defaults import DEFAULT_STATIONS, DEFAULT_VAR_ALL_DICT,
DEFAULT_HPC_LOGIN_LIST, DEFAULT_HPC_HOST_LIST, DEFAULT_CREATE_NEW_MODEL, DEFAULT_TRAIN_MODEL, \
DEFAULT_FRACTION_OF_TRAINING, DEFAULT_EXTREME_VALUES, DEFAULT_EXTREMES_ON_RIGHT_TAIL_ONLY, DEFAULT_PERMUTE_DATA, \
DEFAULT_BATCH_SIZE, DEFAULT_EPOCHS, DEFAULT_TARGET_VAR, DEFAULT_TARGET_DIM, DEFAULT_WINDOW_LEAD_TIME, \
- DEFAULT_DIMENSIONS, DEFAULT_TIME_DIM, DEFAULT_INTERPOLATION_METHOD, DEFAULT_INTERPOLATION_LIMIT, \
+ DEFAULT_WINDOW_DIM, DEFAULT_DIMENSIONS, DEFAULT_TIME_DIM, DEFAULT_INTERPOLATION_METHOD, DEFAULT_INTERPOLATION_LIMIT, \
DEFAULT_TRAIN_START, DEFAULT_TRAIN_END, DEFAULT_TRAIN_MIN_LENGTH, DEFAULT_VAL_START, DEFAULT_VAL_END, \
DEFAULT_VAL_MIN_LENGTH, DEFAULT_TEST_START, DEFAULT_TEST_END, DEFAULT_TEST_MIN_LENGTH, DEFAULT_TRAIN_VAL_MIN_LENGTH, \
DEFAULT_USE_ALL_STATIONS_ON_ALL_DATA_SETS, DEFAULT_EVALUATE_BOOTSTRAPS, DEFAULT_CREATE_NEW_BOOTSTRAPS, \
@@ -211,6 +211,7 @@ class ExperimentSetup(RunEnvironment):
target_var="o3",
target_dim=None,
window_lead_time: int = None,
+ window_dim=None,
dimensions=None,
time_dim=None,
iter_dim=None,
@@ -299,8 +300,9 @@ class ExperimentSetup(RunEnvironment):
self._set_param("transformation", None, scope="preprocessing")
self._set_param("data_handler", data_handler, default=DefaultDataHandler)
- # iteration
+ # iter and window dimension
self._set_param("iter_dim", iter_dim, default=DEFAULT_ITER_DIM)
+ self._set_param("window_dim", window_dim, default=DEFAULT_WINDOW_DIM)
# target
self._set_param("target_var", target_var, default=DEFAULT_TARGET_VAR)
diff --git a/mlair/run_modules/post_processing.py b/mlair/run_modules/post_processing.py
index d223858c..c810d3c5 100644
--- a/mlair/run_modules/post_processing.py
+++ b/mlair/run_modules/post_processing.py
@@ -89,8 +89,9 @@ class PostProcessing(RunEnvironment):
# ols model
self.train_ols_model()
- # forecasts
- self.make_prediction()
+ # forecasts on test data
+ self.make_prediction(self.test_data)
+ self.make_prediction(self.train_val_data)
# calculate error metrics on test data
self.calculate_test_score()
@@ -397,7 +398,7 @@ class PostProcessing(RunEnvironment):
"""Train ordinary least squared model on train data."""
self.ols_model = OrdinaryLeastSquaredModel(self.train_data)
- def make_prediction(self):
+ def make_prediction(self, subset):
"""
Create predictions for NN, OLS, and persistence and add true observation as reference.
@@ -405,10 +406,12 @@ class PostProcessing(RunEnvironment):
predictions for a single station are stored locally under `<forecast/forecast_norm>_<station>_test.nc` and can
be found inside `forecast_path`.
"""
- logging.debug("start make_prediction")
+ subset_type = subset.name
+ logging.debug(f"start make_prediction for {subset_type}")
time_dimension = self.data_store.get("time_dim")
window_dim = self.data_store.get("window_dim")
- for i, data in enumerate(self.test_data):
+ subset_type = subset.name
+ for i, data in enumerate(subset):
input_data = data.get_X()
target_data = data.get_Y(as_numpy=False)
observation_data = data.get_observation()
@@ -446,7 +449,7 @@ class PostProcessing(RunEnvironment):
# save all forecasts locally
path = self.data_store.get("forecast_path")
prefix = "forecasts_norm" if normalised is True else "forecasts"
- file = os.path.join(path, f"{prefix}_{str(data)}_test.nc")
+ file = os.path.join(path, f"{prefix}_{str(data)}_{subset_type}.nc")
all_predictions.to_netcdf(file)
def _get_frequency(self) -> str:
@@ -612,21 +615,33 @@ class PostProcessing(RunEnvironment):
res.loc[match_index, :, k] = v.loc[match_index]
return res
- def _get_external_data(self, station: str) -> Union[xr.DataArray, None]:
+ def _get_internal_data(self, station: str, path: str) -> Union[xr.DataArray, None]:
"""
- Get external data for given station.
+ Get internal data for given station.
- External data is defined as data that is not part of the observed period. From an evaluation perspective, this
+ Internal data is defined as data that is already known to the model. From an evaluation perspective, this
refers to data, that is no test data, and therefore to train and val data.
+ :param station: name of station to load internal data.
+ """
+ try:
+ file = os.path.join(path, f"forecasts_{str(station)}_train_val.nc")
+ return xr.open_dataarray(file)
+ except (IndexError, KeyError):
+ return None
+
+ def _get_external_data(self, station: str, path: str) -> Union[xr.DataArray, None]:
+ """
+ Get external data for given station.
+
+ External data is defined as data that is not known to the model. From an evaluation perspective, this refers to
+ data, that is not train or val data, and therefore to test data.
+
:param station: name of station to load external data.
"""
try:
- data = self.train_val_data[station]
- observation = data.get_observation()
- transformation_func = data.apply_transformation
- external_data = self._create_observation(observation, None, transformation_func, normalised=False)
- return external_data.rename({external_data.dims[0]: 'index'})
+ file = os.path.join(path, f"forecasts_{str(station)}_test.nc")
+ return xr.open_dataarray(file)
except (IndexError, KeyError):
return None
@@ -644,13 +659,13 @@ class PostProcessing(RunEnvironment):
skill_score_competitive = {}
skill_score_climatological = {}
for station in self.test_data:
- file = os.path.join(path, f"forecasts_{str(station)}_test.nc")
- data = xr.open_dataarray(file)
+ external_data = self._get_external_data(station, path)
competitor = self.load_competitors(str(station))
- combined = xr.concat([data, competitor], dim="type") if competitor is not None else data
+ combined = xr.concat([external_data, competitor], dim="type") if competitor is not None else external_data
skill_score = statistics.SkillScores(combined, models=remove_items(list(combined.type.values), "obs"))
- external_data = self._get_external_data(station) # ToDo: check if external is still right?
skill_score_competitive[station] = skill_score.skill_scores(self.window_lead_time)
+
+ internal_data = self._get_internal_data(station, path) # ToDo: check if external is still right?
skill_score_climatological[station] = skill_score.climatological_skill_scores(
- external_data, self.window_lead_time, forecast_name=self.forecast_indicator)
+ internal_data, self.window_lead_time, forecast_name=self.forecast_indicator)
return skill_score_competitive, skill_score_climatological
--
GitLab