diff --git a/src/plotting/postprocessing_plotting.py b/src/plotting/postprocessing_plotting.py index 8c8ea98e9f356be0b9064afcfcab73d00df67311..eb3f7f8c058fae47c2703e6e53ee22bdc013f7a2 100644 --- a/src/plotting/postprocessing_plotting.py +++ b/src/plotting/postprocessing_plotting.py @@ -22,7 +22,6 @@ logging.getLogger('matplotlib').setLevel(logging.WARNING) def plot_monthly_summary(stations, data_path, name: str, window_lead_time, target_var, plot_folder="."): - logging.debug("run plot_monthly_summary()") forecasts = None @@ -40,17 +39,18 @@ def plot_monthly_summary(stations, data_path, name: str, window_lead_time, targe data_concat = xr.concat([data_orig, data_cnn], dim="ahead") data_concat = data_concat.drop("type") - data_concat.index.values = data_concat.index.values.astype("datetime64[M]").astype(int) %12 +1 + data_concat.index.values = data_concat.index.values.astype("datetime64[M]").astype(int) % 12 + 1 data_concat = data_concat.clip(min=0) forecasts = xr.concat([forecasts, data_concat], 'index') if forecasts is not None else data_concat forecasts = forecasts.to_dataset(name='values').to_dask_dataframe() logging.debug("... start plotting") - ax = sns.boxplot(x='index', y='values', hue='ahead', data=forecasts.compute(), whis=1., - palette=[matplotlib.colors.cnames["green"]]+sns.color_palette("Blues_d", window_lead_time).as_hex(), - flierprops={'marker': '.', 'markersize': 1}, - showmeans=True, meanprops={'markersize': 1, 'markeredgecolor': 'k'}) + ax = sns.boxplot(x='index', y='values', hue='ahead', data=forecasts.compute(), whis=1., + palette=[matplotlib.colors.cnames["green"]] + sns.color_palette("Blues_d", + window_lead_time).as_hex(), + flierprops={'marker': '.', 'markersize': 1}, showmeans=True, + meanprops={'markersize': 1, 'markeredgecolor': 'k'}) ax.set(xlabel='month', ylabel=f'{target_var}') plt.tight_layout() plot_path = os.path.join(os.path.abspath(plot_folder), 'test_monthly_box.pdf') @@ -64,7 +64,6 @@ def plot_climsum_boxplot(): def plot_station_map(generators, plot_folder="."): - logging.debug("run station_map()") fig = plt.figure(figsize=(10, 5)) ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) @@ -80,7 +79,8 @@ def plot_station_map(generators, plot_folder="."): for k, v in enumerate(gen): station_coords = gen.get_data_generator(k).meta.loc[['station_lon', 'station_lat']] station_names = gen.get_data_generator(k).meta.loc[['station_id']] - IDx, IDy = float(station_coords.loc['station_lon'].values), float(station_coords.loc['station_lat'].values) + IDx, IDy = float(station_coords.loc['station_lon'].values), float( + station_coords.loc['station_lat'].values) ax.plot(IDx, IDy, mfc=color, mec='k', marker='s', markersize=6, transform=ccrs.PlateCarree()) plot_path = os.path.join(os.path.abspath(plot_folder), 'test_map_plot.pdf') @@ -89,7 +89,8 @@ def plot_station_map(generators, plot_folder="."): def plot_conditional_quantiles(stations: list, plot_folder: str = ".", rolling_window: int = 3, ref_name: str = 'orig', - pred_name: str = 'CNN', season: str = "", forecast_path: str = None): + 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 @@ -102,6 +103,8 @@ def plot_conditional_quantiles(stations: list, plot_folder: str = ".", rolling_w :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()") @@ -146,15 +149,24 @@ def plot_conditional_quantiles(stations: list, plot_folder: str = ".", rolling_w return quantile_panel - opts = {"q": [.1, .25, .5, .75, .9], - "linetype": [':', '-.', '--', '-.', ':'], + def labels(plot_type, data_unit="ppb"): + names = (f"forecast concentration (in {data_unit})", f"observed concentration (in {data_unit})") + if plot_type == "orig": + 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": "forecast concentration (in ppb)", - "ylabel": "observed concentration (in ppb)"} + "xlabel": xlabel, "ylabel": ylabel} # set name and path of the plot - plot_name = f"test_conditional_quantiles{f'_{season}' if len(season) > 0 else ''}" - plot_path = os.path.join(os.path.abspath(plot_folder), f"{plot_name}_plot.pdf") + base_name = "conditional_quantiles" + 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: @@ -182,10 +194,10 @@ def plot_conditional_quantiles(stations: list, plot_folder: str = ".", rolling_w 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) + 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') + 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') @@ -197,8 +209,8 @@ def plot_conditional_quantiles(stations: list, plot_folder: str = ".", rolling_w 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)) + 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 diff --git a/src/run_modules/post_processing.py b/src/run_modules/post_processing.py index b87c637055b594005cecbfa391e137ce103e0f32..d773c9789619c802f05f3d922d554fb1b246ffe0 100644 --- a/src/run_modules/post_processing.py +++ b/src/run_modules/post_processing.py @@ -57,7 +57,10 @@ class PostProcessing(RunEnvironment): window_lead_time = self.data_store.get("window_lead_time", "general") target_var = self.data_store.get("target_var", "general") - plot_conditional_quantiles(self.test_data.stations, plot_folder=self.plot_path, forecast_path=self.data_store.get("forecast_path", "general")) + plot_conditional_quantiles(self.test_data.stations, plot_folder=self.plot_path, pred_name="CNN", ref_name="orig", + forecast_path=self.data_store.get("forecast_path", "general"), plot_name_affix="cali-ref") + plot_conditional_quantiles(self.test_data.stations, plot_folder=self.plot_path, pred_name="orig", ref_name="CNN", + forecast_path=self.data_store.get("forecast_path", "general"), plot_name_affix="like-bas") plot_station_map(generators={'b': self.test_data}, plot_folder=self.plot_path) plot_monthly_summary(self.test_data.stations, path, r"forecasts_%s_test.nc", window_lead_time, target_var, plot_folder=self.plot_path)