cali-ref and like-bas plots can be created now

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

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)