updated the python codes to modify the user settings and also the map…

…plots format

ush/SpatialTemporalStatsTool/SpatialTemporalStats.py

@@ -1,6 +1,8 @@
 import os
 from datetime import datetime
 
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
 import geopandas as gpd
 import matplotlib.pyplot as plt
 import numpy as np
@@ -177,55 +179,72 @@ def plot_obs(self, selected_var_gdf, var_name, region, resolution, output_path):
 
         for _, item in enumerate(var_names):
             plt.figure(figsize=(12, 8))
-            ax = plt.subplot(1, 1, 1)
+            ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree())
+            # Add global map coastlines
+            ax.add_feature(cfeature.GSHHSFeature(scale="auto"))
+            filtered_gdf = selected_var_gdf.copy()
 
             if region == 1:
                 # Plotting global region (no need for filtering)
-                title = "Global Region"
-                filtered_gdf = selected_var_gdf
+                title = "Global"
+                # filtered_gdf = selected_var_gdf
 
             elif region == 2:
                 # Plotting polar region (+60 latitude and above)
                 title = "Polar Region (+60 latitude and above)"
-                filtered_gdf = selected_var_gdf[
-                    selected_var_gdf.geometry.apply(
+                filtered_gdf[item] = np.where(
+                    filtered_gdf.geometry.apply(
                         lambda geom: self.is_polygon_in_polar_region(geom, 60)
-                    )
-                ]
+                    ),
+                    filtered_gdf[item],
+                    np.nan,
+                )
 
             elif region == 3:
                 # Plotting northern mid-latitudes region (20 to 60 latitude)
                 title = "Northern Mid-latitudes Region (20 to 60 latitude)"
-                filtered_gdf = selected_var_gdf[
-                    selected_var_gdf.geometry.apply(
+                filtered_gdf[item] = np.where(
+                    filtered_gdf.geometry.apply(
                         lambda geom: self.is_polygon_in_latitude_range(geom, 20, 60)
-                    )
-                ]
+                    ),
+                    filtered_gdf[item],
+                    np.nan,
+                )
 
             elif region == 4:
                 # Plotting tropics region (-20 to 20 latitude)
                 title = "Tropics Region (-20 to 20 latitude)"
-                filtered_gdf = selected_var_gdf[
-                    selected_var_gdf.geometry.apply(
+
+                filtered_gdf[item] = np.where(
+                    filtered_gdf.geometry.apply(
                         lambda geom: self.is_polygon_in_latitude_range(geom, -20, 20)
-                    )
-                ]
+                    ),
+                    filtered_gdf[item],
+                    np.nan,
+                )
 
             elif region == 5:
                 # Plotting southern mid-latitudes region (-60 to -20 latitude)
                 title = "Southern Mid-latitudes Region (-60 to -20 latitude)"
-                filtered_gdf = selected_var_gdf[
-                    selected_var_gdf.geometry.apply(
+                filtered_gdf[item] = np.where(
+                    filtered_gdf.geometry.apply(
                         lambda geom: self.is_polygon_in_latitude_range(geom, -60, -20)
-                    )
-                ]
+                    ),
+                    filtered_gdf[item],
+                    np.nan,
+                )
 
             elif region == 6:
                 # Plotting southern polar region (less than -60 latitude)
                 title = "Southern Polar Region (less than -60 latitude)"
-                filtered_gdf = selected_var_gdf[
-                    selected_var_gdf.geometry.apply(lambda geom: geom.centroid.y < -60)
-                ]
+                filtered_gdf[item] = np.where(
+                    filtered_gdf.geometry.apply(lambda geom: geom.centroid.y < -60),
+                    filtered_gdf[item],
+                    np.nan,
+                )
+                # filtered_gdf = selected_var_gdf[
+                #     selected_var_gdf.geometry.apply(lambda geom: geom.centroid.y < -60)
+                # ]
 
             min_val, max_val, std_val, avg_val = (
                 filtered_gdf[item].min(),
@@ -270,6 +289,14 @@ def plot_obs(self, selected_var_gdf, var_name, region, resolution, output_path):
                 },
             )
 
+            filtered_gdf.to_file(
+                os.path.join(
+                    output_path,
+                    "%s_ch%d_%s_region_%d.gpkg"
+                    % (self.sensor, self.channel_no, item, region),
+                )
+            )
+
             plt.title("%s\n%s ch:%d %s" % (title, self.sensor, self.channel_no, item))
             plt.savefig(
                 os.path.join(
@@ -364,7 +391,11 @@ def make_summary_plots(
         for this_cycle_obs_file in studied_cycle_files:
             ds = xarray.open_dataset(this_cycle_obs_file)
             if QC_filter:
-                QC_bool = ds["QC_Flag"].data == 0
+                QC_bool = ds["QC_Flag"].data == 0.0
+            else:
+                QC_bool = np.ones(
+                    ds["QC_Flag"].data.shape, dtype=bool
+                )  # this selects all obs as True
             for this_channel in unique_channels:
                 channel_bool = ds["Channel_Index"].data == this_channel
 
@@ -379,7 +410,7 @@ def make_summary_plots(
             this_channel_values = Allchannels_data[this_channel]
             squared_values = [x**2 for x in this_channel_values]
             mean_of_squares = sum(squared_values) / len(squared_values)
-            rms_value = mean_of_squares ** 0.5
+            rms_value = mean_of_squares**0.5
             Summary_results.append(
                 [
                     this_channel,
@@ -391,7 +422,8 @@ def make_summary_plots(
             )
 
         Summary_resultsDF = pd.DataFrame(
-            Summary_results, columns=["channel", "count", "std", "mean", "rms"])
+            Summary_results, columns=["channel", "count", "std", "mean", "rms"]
+        )
         # Plotting
         plt.figure(figsize=(10, 6))
         plt.scatter(Summary_resultsDF["channel"], Summary_resultsDF["count"], s=50)
@@ -442,3 +474,4 @@ def make_summary_plots(
         )
 
         return Summary_resultsDF
+