diff --git a/src/plotting/postprocessing_plotting.py b/src/plotting/postprocessing_plotting.py
index 32a84f25fefb24f40a40681592f1c9a7c5c5afcf..19c511d9b6f8442ea30ec7f7623b0f2f598f5cc2 100644
--- a/src/plotting/postprocessing_plotting.py
+++ b/src/plotting/postprocessing_plotting.py
@@ -191,10 +191,16 @@ class PlotStationMap(AbstractPlotClass):
@TimeTrackingWrapper
class PlotConditionalQuantiles(AbstractPlotClass):
"""
- This class creates conditional quantile plots as originally proposed by Murphy, Brown and Chen (1989)
+ This class creates cond.quantile plots as originally proposed by Murphy, Brown and Chen (1989) [But in log scale]
Link to paper: https://journals.ametsoc.org/doi/pdf/10.1175/1520-0434%281989%29004%3C0485%3ADVOTF%3E2.0.CO%3B2
"""
+ # ignore warnings if nans appear in quantile grouping
+ warnings.filterwarnings("ignore", message="All-NaN slice encountered")
+ # ignore warnings if mean is calculated on nans
+ warnings.filterwarnings("ignore", message="Mean of empty slice")
+ # ignore warnings for y tick = 0 on log scale (instead of 0.00001 or similar)
+ warnings.filterwarnings("ignore", message="Attempted to set non-positive bottom ylim on a log-scaled axis.")
def __init__(self, stations: List, data_pred_path: str, plot_folder: str = ".", plot_per_seasons=True,
rolling_window: int = 3, model_mame: str = "CNN", obs_name: str = "obs", **kwargs):
@@ -210,6 +216,8 @@ class PlotConditionalQuantiles(AbstractPlotClass):
:param kwargs: Some further arguments which are listed in self._opts
"""
+ logging.debug(f"started plot_conditional_quantiles()")
+
super().__init__(plot_folder, "conditional_quantiles")
self._data_pred_path = data_pred_path
@@ -225,7 +233,6 @@ class PlotConditionalQuantiles(AbstractPlotClass):
'reference 1:1']),
"data_unit": kwargs.get("data_unit", "ppb"),
}
- # self._data_unit = kwargs.get("data_unit", "ppb")
if plot_per_seasons is True:
self.seasons = ['DJF', 'MAM', 'JJA', 'SON']
else:
@@ -428,142 +435,6 @@ class PlotConditionalQuantiles(AbstractPlotClass):
plt.close('all')
-@TimeTrackingWrapper
-def plot_conditional_quantiles(stations: list, plot_folder: str = ".", rolling_window: int = 3, ref_name: str = 'obs',
- pred_name: str = 'CNN', season: str = "", forecast_path: str = None,
- plot_name_affix: str = "", units: str = "ppb"):
- """
- This plot was originally taken from Murphy, Brown and Chen (1989):
- https://journals.ametsoc.org/doi/pdf/10.1175/1520-0434%281989%29004%3C0485%3ADVOTF%3E2.0.CO%3B2
-
- :param stations: stations to include in the plot (forecast data needs to be available already)
- :param plot_folder: path to save the plot (default: current directory)
- :param rolling_window: the rolling window mean will smooth the plot appearance (no smoothing in bin calculation,
- this is only a cosmetic step, default: 3)
- :param ref_name: name of the reference data series
- :param pred_name: name of the investigated data series
- :param season: season name to highlight if not empty
- :param forecast_path: path to save the plot file
- :param plot_name_affix: name to specify this plot (e.g. 'cali-ref', default: '')
- :param units: units of the forecasted values (default: ppb)
- """
- # time = TimeTracking()
- logging.debug(f"started plot_conditional_quantiles()")
- # ignore warnings if nans appear in quantile grouping
- warnings.filterwarnings("ignore", message="All-NaN slice encountered")
- # ignore warnings if mean is calculated on nans
- warnings.filterwarnings("ignore", message="Mean of empty slice")
- # ignore warnings for y tick = 0 on log scale (instead of 0.00001 or similar)
- warnings.filterwarnings("ignore", message="Attempted to set non-positive bottom ylim on a log-scaled axis.")
-
- def load_data():
- logging.debug("... load data")
- data_collector = []
- for station in stations:
- file = os.path.join(forecast_path, f"forecasts_{station}_test.nc")
- data_tmp = xr.open_dataarray(file)
- data_collector.append(data_tmp.loc[:, :, ['CNN', 'obs', 'OLS']].assign_coords(station=station))
- return xr.concat(data_collector, dim='station').transpose('index', 'type', 'ahead', 'station')
-
- def segment_data(data):
- logging.debug("... segment data")
- # combine index and station to multi index
- data = data.stack(z=['index', 'station'])
- # replace multi index by simple position index (order is not relevant anymore)
- data.coords['z'] = range(len(data.coords['z']))
- # segment data of pred_name into bins
- data.loc[pred_name, ...] = data.loc[pred_name, ...].to_pandas().T.apply(pd.cut, bins=bins,
- labels=bins[1:]).T.values
- return data
-
- def create_quantile_panel(data, q):
- logging.debug("... create quantile panel")
- # create empty xarray with dims: time steps ahead, quantiles, bin index (numbers create in previous step)
- quantile_panel = xr.DataArray(np.full([data.ahead.shape[0], len(q), bins[1:].shape[0]], np.nan),
- coords=[data.ahead, q, bins[1:]], dims=['ahead', 'quantiles', 'categories'])
- # ensure that the coordinates are in the right order
- quantile_panel = quantile_panel.transpose('ahead', 'quantiles', 'categories')
- # calculate for each bin of the pred_name data the quantiles of the ref_name data
- for bin in bins[1:]:
- mask = (data.loc[pred_name, ...] == bin)
- quantile_panel.loc[..., bin] = data.loc[ref_name, ...].where(mask).quantile(q, dim=['z']).T
-
- return quantile_panel
-
- def labels(plot_type, data_unit="ppb"):
- names = (f"forecast concentration (in {data_unit})", f"observed concentration (in {data_unit})")
- if plot_type == "obs":
- return names
- else:
- return names[::-1]
-
- xlabel, ylabel = labels(ref_name, units)
-
- opts = {"q": [.1, .25, .5, .75, .9], "linetype": [':', '-.', '--', '-.', ':'],
- "legend": ['.10th and .90th quantile', '.25th and .75th quantile', '.50th quantile', 'reference 1:1'],
- "xlabel": xlabel, "ylabel": ylabel}
-
- # set name and path of the plot
- base_name = "conditional_quantiles-orig"
- def add_affix(x): return f"_{x}" if len(x) > 0 else ""
- plot_name = f"{base_name}{add_affix(season)}{add_affix(plot_name_affix)}_plot.pdf"
- plot_path = os.path.join(os.path.abspath(plot_folder), plot_name)
-
- # check forecast path
- if forecast_path is None:
- raise ValueError("Forecast path is not given but required.")
-
- # load data and set data bins
- orig_data = load_data()
- bins = np.arange(0, math.ceil(orig_data.max().max()) + 1, 1).astype(int)
- segmented_data = segment_data(orig_data)
- quantile_panel = create_quantile_panel(segmented_data, q=opts["q"])
-
- # init pdf output
- pdf_pages = matplotlib.backends.backend_pdf.PdfPages(plot_path)
- logging.debug(f"... plot path is {plot_path}")
-
- # create plot for each time step ahead
- y2_max = 0
- for iteration, d in enumerate(segmented_data.ahead):
- logging.debug(f"... plotting {d.values} time step(s) ahead")
- # plot smoothed lines with rolling mean
- smooth_data = quantile_panel.loc[d, ...].rolling(categories=rolling_window, center=True).mean().to_pandas().T
- ax = smooth_data.plot(style=opts["linetype"], color='black', legend=False)
- ax2 = ax.twinx()
- # add reference line
- ax.plot([0, bins.max()], [0, bins.max()], color='k', label='reference 1:1', linewidth=.8)
- # add histogram of the segmented data (pred_name)
- handles, labels = ax.get_legend_handles_labels()
- segmented_data.loc[pred_name, d, :].to_pandas().hist(bins=bins, ax=ax2, color='k', alpha=.3, grid=False,
- rwidth=1)
- # add legend
- plt.legend(handles[:3] + [handles[-1]], opts["legend"], loc='upper left', fontsize='large')
- # adjust limits and set labels
- ax.set(xlim=(0, bins.max()), ylim=(0, bins.max()))
- ax.set_xlabel(opts["xlabel"], fontsize='x-large')
- ax.tick_params(axis='x', which='major', labelsize=15)
- ax.set_ylabel(opts["ylabel"], fontsize='x-large')
- ax.tick_params(axis='y', which='major', labelsize=15)
- ax2.yaxis.label.set_color('gray')
- ax2.tick_params(axis='y', colors='gray')
- ax2.yaxis.labelpad = -15
- ax2.set_yscale('log')
- if iteration == 0:
- y2_max = ax2.get_ylim()[1] + 100
- ax2.set(ylim=(0, y2_max * 10 ** 8), yticks=np.logspace(0, 4, 5))
- ax2.set_ylabel(' sample size', fontsize='x-large')
- ax2.tick_params(axis='y', which='major', labelsize=15)
- # set title and save current figure
- title = f"{d.values} time step(s) ahead{f' ({season})' if len(season) > 0 else ''}"
- plt.title(title)
- pdf_pages.savefig()
- # close all open figures / plots
- pdf_pages.close()
- plt.close('all')
- #logging.info(f"plot_conditional_quantiles() finished after {time}")
-
-
@TimeTrackingWrapper
class PlotClimatologicalSkillScore(AbstractPlotClass):
"""
@@ -970,8 +841,3 @@ if __name__ == "__main__":
con_quan_cls = PlotConditionalQuantiles(stations, path, plt_path)
- # con_quan_cls = PlotConditionalQuantiles(stations, data_pred_path=path, plot_name_affix="", pred_name="CNN",
- # ref_name="obs", plot_folder=plt_path, seasons=None)
- plot_conditional_quantiles(stations, pred_name="CNN", ref_name="obs",
- forecast_path=path, plot_name_affix="cali-ref-orig", plot_folder=plt_path)
-
diff --git a/src/run_modules/post_processing.py b/src/run_modules/post_processing.py
index b1f3afeae4e4c15ffdbbe0da252737256ebe5687..e8353048f2fa982a9d54900fab93142ed4f1ae30 100644
--- a/src/run_modules/post_processing.py
+++ b/src/run_modules/post_processing.py
@@ -21,7 +21,7 @@ from src.model_modules.linear_model import OrdinaryLeastSquaredModel
from src.model_modules.model_class import AbstractModelClass
from src.plotting.postprocessing_plotting import PlotMonthlySummary, PlotStationMap, PlotClimatologicalSkillScore, \
PlotCompetitiveSkillScore, PlotTimeSeries, PlotBootstrapSkillScore, PlotAvailability, PlotConditionalQuantiles
-from src.plotting.postprocessing_plotting import plot_conditional_quantiles
+# from src.plotting.postprocessing_plotting import plot_conditional_quantiles
from src.run_modules.run_environment import RunEnvironment
from typing import Dict
@@ -196,10 +196,6 @@ class PostProcessing(RunEnvironment):
if self.bootstrap_skill_scores is not None and "PlotBootstrapSkillScore" in plot_list:
PlotBootstrapSkillScore(self.bootstrap_skill_scores, plot_folder=self.plot_path, model_setup="CNN")
if "plot_conditional_quantiles" in plot_list:
- plot_conditional_quantiles(self.test_data.stations, pred_name="CNN", ref_name="obs",
- forecast_path=path, plot_name_affix="cali-ref", plot_folder=self.plot_path)
- plot_conditional_quantiles(self.test_data.stations, pred_name="obs", ref_name="CNN",
- forecast_path=path, plot_name_affix="like-bas", plot_folder=self.plot_path)
PlotConditionalQuantiles(self.test_data.stations, data_pred_path=path, plot_folder=self.plot_path)
if "PlotStationMap" in plot_list:
PlotStationMap(generators={'b': self.test_data}, plot_folder=self.plot_path)