From 07f4f2d61fec6e60b679a96722e94bd9668a8552 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Mon, 13 Nov 2023 15:15:33 -0500
Subject: [PATCH] Update plot_feature_detection.py

---
 examples/retrieve/plot_feature_detection.py | 62 +++++++++++++--------
 1 file changed, 40 insertions(+), 22 deletions(-)

diff --git a/examples/retrieve/plot_feature_detection.py b/examples/retrieve/plot_feature_detection.py
index f4489302b8..5122ce6928 100644
--- a/examples/retrieve/plot_feature_detection.py
+++ b/examples/retrieve/plot_feature_detection.py
@@ -277,7 +277,7 @@
     [
         (33 / 255, 145 / 255, 140 / 255),
         (253 / 255, 231 / 255, 37 / 255),
-        (210 / 255, 180 / 255, 140 / 255)
+        (210 / 255, 180 / 255, 140 / 255),
     ],
     N=3,
 )
@@ -351,7 +351,7 @@
     weakecho=3,
     bkgd_val=1,
     feat_val=2,
-    field='maxdz',
+    field="maxdz",
     estimate_flag=True,
     estimate_offset=5,
 )
@@ -430,7 +430,9 @@
 )
 
 # rescale reflectivity to snow rate
-grid = pyart.retrieve.qpe.ZtoR(grid, ref_field="reflectivity", a=57.3, b=1.67, save_name="snow_rate")
+grid = pyart.retrieve.qpe.ZtoR(
+    grid, ref_field="reflectivity", a=57.3, b=1.67, save_name="snow_rate"
+)
 
 # get dx dy
 dx = grid.x["data"][1] - grid.x["data"][0]
@@ -521,12 +523,8 @@
 ref_field_under["long_name"] = ref_field["long_name"] + " underestimate"
 
 # add to grid object
-grid.add_field(
-    "reflectivity_over", ref_field_over, replace_existing=True
-)
-grid.add_field(
-    "reflectivity_under", ref_field_under, replace_existing=True
-)
+grid.add_field("reflectivity_over", ref_field_over, replace_existing=True)
+grid.add_field("reflectivity_under", ref_field_under, replace_existing=True)
 
 # convert to snow rate
 grid = pyart.retrieve.qpe.ZtoR(
@@ -564,15 +562,23 @@
 fig, axs = plt.subplots(2, 2, figsize=(12, 10))
 # snow rate
 rpm = axs[0, 0].pcolormesh(
-    x / 1000, y / 1000, grid.fields["snow_rate"]["data"][0, :, :], vmin=0, vmax=10,
+    x / 1000,
+    y / 1000,
+    grid.fields["snow_rate"]["data"][0, :, :],
+    vmin=0,
+    vmax=10,
 )
 axs[0, 0].set_aspect("equal")
 axs[0, 0].set_title("Reflectivity [dBZ]")
 plt.colorbar(rpm, ax=axs[0, 0])
 # features best
 bpm = axs[0, 1].pcolormesh(
-    x / 1000, y / 1000, feature_estimate_dict["feature_detection"]["data"][0, :, :],
-    vmin=0, vmax=3, cmap=plt.get_cmap("viridis", 3),
+    x / 1000,
+    y / 1000,
+    feature_estimate_dict["feature_detection"]["data"][0, :, :],
+    vmin=0,
+    vmax=3,
+    cmap=plt.get_cmap("viridis", 3),
 )
 axs[0, 1].set_aspect("equal")
 axs[0, 1].set_title("Feature detection (best)")
@@ -580,8 +586,12 @@
 cb.ax.set_yticklabels(["Background", "Feature"])
 # features underestimate
 upm = axs[1, 0].pcolormesh(
-    x / 1000, y / 1000, feature_estimate_dict["feature_under"]["data"][0, :, :],
-    vmin=0, vmax=3, cmap=plt.get_cmap("viridis", 3),
+    x / 1000,
+    y / 1000,
+    feature_estimate_dict["feature_under"]["data"][0, :, :],
+    vmin=0,
+    vmax=3,
+    cmap=plt.get_cmap("viridis", 3),
 )
 axs[1, 0].set_aspect("equal")
 axs[1, 0].set_title("Feature detection (underestimate)")
@@ -589,8 +599,12 @@
 cb.ax.set_yticklabels(["Background", "Feature"])
 # features overestimate
 opm = axs[1, 1].pcolormesh(
-    x / 1000, y / 1000, feature_estimate_dict["feature_over"]["data"][0, :, :],
-    vmin=0, vmax=3, cmap=plt.get_cmap("viridis", 3),
+    x / 1000,
+    y / 1000,
+    feature_estimate_dict["feature_over"]["data"][0, :, :],
+    vmin=0,
+    vmax=3,
+    cmap=plt.get_cmap("viridis", 3),
 )
 axs[1, 1].set_aspect("equal")
 axs[1, 1].set_title("Feature detection (overestimate)")
@@ -636,8 +650,12 @@
 faint_features = scalar_features_masked - cosine_features_masked
 faint_features_masked = np.ma.masked_less_equal(faint_features, 0)
 # add strong and faint features
-features_faint_strong = 2 * faint_features_masked.filled(0) + cosine_features_masked.filled(0)
-features_faint_strong_masked = np.ma.masked_where(cosine_features_masked.mask, features_faint_strong)
+features_faint_strong = 2 * faint_features_masked.filled(
+    0
+) + cosine_features_masked.filled(0)
+features_faint_strong_masked = np.ma.masked_where(
+    cosine_features_masked.mask, features_faint_strong
+)
 
 # add to grid object
 new_dict = {
@@ -669,7 +687,7 @@
     [
         (32 / 255, 144 / 255, 140 / 255),
         (254 / 255, 238 / 255, 93 / 255),
-        (254 / 255, 152 / 255, 93 / 255)
+        (254 / 255, 152 / 255, 93 / 255),
     ],
     N=3,
 )
@@ -706,6 +724,7 @@
 # This last section shows how we can image mute the fields (Tomkins et al. 2022) and reduce the visual prominence of
 # melting and mixed precipitation.
 
+
 # create a function to quickly apply image muting to other fields
 def quick_image_mute(field, muted_ref):
     nonmuted_field = np.ma.masked_where(~muted_ref.mask, field)
@@ -713,6 +732,7 @@ def quick_image_mute(field, muted_ref):
 
     return nonmuted_field, muted_field
 
+
 # get fields
 snow_rate = grid.fields["snow_rate"]["data"][0, :, :]
 muted_ref = grid.fields["muted_reflectivity"]["data"][0, :, :]
@@ -734,9 +754,7 @@ def quick_image_mute(field, muted_ref):
 # plot
 fig, axs = plt.subplots(1, 2, figsize=(10, 4))
 # snow rate
-srpm = axs[0].pcolormesh(
-    x / 1000, y / 1000, snow_rate_nonmuted, vmin=0, vmax=10
-)
+srpm = axs[0].pcolormesh(x / 1000, y / 1000, snow_rate_nonmuted, vmin=0, vmax=10)
 srpmm = axs[0].pcolormesh(
     x / 1000, y / 1000, snow_rate_muted, vmin=0, vmax=10, cmap=muted_cmap
 )