From 7af7fe0f2391d8bfac6e184b369b13fb56115546 Mon Sep 17 00:00:00 2001
From: Azadeh Gholoubi <51101867+azadeh-gh@users.noreply.github.com>
Date: Tue, 16 Apr 2024 11:36:50 -0400
Subject: [PATCH] (1) Make O-F bias plot using symmetric colorbar and set cmap
 to bwr. (2) The index for studied_cycle_files (data frame subset from parent)
 does not alway start from zer0; therefore, assign index variable to get the
 actual index from parent dataframe (#148)

Co-authored-by: Emily Liu <emily.liu@noaa.gov>
---
 .../SpatialTemporalStats.py                   | 26 ++++++++++++++++---
 1 file changed, 22 insertions(+), 4 deletions(-)

diff --git a/ush/SpatialTemporalStatsTool/SpatialTemporalStats.py b/ush/SpatialTemporalStatsTool/SpatialTemporalStats.py
index c2a7244..ed28cb1 100644
--- a/ush/SpatialTemporalStatsTool/SpatialTemporalStats.py
+++ b/ush/SpatialTemporalStatsTool/SpatialTemporalStats.py
@@ -227,22 +227,36 @@ def plot_obs(self, selected_var_gdf, var_name, region, resolution, output_path):
                     selected_var_gdf.geometry.apply(lambda geom: geom.centroid.y < -60)
                 ]
 
-            min_val, max_val, std_val = (
+            min_val, max_val, std_val, avg_val = (
                 filtered_gdf[item].min(),
                 filtered_gdf[item].max(),
                 filtered_gdf[item].std(),
+                filtered_gdf[item].mean(),
             )
-            cbar_label = "grid=%dX%d, min=%.3lf, max=%.3lf, std=%.3lf\n" % (
+
+            if item == "Obs_Minus_Forecast_adjusted_Average":
+                max_val_cbar = 5.0*std_val
+                min_val_cbar = -5.0*std_val
+                cmap = "bwr"
+            else:
+                max_val_cbar = max_val
+                min_val_cbar = min_val
+                cmap = "jet"
+
+            cbar_label = "grid=%dx%d,   min=%.3lf,   max=%.3lf,   bias=%.3lf,   std=%.3lf\n" % (
                 resolution,
                 resolution,
                 min_val,
                 max_val,
+                avg_val,
                 std_val,
             )
 
             filtered_gdf.plot(
                 ax=ax,
-                cmap="jet",
+                cmap=cmap,
+                vmin=min_val_cbar,
+                vmax=max_val_cbar,
                 column=item,
                 legend=True,
                 missing_kwds={"color": "lightgrey"},
@@ -333,11 +347,15 @@ def make_summary_plots(
                 & ((files_date_times_df["date"] <= end_date))
             )
         ]["file_name"]
+        index = studied_cycle_files.index
+
         Summary_results = []
 
         # get unique channels from one of the files
-        ds = xarray.open_dataset(studied_cycle_files[0])
+        ds = xarray.open_dataset(studied_cycle_files[index[0]])
         unique_channels = np.unique(ds["Channel_Index"].data).tolist()
+        print('Total Number of Channels ', len(unique_channels))
+        print('Channels ', unique_channels)
 
         for this_channel in unique_channels:
             this_channel_values = np.empty(shape=(0,))