From 7a396822cc99894284c125264b6d51c9d363e06b Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Wed, 4 Oct 2023 11:48:50 -0400
Subject: [PATCH 01/21] Rename functions and variables to remove
 "convective-stratiform"

We want to make this function more generic so it can be used on a wide variety of fields and not just to detect convective and stratiform elements. I renamed most of the functions and some of the variables to reflect this and hopefully make it easier to use.
---
 pyart/retrieve/__init__.py    |   2 +-
 pyart/retrieve/_echo_class.py | 445 ++++++++++++++++------------------
 pyart/retrieve/echo_class.py  | 198 +++++++--------
 3 files changed, 304 insertions(+), 341 deletions(-)

diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py
index 4666f6913d..41811f33d7 100644
--- a/pyart/retrieve/__init__.py
+++ b/pyart/retrieve/__init__.py
@@ -5,7 +5,7 @@
 
 from .advection import grid_displacement_pc, grid_shift  # noqa
 from .comp_z import composite_reflectivity  # noqa
-from .echo_class import conv_strat_yuter  # noqa
+from .echo_class import feature_detection  # noqa
 from .echo_class import get_freq_band  # noqa
 from .echo_class import hydroclass_semisupervised  # noqa
 from .echo_class import steiner_conv_strat  # noqa
diff --git a/pyart/retrieve/_echo_class.py b/pyart/retrieve/_echo_class.py
index ca9ff8209a..b3dadce2a3 100644
--- a/pyart/retrieve/_echo_class.py
+++ b/pyart/retrieve/_echo_class.py
@@ -239,8 +239,8 @@ def steiner_class_buff(
     return sclass
 
 
-def _revised_conv_strat(
-    refl,
+def _feature_detection(
+    field,
     dx,
     dy,
     always_core_thres=42,
@@ -252,48 +252,47 @@ def _revised_conv_strat(
     use_addition=True,
     calc_thres=0.75,
     weak_echo_thres=5.0,
-    min_dBZ_used=5.0,
+    min_val_used=5.0,
     dB_averaging=True,
     remove_small_objects=True,
     min_km2_size=10,
-    val_for_max_conv_rad=30,
-    max_conv_rad_km=5.0,
-    cs_core=3,
+    val_for_max_rad=30,
+    max_rad_km=5.0,
+    core=3,
     nosfcecho=0,
     weakecho=3,
-    sf=1,
-    conv=2,
+    bkgd_val=1,
+    feat_val=2,
 ):
     """
-    We perform the Yuter and Houze (1997) algorithm for echo classification
-    using only the reflectivity field in order to classify each grid point
-    as either convective, stratiform or undefined. Grid points are
+    These functions are used to detect features in fields based on how distinct they are from the background
+    average. Methodology described in Tomkins et al. (2023), originally based on convective-stratiform algorithms
+    developed by Steiner et al. (1995), Yuter and Houze (1997), and Yuter et al. (2005) Grid points are
     classified as follows,
 
     nosfcecho = No Surface Echo/ Undefined
-    sf = Stratiform
-    conv = Convective
+    bkgd_val = Background echo (e.g. Stratiform)
+    feat_val = Feature (e.g. Convective)
     weakecho = Weak Echo
 
-    refl : array
-        array of reflectivity values
+    field : array
+        array of values to find features
     x, y : array
-        x and y coordinates of reflectivity array, respectively
+        x and y coordinates of field array, respectively
     dx, dy : float
         The x- and y-dimension resolutions in meters, respectively.
     always_core_thres : float, optional
-        Threshold for points that are always convective. All values above the threshold are classifed as convective
+        Threshold for points that are always features. All values above the threshold are classified as features.
     bkg_rad_km : float, optional
-        Radius to compute background reflectivity in kilometers. Default is 11 km. Recommended to be at least 3 x
-        grid spacing
+        Radius to compute background field in kilometers. Default is 11 km. Recommended to be at least 3 x grid spacing
     use_cosine : bool, optional
-        Boolean used to determine if cosine scheme should be used for identifying convective cores (True) or a scalar
-        scheme (False)
+        Boolean used to determine if a cosine scheme (see Yuter and Houze (1997)) should be used for identifying
+        cores (True) or if a simpler scalar scheme (False) should be used.
     max_diff : float, optional
-        Maximum difference between background average and reflectivity in order to be classified as convective.
+        Maximum difference between background average and grid value in order to be classified as features.
         "a" value in Eqn. B1 in Yuter and Houze (1997)
     zero_diff_cos_val : float, optional
-        Value where difference between background average and reflectivity is zero in the cosine function
+        Value where difference between background average and grid value is zero in the cosine function
         "b" value in Eqn. B1 in Yuter and Houze (1997)
     scalar_diff : float, optional
         If using a scalar difference scheme, this value is the multiplier or addition to the background average
@@ -303,50 +302,48 @@ def _revised_conv_strat(
         Minimum percentage of points needed to be considered in background average calculation
     weak_echo_thres : float, optional
         Threshold for determining weak echo. All values below this threshold will be considered weak echo
-    min_dBZ_used : float, optional
-        Minimum dBZ value used for classification. All values below this threshold will be considered no surface echo
+    min_val_used : float, optional
+        Minimum value used for classification. All values below this threshold will be considered no surface echo
+        See Yuter and Houze (1997) and Yuter et al. (2005) for more detail. Units based on input field
     dB_averaging : bool, optional
         True if using dBZ values that need to be converted to linear Z before averaging. False for other types of values
     remove_small_objects : bool, optional
-        Determines if small objects should be removed from convective core array. Default is True.
+        Determines if small objects should be removed from core array. Default is True.
     min_km2_size : float, optional
-        Minimum size of convective cores to be considered. Cores less than this size will be removed. Default is 10
-        km^2.
-    val_for_max_conv_rad : float, optional
-        dBZ for maximum convective radius. Convective cores with values above this will have the maximum convective
-        radius
-    max_conv_rad_km : float, optional
-        Maximum radius around convective cores to classify as convective. Default is 5 km.
-    cs_core : int, optional
-        Value for points classified as convective cores
+        Minimum size of Cores to be considered. Cores less than this size will be removed. Default is 10 km^2.
+    val_for_max_rad : float, optional
+        value used for maximum radius. Cores with values above this will have the maximum radius incorporated.
+    max_rad_km : float, optional
+        Maximum radius around cores to classify as feature. Default is 5 km
+    core : int, optional
+        Value for points classified as cores
     nosfcecho : int, optional
-        Value for points classified as no surface echo, based on min_dBZ_used
+        Value for points classified as no surface echo, based on min_val_used
     weakecho : int, optional
-        Value for points classified as weak echo, based on weak_echo_thres
-    sf : int, optional
-        Value for points classified as stratiform
-    conv : int, optional
-        Value for points classified as convective
+        Value for points classified as weak echo, based on weak_echo_thres.
+    bkgd_val : int, optional
+        Value for points classified as background echo.
+    feat_val : int, optional
+        Value for points classified as features.
 
     Returns
     -------
-    refl_bkg : array
+    field_bkg : array
         Array of background values
-    conv_core_array : array
-        Array of initial convective cores (identified convective elements without convective radii applied)
-    conv_strat_array : array
-        Array of convective stratiform classifcation with convective radii applied
+    core_array : array
+        Array of initial cores (identified features without radii applied)
+    feature_array : array
+        Array of feature detection with radii applied
     """
 
-    # Set up mask arrays for background average and
-    # prepare for convective mask arrays
-    # calculate maximum convective diameter from max. convective radius (input)
-    max_conv_diameter = int(np.floor((max_conv_rad_km / (dx / 1000)) * 2))
+    # Set up mask arrays for background average and prepare for mask arrays
+    # calculate maximum diameter from max. radius (input)
+    max_diameter = int(np.floor((max_rad_km / (dx / 1000)) * 2))
     # if diameter is even, make odd
-    if max_conv_diameter % 2 == 0:
-        max_conv_diameter = max_conv_diameter + 1
+    if max_diameter % 2 == 0:
+        max_diameter = max_diameter + 1
     # find center point
-    center_conv_mask_x = int(np.floor(max_conv_diameter / 2))
+    center_mask_x = int(np.floor(max_diameter / 2))
 
     # prepare background mask array for computing background average
     # calculate number of pixels for background array given requested background radius and dx
@@ -370,42 +367,30 @@ def _revised_conv_strat(
         circular=True,
     )
 
-    # Convective stratiform detection
-    # start by making reflectivity a masked array
-    refl = np.ma.masked_invalid(refl)
+    # Feature detection
+    # start by making field a masked array
+    field = np.ma.masked_invalid(field)
     # Compute background radius
-    refl_bkg = calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres)
+    field_bkg = calc_bkg_intensity(field, bkg_mask_array, dB_averaging, calc_thres)
 
-    # Get convective core array from cosine scheme, or scalar scheme
+    # Get core array from cosine scheme, or scalar scheme
     if use_cosine:
-        conv_core_array = convcore_cos_scheme(
-            refl, refl_bkg, max_diff, zero_diff_cos_val, always_core_thres, cs_core
-        )
+        core_array = core_cos_scheme(field, field_bkg, max_diff, zero_diff_cos_val, always_core_thres, core)
     else:
-        conv_core_array = convcore_scalar_scheme(
-            refl,
-            refl_bkg,
-            scalar_diff,
-            always_core_thres,
-            cs_core,
-            use_addition=use_addition,
-        )
+        core_array = core_scalar_scheme(field, field_bkg, scalar_diff, always_core_thres, core,
+                                        use_addition=use_addition)
 
-    # Assign convective radii based on background reflectivity
-    conv_radius_km = assign_conv_radius_km(
-        refl_bkg,
-        val_for_max_conv_rad=val_for_max_conv_rad,
-        max_conv_rad=max_conv_rad_km,
-    )
+    # Assign radii based on background field
+    radius_array_km = assign_feature_radius_km(field_bkg, val_for_max_rad=val_for_max_rad, max_rad=max_rad_km)
 
-    # remove small objects in convective core array
+    # remove small objects in core array
     if remove_small_objects:
         # calculate minimum pixel size given dx and dy
         min_pix_size = min_km2_size / ((dx / 1000) * (dy / 1000))
-        # label connected objects in convective core array
-        cc_labels, _ = scipy.ndimage.label(conv_core_array)
-        # mask labels where convective core array is masked
-        cc_labels = np.ma.masked_where(conv_core_array.mask, cc_labels)
+        # label connected objects in core array
+        cc_labels, _ = scipy.ndimage.label(core_array)
+        # mask labels where core array is masked
+        cc_labels = np.ma.masked_where(core_array.mask, cc_labels)
 
         # loop through each unique label
         for lab in np.unique(cc_labels):
@@ -413,51 +398,41 @@ def _revised_conv_strat(
             size_lab = np.count_nonzero(cc_labels == lab)
             # if number of pixels is less than minimum, then remove core
             if size_lab < min_pix_size:
-                conv_core_array[cc_labels == lab] = 0
+                core_array[cc_labels == lab] = 0
 
-    # Incorporate convective radius using binary dilation
+    # Incorporate radius using binary dilation
     # Create empty array for assignment
-    temp_assignment = np.zeros_like(conv_core_array)
+    temp_assignment = np.zeros_like(core_array)
 
     # Loop through radii
-    for radius in np.arange(1, max_conv_rad_km + 1):
+    for radius in np.arange(1, max_rad_km + 1):
         # create mask array for radius incorporation
-        conv_mask_array = create_conv_radius_mask(
-            max_conv_diameter, radius, dx / 1000, dy / 1000, center_conv_mask_x
+        radius_mask_array = create_radius_mask(
+            max_diameter, radius, dx / 1000, dy / 1000, center_mask_x
         )
         # find location of radius
-        temp = conv_radius_km == radius
+        temp = radius_array_km == radius
         # get cores for given radius
-        temp_core = np.ma.masked_where(~temp, conv_core_array)
+        temp_core = np.ma.masked_where(~temp, core_array)
         # dilate cores
         temp_dilated = scipy.ndimage.binary_dilation(
-            temp_core.filled(0), conv_mask_array
+            temp_core.filled(0), radius_mask_array
         )
         # add to assignment array
         temp_assignment = temp_assignment + temp_dilated
 
     # add dilated cores to original array
-    conv_core_copy = np.ma.copy(conv_core_array)
-    conv_core_copy[temp_assignment >= 1] = cs_core
+    core_copy = np.ma.copy(core_array)
+    core_copy[temp_assignment >= 1] = core
 
-    # Now do convective stratiform classification
-    conv_strat_array = np.zeros_like(refl)
-    conv_strat_array = classify_conv_strat_array(
-        refl,
-        conv_strat_array,
-        conv_core_copy,
-        nosfcecho,
-        conv,
-        sf,
-        weakecho,
-        cs_core,
-        min_dBZ_used,
-        weak_echo_thres,
-    )
-    # mask where reflectivity is masked
-    conv_strat_array = np.ma.masked_where(refl.mask, conv_strat_array)
+    # Now do feature detection
+    feature_array = np.zeros_like(field)
+    feature_array = classify_feature_array(field, feature_array, core_copy, nosfcecho, feat_val, bkgd_val, weakecho,
+                                           core, min_val_used, weak_echo_thres)
+    # mask where field is masked
+    feature_array = np.ma.masked_where(field.mask, feature_array)
 
-    return refl_bkg, conv_core_array, conv_strat_array
+    return field_bkg, core_array, feature_array
 
 
 # functions
@@ -524,15 +499,15 @@ def create_radial_mask(
     return mask_array
 
 
-def calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres=None):
+def calc_bkg_intensity(field, bkg_mask_array, dB_averaging, calc_thres=None):
     """
-    Calculate the background of the given refl array. The footprint used to
+    Calculate the background of the given field. The footprint used to
     calculate the average for each pixel is given by bkg_mask_array
 
     Parameters
     ----------
-    refl : array
-        Reflectivity array to compute average
+    field : array
+        Field array to compute average
     bkg_mask_array : array
         Array of radial points to use for average
     dB_averaging : bool
@@ -542,17 +517,17 @@ def calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres=None):
 
     Returns
     -------
-    refl_bkg : array
+    field_bkg : array
         Array of average values
     """
 
     # if dBaverage is true, convert reflectivity to linear Z
     if dB_averaging:
-        refl = 10 ** (refl / 10)
+        field = 10 ** (field / 10)
 
-    # calculate background reflectivity with circular footprint
-    refl_bkg = scipy.ndimage.generic_filter(
-        refl.filled(np.nan),
+    # calculate background field with circular footprint
+    field_bkg = scipy.ndimage.generic_filter(
+        field.filled(np.nan),
         function=np.nanmean,
         mode="constant",
         footprint=bkg_mask_array.astype(bool),
@@ -562,8 +537,8 @@ def calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres=None):
     # if calc_thres is not none, then calculate the number of points used to calculate average
     if calc_thres is not None:
         # count valid points
-        refl_count = scipy.ndimage.generic_filter(
-            refl.filled(0),
+        field_count = scipy.ndimage.generic_filter(
+            field.filled(0),
             function=np.count_nonzero,
             mode="constant",
             footprint=bkg_mask_array.astype(bool),
@@ -572,223 +547,223 @@ def calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres=None):
         # find threshold number of points
         val = calc_thres * np.count_nonzero(bkg_mask_array)
         # mask out values
-        refl_bkg = np.ma.masked_where(refl_count < val, refl_bkg)
+        field_bkg = np.ma.masked_where(field_count < val, field_bkg)
 
-    # mask where original reflectivity is invalid
-    refl_bkg = np.ma.masked_where(refl.mask, refl_bkg)
+    # mask where original field is invalid
+    field_bkg = np.ma.masked_where(field.mask, field_bkg)
 
     # if dBaveraging is true, convert background reflectivity to dBZ
     if dB_averaging:
-        refl_bkg = 10 * np.log10(refl_bkg)
-        # mask where original reflectivity is invalid
-        refl_bkg = np.ma.masked_where(refl.mask, refl_bkg)
+        field_bkg = 10 * np.log10(field_bkg)
+        # mask where original field is invalid
+        field_bkg = np.ma.masked_where(field.mask, field_bkg)
 
-    return refl_bkg
+    return field_bkg
 
 
-def convcore_cos_scheme(
-    refl, refl_bkg, max_diff, zero_diff_cos_val, always_core_thres, CS_CORE
+def core_cos_scheme(
+    field, field_bkg, max_diff, zero_diff_cos_val, always_core_thres, CS_CORE
 ):
     """
-    Function for assigning convective cores based on a cosine function
+    Function for assigning cores based on a cosine function
 
     Parameters
     ----------
-    refl : array
-        Reflectivity values
-    refl_bkg : array
-        Background average of reflectivity values
+    field : array
+        Field values
+    field_bkg : array
+        Background average of field values
     max_diff : float
-        Maximum difference between refl and refl_bkg needed for convective classification
+        Maximum difference between field and field_bkg needed for feature detection
     zero_diff_cos_val : float
         Value where the cosine function returns a zero difference
     always_core_thres : float
-        All values above this threshold considered to be convective
+        All values above this threshold considered to be features
     CS_CORE : int
-        Value assigned to convective pixels
+        Value assigned to features
 
     Returns
     -------
-    conv_core_array : array
-        Array of booleans if point is convective (1) or not (0)
+    core_array : array
+        Array of booleans if point is feature (1) or not (0)
     """
 
-    # initialize entire array to not a convective core
-    conv_core_array = np.zeros_like(refl)
+    # initialize entire array to not a core
+    core_array = np.zeros_like(field)
 
     # calculate zeDiff for entire array
-    zDiff = max_diff * np.cos((np.pi * refl_bkg) / (2 * zero_diff_cos_val))
+    zDiff = max_diff * np.cos((np.pi * field_bkg) / (2 * zero_diff_cos_val))
     zDiff[zDiff < 0] = 0  # where difference less than zero, set to zero
-    zDiff[refl_bkg < 0] = max_diff  # where background less than zero, set to max. diff
+    zDiff[field_bkg < 0] = max_diff  # where background less than zero, set to max. diff
 
     # set core values
-    # where refl >= always_core_thres and where difference exceeds minimum
+    # where field >= always_core_thres and where difference exceeds minimum
     core_elements = np.logical_or(
-        (refl >= always_core_thres), (refl - refl_bkg) >= zDiff
+        (field >= always_core_thres), (field - field_bkg) >= zDiff
     )
     core_elements = core_elements.filled(0)
-    conv_core_array[core_elements] = CS_CORE
+    core_array[core_elements] = CS_CORE
 
-    # mask by refl array
-    conv_core_array = np.ma.masked_where(refl.mask, conv_core_array)
+    # mask by field array
+    core_array = np.ma.masked_where(field.mask, core_array)
 
-    return conv_core_array
+    return core_array
 
 
-def convcore_scalar_scheme(
-    refl, refl_bkg, max_diff, always_core_thres, CS_CORE, use_addition=False
+def core_scalar_scheme(
+    field, field_bkg, max_diff, always_core_thres, CORE, use_addition=False
 ):
     """
-    Function for assigning convective cores based on a scalar difference
+    Function for assigning cores based on a scalar difference
 
     Parameters
     ----------
-    refl : array
-        Reflectivity values
-    refl_bkg : array
-        Background average of reflectivity values
+    field : array
+        Field values
+    field_bkg : array
+        Background average of field values
     max_diff : float
-        Maximum difference between refl and refl_bkg needed for convective classification
+        Maximum difference between field and field_bkg needed for feature detection
     always_core_thres : float
-        All values above this threshold considered to be convective
-    CS_CORE : int
-        Value assigned to convective pixels
+        All values above this threshold considered to be a feature
+    CORE : int
+        Value assigned to features
     use_addition : bool
         Boolean to determine if scalar should be added (True) or multiplied (False)
 
     Returns
     -------
-    conv_core_array : array
-        Array of booleans if point is convective (1) or not (0)
+    core_array : array
+        Array of booleans if point is feature (1) or not (0)
     """
 
-    # initialize entire array to not a convective core
-    conv_core_array = np.zeros_like(refl)
+    # initialize entire array to not a core
+    core_array = np.zeros_like(field)
 
     # calculate zDiff for entire array
     # if addition, add difference. Else, multiply difference
     if use_addition:
-        zDiff = (max_diff + refl_bkg) - refl_bkg
+        zDiff = (max_diff + field_bkg) - field_bkg
     else:
-        zDiff = (max_diff * refl_bkg) - refl_bkg
+        zDiff = (max_diff * field_bkg) - field_bkg
 
     zDiff[zDiff < 0] = 0  # where difference less than zero, set to zero
-    zDiff[refl_bkg < 0] = 0  # where background less than zero, set to zero
+    zDiff[field_bkg < 0] = 0  # where background less than zero, set to zero
 
     # set core values
-    # where refl >= always_core_thres and where difference exceeds minimum
+    # where field >= always_core_thres and where difference exceeds minimum
     core_elements = np.logical_or(
-        (refl >= always_core_thres), (refl - refl_bkg) >= zDiff
+        (field >= always_core_thres), (field - field_bkg) >= zDiff
     )
     core_elements = core_elements.filled(0)
-    conv_core_array[core_elements] = CS_CORE
+    core_array[core_elements] = CORE
 
-    # mask by refl array
-    conv_core_array = np.ma.masked_where(refl.mask, conv_core_array)
+    # mask by field array
+    core_array = np.ma.masked_where(field.mask, core_array)
 
-    return conv_core_array
+    return core_array
 
 
-def create_conv_radius_mask(
-    max_conv_diameter, radius_km, x_spacing, y_spacing, center_conv_mask_x
+def create_radius_mask(
+    max_diameter, radius_km, x_spacing, y_spacing, center_mask_x
 ):
     """
-    Does and initial convective stratiform classification
+    Creates a circular mask based on input diameter
 
     Parameters
     ----------
-    max_conv_diameter : int
-        maximum convective diameter in kilometers
+    max_diameter : int
+        maximum diameter in kilometers
     radius_km : int
-        convective radius in kilometers
+        radius in kilometers
     x_spacing, y_spacing : float
         x- and y-dimension pixel size in meters, respectively
-    center_conv_mask_x : int
+    center_mask_x : int
         index of center point
 
     Returns
     -------
-    conv_mask_array : array
-        array masked based on distance of convective diameter
+    feature_mask_array : array
+        array masked based on distance of diameter to incorporate
     """
 
-    conv_mask_array = np.zeros((max_conv_diameter, max_conv_diameter))
-    conv_mask_array = create_radial_mask(
-        conv_mask_array,
+    feature_mask_array = np.zeros((max_diameter, max_diameter))
+    feature_mask_array = create_radial_mask(
+        feature_mask_array,
         0,
         radius_km,
         x_spacing,
         y_spacing,
-        center_conv_mask_x,
-        center_conv_mask_x,
+        center_mask_x,
+        center_mask_x,
         True,
     )
 
-    return conv_mask_array
+    return feature_mask_array
 
 
-def assign_conv_radius_km(refl_bkg, val_for_max_conv_rad, max_conv_rad=5):
+def assign_feature_radius_km(field_bkg, val_for_max_rad, max_rad=5):
     """
-    Assigns the convective radius in kilometers based on the background values
+    Assigns the radius in kilometers based on the background values
 
     Parameters
     ----------
-    refl_bkg : array
-        array of background reflectivity values
-    val_for_max_conv_rad : float
-        reflectivity value for maximum convective radius (5 km)
-    max_conv_rad : float, optional
-        maximum convective radius in kilometers
+    field_bkg : array
+        array of background field values
+    val_for_max_rad : float
+        field value for maximum radius (5 km)
+    max_rad : float, optional
+        maximum radius in kilometers
 
     Returns
     -------
-    convRadiuskm : array
-        array of convective radii based on background values and val for max. conv radius
+    radius_array_km : array
+        array of radii based on background values and val for max. radius
     """
 
-    convRadiuskm = np.ones_like(refl_bkg)
+    radius_array_km = np.ones_like(field_bkg)
 
-    convRadiuskm[refl_bkg >= (val_for_max_conv_rad - 15)] = max_conv_rad - 3
-    convRadiuskm[refl_bkg >= (val_for_max_conv_rad - 10)] = max_conv_rad - 2
-    convRadiuskm[refl_bkg >= (val_for_max_conv_rad - 5)] = max_conv_rad - 1
-    convRadiuskm[refl_bkg >= val_for_max_conv_rad] = max_conv_rad
+    radius_array_km[field_bkg >= (val_for_max_rad - 15)] = max_rad - 3
+    radius_array_km[field_bkg >= (val_for_max_rad - 10)] = max_rad - 2
+    radius_array_km[field_bkg >= (val_for_max_rad - 5)] = max_rad - 1
+    radius_array_km[field_bkg >= val_for_max_rad] = max_rad
 
-    return convRadiuskm
+    return radius_array_km
 
 
-def classify_conv_strat_array(
-    refl,
-    conv_strat_array,
-    conv_core_array,
+def classify_feature_array(
+    field,
+    feature_array,
+    core_array,
     NOSFCECHO,
-    CONV,
-    SF,
+    FEAT_VAL,
+    BKGD_VAL,
     WEAKECHO,
-    CS_CORE,
+    CORE,
     MINDBZUSE,
     WEAKECHOTHRES,
 ):
     """
-    Does and initial convective stratiform classification
+    Does an initial feature detection
 
     Parameters
     ----------
-    refl : array
-        Array of reflectivity values
-    conv_strat_array : array
-        Array with convective stratiform classifications
-    conv_core_array : array
-        Array with convective cores
+    field : array
+        Array of field values
+    feature_array : array
+        Array with feature detection
+    core_array : array
+        Array with cores
     NOSFCECHO : int
         Value to assign points classified as no surface echo
-    CONV : int
-        Value to assign points classified as convective
-    SF : int
-        Value to assign points classified as stratiform
+    FEAT_VAL : int
+        Value to assign points classified as features
+    BKGD_VAL : int
+        Value to assign points classified as background echo
     WEAKECHO : int
         Value to assign points classfied as weak echo
-    CS_CORE : int
-        Value assigned to convective cores in conv_core_array
+    CORE : int
+        Value assigned to cores in core_array
     MINDBZUSE : float
         Minimum dBZ value to consider in classification, all values below this will be set to NOSFCECHO
     WEAKECHOTHRES : float
@@ -796,20 +771,20 @@ def classify_conv_strat_array(
 
     Returns
     -------
-    conv_strat_array : array
+    feature_array : array
         array of classifications
     """
 
     # assuming order so that each point is only assigned one time, no overlapping assignment
     # initially, assign every point to stratiform
-    conv_strat_array[:] = SF
-    # where reflectivity is masked, set to no surface echo
-    conv_strat_array[refl.mask] = NOSFCECHO
-    # assign convective cores to CONV
-    conv_strat_array[conv_core_array == CS_CORE] = CONV
-    # assign reflectivity less than weakechothres to weak echo
-    conv_strat_array[refl < WEAKECHOTHRES] = WEAKECHO
-    # assign reflectivity less than minimum to no surface echo
-    conv_strat_array[refl < MINDBZUSE] = NOSFCECHO
-
-    return conv_strat_array
+    feature_array[:] = BKGD_VAL
+    # where field is masked, set to no surface echo
+    feature_array[field.mask] = NOSFCECHO
+    # assign cores to FEAT
+    feature_array[core_array == CORE] = FEAT_VAL
+    # assign field values less than weakechothres to weak echo
+    feature_array[field < WEAKECHOTHRES] = WEAKECHO
+    # assign field values less than minimum to no surface echo
+    feature_array[field < MINDBZUSE] = NOSFCECHO
+
+    return feature_array
diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index f118a2e99a..469bef6e4f 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -8,7 +8,7 @@
 import numpy as np
 
 from ..config import get_field_name, get_fillvalue, get_metadata
-from ._echo_class import _revised_conv_strat, steiner_class_buff
+from ._echo_class import _feature_detection, steiner_class_buff
 
 
 def steiner_conv_strat(
@@ -125,7 +125,7 @@ def steiner_conv_strat(
     }
 
 
-def conv_strat_yuter(
+def feature_detection(
     grid,
     dx=None,
     dy=None,
@@ -139,24 +139,25 @@ def conv_strat_yuter(
     use_addition=True,
     calc_thres=0.75,
     weak_echo_thres=5.0,
-    min_dBZ_used=5.0,
+    min_val_used=5.0,
     dB_averaging=True,
     remove_small_objects=True,
     min_km2_size=10,
-    val_for_max_conv_rad=30,
-    max_conv_rad_km=5.0,
-    cs_core=3,
+    val_for_max_rad=30,
+    max_rad_km=5.0,
+    core=3,
     nosfcecho=0,
     weakecho=3,
-    sf=1,
-    conv=2,
-    refl_field=None,
+    bkgd_val=1,
+    feat_val=2,
+    field=None,
     estimate_flag=True,
     estimate_offset=5,
 ):
     """
-    Partition reflectivity into convective-stratiform using the Yuter et al. (2005)
-    and Yuter and Houze (1997) algorithm.
+    This function can be used to detect features in a field (e.g. reflectivity, rain rate, snow rate,
+    etc.) described by Tomkins et al. (2023) and based on original convective-stratiform algorithms developed by
+    Steiner et al. (1995), Yuter et al. (2005) and Yuter and Houze (1997) algorithm.
 
     Parameters
     ----------
@@ -166,18 +167,17 @@ def conv_strat_yuter(
         The x- and y-dimension resolutions in meters, respectively. If None
         the resolution is determined from the first two axes values parsed from grid object.
     level_m : float, optional
-        Desired height in meters to classify with convective stratiform algorithm.
+        Desired height in meters to run feature detection algorithm.
     always_core_thres : float, optional
-        Threshold for points that are always convective. All values above the threshold are classifed as convective
-        See Yuter et al. (2005) for more detail.
+        Threshold for points that are always features. All values above the threshold are classified as features.
     bkg_rad_km : float, optional
         Radius to compute background reflectivity in kilometers. Default is 11 km. Recommended to be at least 3 x
         grid spacing
     use_cosine : bool, optional
         Boolean used to determine if a cosine scheme (see Yuter and Houze (1997)) should be used for identifying
-        convective cores (True) or if a simpler scalar scheme (False) should be used.
+        cores (True) or if a simpler scalar scheme (False) should be used.
     max_diff : float, optional
-        Maximum difference between background average and reflectivity in order to be classified as convective.
+        Maximum difference between background average and reflectivity in order to be classified as features.
         "a" value in Eqn. B1 in Yuter and Houze (1997)
     zero_diff_cos_val : float, optional
         Value where difference between background average and reflectivity is zero in the cosine function
@@ -190,50 +190,51 @@ def conv_strat_yuter(
         Minimum percentage of points needed to be considered in background average calculation
     weak_echo_thres : float, optional
         Threshold for determining weak echo. All values below this threshold will be considered weak echo
-    min_dBZ_used : float, optional
-        Minimum dBZ value used for classification. All values below this threshold will be considered no surface echo
-        See Yuter and Houze (1997) and Yuter et al. (2005) for more detail.
+    min_val_used : float, optional
+        Minimum value used for classification. All values below this threshold will be considered no surface echo
+        See Yuter and Houze (1997) and Yuter et al. (2005) for more detail. Units based on input field
     dB_averaging : bool, optional
         True if using dBZ reflectivity values that need to be converted to linear Z before averaging. False for
         other non-dBZ values (i.e. snow rate)
     remove_small_objects : bool, optional
-        Determines if small objects should be removed from convective core array. Default is True.
+        Determines if small objects should be removed from core array. Default is True.
     min_km2_size : float, optional
-        Minimum size of convective cores to be considered. Cores less than this size will be removed. Default is 10
-        km^2.
-    val_for_max_conv_rad : float, optional
-        dBZ for maximum convective radius. Convective cores with values above this will have the maximum convective
-        radius
-    max_conv_rad_km : float, optional
-        Maximum radius around convective cores to classify as convective. Default is 5 km
-    cs_core : int, optional
-        Value for points classified as convective cores
+        Minimum size of Cores to be considered. Cores less than this size will be removed. Default is 10 km^2.
+    val_for_max_rad : float, optional
+        value used for maximum radius. Cores with values above this will have the maximum radius incorporated.
+    max_rad_km : float, optional
+        Maximum radius around cores to classify as feature. Default is 5 km
+    core : int, optional
+        Value for points classified as cores
     nosfcecho : int, optional
-        Value for points classified as no surface echo, based on min_dBZ_used
+        Value for points classified as no surface echo, based on min_val_used
     weakecho : int, optional
-        Value for points classified as weak echo, based on weak_echo_thres
-    sf : int, optional
-        Value for points classified as stratiform
-    conv : int, optional
-        Value for points classified as convective
-    refl_field : str, optional
-        Field in grid to use as the reflectivity during partitioning. None will use the default reflectivity
-        field name from the Py-ART configuration file.
+        Value for points classified as weak echo, based on weak_echo_thres.
+    bkgd_val : int, optional
+        Value for points classified as background echo.
+    feat_val : int, optional
+        Value for points classified as features.
+    field : str, optional
+        Field in grid to find objects in. None will use the default reflectivity field name from the Py-ART
+        configuration file.
     estimate_flag : bool, optional
         Determines if over/underestimation should be applied. If true, the algorithm will also be run on the same field
         wih the estimate_offset added and the same field with the estimate_offset subtracted.
         Default is True (recommended)
     estimate_offset : float, optional
-        Value used to offset the reflectivity values by for the over/underestimation application. Default value is 5
-        dBZ.
+        Value used to offset the field values by for the over/underestimation application. Default value is 5 dBZ.
 
     Returns
     -------
-    convsf_dict : dict
-        Convective-stratiform classification dictionary.
+    feature_dict : dict
+        Feature detection classification dictionary.
 
     References
     ----------
+    Steiner, M. R., R. A. Houze Jr., and S. E. Yuter, 1995: Climatological
+    Characterization of Three-Dimensional Storm Structure from Operational
+    Radar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007.
+
     Yuter, S. E., and R. A. Houze, Jr., 1997: Measurements of raindrop size
     distributions over the Pacific warm pool and implications for Z-R relations.
     J. Appl. Meteor., 36, 847-867.
@@ -245,14 +246,14 @@ def conv_strat_yuter(
 
     """
 
-    # Maxmimum convective radius must be less than 5 km
-    if max_conv_rad_km > 5:
-        print("Max conv radius must be less than 5 km, exiting")
+    # Maxmimum radius must be less than 5 km
+    if max_rad_km > 5:
+        print("Max radius must be less than 5 km, exiting")
         raise
 
     # Parse field parameters
-    if refl_field is None:
-        refl_field = get_field_name("reflectivity")
+    if field is None:
+        field = get_field_name("reflectivity")
         dB_averaging = True
 
     # parse dx and dy if None
@@ -274,15 +275,15 @@ def conv_strat_yuter(
     # Get reflectivity data at desired level
     if level_m is None:
         try:
-            ze = np.ma.copy(grid.fields[refl_field]["data"][0, :, :])
+            ze = np.ma.copy(grid.fields[field]["data"][0, :, :])
         except:
-            ze = np.ma.copy(grid.fields[refl_field]["data"][:, :])
+            ze = np.ma.copy(grid.fields[field]["data"][:, :])
     else:
         zslice = np.argmin(np.abs(z - level_m))
-        ze = np.ma.copy(grid.fields[refl_field]["data"][zslice, :, :])
+        ze = np.ma.copy(grid.fields[field]["data"][zslice, :, :])
 
-    # run convective stratiform algorithm
-    _, _, convsf_best = _revised_conv_strat(
+    # run feature detection algorithm
+    _, _, feature_best = _feature_detection(
         ze,
         dx,
         dy,
@@ -295,42 +296,37 @@ def conv_strat_yuter(
         use_addition=use_addition,
         calc_thres=calc_thres,
         weak_echo_thres=weak_echo_thres,
-        min_dBZ_used=min_dBZ_used,
+        min_val_used=min_val_used,
         dB_averaging=dB_averaging,
         remove_small_objects=remove_small_objects,
         min_km2_size=min_km2_size,
-        val_for_max_conv_rad=val_for_max_conv_rad,
-        max_conv_rad_km=max_conv_rad_km,
-        cs_core=cs_core,
+        val_for_max_rad=val_for_max_rad,
+        max_rad_km=max_rad_km,
+        core=core,
         nosfcecho=nosfcecho,
         weakecho=weakecho,
-        sf=sf,
-        conv=conv,
+        bkgd_val=bkgd_val,
+        feat_val=feat_val,
     )
 
     # put data into a dictionary to be added as a field
-    convsf_dict = {
-        "convsf": {
-            "data": convsf_best,
-            "standard_name": "convsf",
-            "long_name": "Convective stratiform classification",
+    feature_dict = {
+        "feature_detection": {
+            "data": feature_best,
+            "standard_name": "feature_detection",
+            "long_name": "Feature Detection",
             "valid_min": 0,
             "valid_max": 3,
             "comment_1": (
-                "Convective-stratiform echo "
-                "classification based on "
-                "Yuter and Houze (1997) and Yuter et al. (2005)"
-            ),
-            "comment_2": (
-                "0 = No surface echo/Undefined, 1 = Stratiform, "
-                "2 = Convective, 3 = weak echo"
+                "{0} = No surface echo/Undefined, {1} = Background echo, {2} = Features, {3} = weak echo".format(
+                    nosfcecho, bkgd_val, feat_val, weakecho)
             ),
         }
     }
 
     # If estimation is True, run the algorithm on the field with offset subtracted and the field with the offset added
     if estimate_flag:
-        _, _, convsf_under = _revised_conv_strat(
+        _, _, feature_under = _feature_detection(
             ze - estimate_offset,
             dx,
             dy,
@@ -343,20 +339,20 @@ def conv_strat_yuter(
             use_addition=use_addition,
             calc_thres=calc_thres,
             weak_echo_thres=weak_echo_thres,
-            min_dBZ_used=min_dBZ_used,
+            min_val_used=min_val_used,
             dB_averaging=dB_averaging,
             remove_small_objects=remove_small_objects,
             min_km2_size=min_km2_size,
-            val_for_max_conv_rad=val_for_max_conv_rad,
-            max_conv_rad_km=max_conv_rad_km,
-            cs_core=cs_core,
+            val_for_max_rad=val_for_max_rad,
+            max_rad_km=max_rad_km,
+            core=core,
             nosfcecho=nosfcecho,
             weakecho=weakecho,
-            sf=sf,
-            conv=conv,
+            bkgd_val=bkgd_val,
+            feat_val=feat_val,
         )
 
-        _, _, convsf_over = _revised_conv_strat(
+        _, _, feature_over = _feature_detection(
             ze + estimate_offset,
             dx,
             dy,
@@ -369,53 +365,45 @@ def conv_strat_yuter(
             use_addition=use_addition,
             calc_thres=calc_thres,
             weak_echo_thres=weak_echo_thres,
-            min_dBZ_used=min_dBZ_used,
+            min_val_used=min_val_used,
             dB_averaging=dB_averaging,
             remove_small_objects=remove_small_objects,
             min_km2_size=min_km2_size,
-            val_for_max_conv_rad=val_for_max_conv_rad,
-            max_conv_rad_km=max_conv_rad_km,
-            cs_core=cs_core,
+            val_for_max_rad=val_for_max_rad,
+            max_rad_km=max_rad_km,
+            core=core,
             nosfcecho=nosfcecho,
             weakecho=weakecho,
-            sf=sf,
-            conv=conv,
+            bkgd_val=bkgd_val,
+            feat_val=feat_val,
         )
 
         # save into dictionaries
-        convsf_dict["convsf_under"] = {
-            "data": convsf_under,
-            "standard_name": "convsf_under",
-            "long_name": "Convective stratiform classification (Underestimate)",
+        feature_dict["feature_under"] = {
+            "data": feature_under,
+            "standard_name": "feature_detection_under",
+            "long_name": "Feature Detection Underestimate",
             "valid_min": 0,
             "valid_max": 3,
             "comment_1": (
-                "Convective-stratiform echo "
-                "classification based on "
-                "Yuter and Houze (1997) and Yuter et al. (2005)"
-            ),
-            "comment_2": (
-                "0 = Undefined, 1 = Stratiform, " "2 = Convective, 3 = weak echo"
+                "{0} = No surface echo/Undefined, {1} = Background echo, {2} = Features, {3} = weak echo".format(
+                    nosfcecho, bkgd_val, feat_val, weakecho)
             ),
         }
 
-        convsf_dict["convsf_over"] = {
-            "data": convsf_over,
-            "standard_name": "convsf_under",
-            "long_name": "Convective stratiform classification (Overestimate)",
+        feature_dict["feature_over"] = {
+            "data": feature_over,
+            "standard_name": "feature_detection_over",
+            "long_name": "Feature Detection Overestimate",
             "valid_min": 0,
             "valid_max": 3,
             "comment_1": (
-                "Convective-stratiform echo "
-                "classification based on "
-                "Yuter and Houze (1997)"
-            ),
-            "comment_2": (
-                "0 = Undefined, 1 = Stratiform, " "2 = Convective, 3 = weak echo"
+                "{0} = No surface echo/Undefined, {1} = Background echo, {2} = Features, {3} = weak echo".format(
+                    nosfcecho, bkgd_val, feat_val, weakecho)
             ),
         }
 
-    return convsf_dict
+    return feature_dict
 
 
 def hydroclass_semisupervised(

From 09397917b04d1f8b5a6fe62cd0878ff1085a0941 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Wed, 4 Oct 2023 14:09:52 -0400
Subject: [PATCH 02/21] Add new reference

---
 pyart/retrieve/echo_class.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index 469bef6e4f..09b5cfeb25 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -244,6 +244,9 @@ def feature_detection(
     2005: Physical characterization of tropical oceanic convection observed in
     KWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1
 
+    Tomkins, L. M., S. E. Yuter, and M. A. Miller, 2023: Objective identification
+    of faint and strong snow bands in radar observations of winter storms. in prep.
+
     """
 
     # Maxmimum radius must be less than 5 km

From 79c598e2d79ee615910fd0452191d219459758ee Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Wed, 4 Oct 2023 14:25:12 -0400
Subject: [PATCH 03/21] Incorporate binary closing option to functions

---
 pyart/retrieve/_echo_class.py |  22 +++++--
 pyart/retrieve/echo_class.py  | 114 +++++++++++-----------------------
 2 files changed, 52 insertions(+), 84 deletions(-)

diff --git a/pyart/retrieve/_echo_class.py b/pyart/retrieve/_echo_class.py
index b3dadce2a3..24daf851f8 100644
--- a/pyart/retrieve/_echo_class.py
+++ b/pyart/retrieve/_echo_class.py
@@ -256,9 +256,10 @@ def _feature_detection(
     dB_averaging=True,
     remove_small_objects=True,
     min_km2_size=10,
+    binary_close=False,
     val_for_max_rad=30,
     max_rad_km=5.0,
-    core=3,
+    core_val=3,
     nosfcecho=0,
     weakecho=3,
     bkgd_val=1,
@@ -311,11 +312,13 @@ def _feature_detection(
         Determines if small objects should be removed from core array. Default is True.
     min_km2_size : float, optional
         Minimum size of Cores to be considered. Cores less than this size will be removed. Default is 10 km^2.
+    binary_close : bool, optional
+        Determines if a binary closing should be performed on the cores. Default is False.
     val_for_max_rad : float, optional
         value used for maximum radius. Cores with values above this will have the maximum radius incorporated.
     max_rad_km : float, optional
         Maximum radius around cores to classify as feature. Default is 5 km
-    core : int, optional
+    core_val : int, optional
         Value for points classified as cores
     nosfcecho : int, optional
         Value for points classified as no surface echo, based on min_val_used
@@ -375,9 +378,9 @@ def _feature_detection(
 
     # Get core array from cosine scheme, or scalar scheme
     if use_cosine:
-        core_array = core_cos_scheme(field, field_bkg, max_diff, zero_diff_cos_val, always_core_thres, core)
+        core_array = core_cos_scheme(field, field_bkg, max_diff, zero_diff_cos_val, always_core_thres, core_val)
     else:
-        core_array = core_scalar_scheme(field, field_bkg, scalar_diff, always_core_thres, core,
+        core_array = core_scalar_scheme(field, field_bkg, scalar_diff, always_core_thres, core_val,
                                         use_addition=use_addition)
 
     # Assign radii based on background field
@@ -400,6 +403,13 @@ def _feature_detection(
             if size_lab < min_pix_size:
                 core_array[cc_labels == lab] = 0
 
+    # perform binary closing
+    if binary_close:
+        # binary closing - returns binary array
+        close_core = scipy.ndimage.binary_closing(core_array).astype(int)
+        # set values to core values
+        core_array = close_core * core_val
+
     # Incorporate radius using binary dilation
     # Create empty array for assignment
     temp_assignment = np.zeros_like(core_array)
@@ -423,12 +433,12 @@ def _feature_detection(
 
     # add dilated cores to original array
     core_copy = np.ma.copy(core_array)
-    core_copy[temp_assignment >= 1] = core
+    core_copy[temp_assignment >= 1] = core_val
 
     # Now do feature detection
     feature_array = np.zeros_like(field)
     feature_array = classify_feature_array(field, feature_array, core_copy, nosfcecho, feat_val, bkgd_val, weakecho,
-                                           core, min_val_used, weak_echo_thres)
+                                           core_val, min_val_used, weak_echo_thres)
     # mask where field is masked
     feature_array = np.ma.masked_where(field.mask, feature_array)
 
diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index 09b5cfeb25..ba125b24d9 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -143,9 +143,10 @@ def feature_detection(
     dB_averaging=True,
     remove_small_objects=True,
     min_km2_size=10,
+    binary_close=False,
     val_for_max_rad=30,
     max_rad_km=5.0,
-    core=3,
+    core_val=3,
     nosfcecho=0,
     weakecho=3,
     bkgd_val=1,
@@ -200,11 +201,13 @@ def feature_detection(
         Determines if small objects should be removed from core array. Default is True.
     min_km2_size : float, optional
         Minimum size of Cores to be considered. Cores less than this size will be removed. Default is 10 km^2.
+    binary_close : bool, optional
+        Determines if a binary closing should be performed on the cores. Default is False.
     val_for_max_rad : float, optional
         value used for maximum radius. Cores with values above this will have the maximum radius incorporated.
     max_rad_km : float, optional
         Maximum radius around cores to classify as feature. Default is 5 km
-    core : int, optional
+    core_val : int, optional
         Value for points classified as cores
     nosfcecho : int, optional
         Value for points classified as no surface echo, based on min_val_used
@@ -286,31 +289,16 @@ def feature_detection(
         ze = np.ma.copy(grid.fields[field]["data"][zslice, :, :])
 
     # run feature detection algorithm
-    _, _, feature_best = _feature_detection(
-        ze,
-        dx,
-        dy,
-        always_core_thres=always_core_thres,
-        bkg_rad_km=bkg_rad_km,
-        use_cosine=use_cosine,
-        max_diff=max_diff,
-        zero_diff_cos_val=zero_diff_cos_val,
-        scalar_diff=scalar_diff,
-        use_addition=use_addition,
-        calc_thres=calc_thres,
-        weak_echo_thres=weak_echo_thres,
-        min_val_used=min_val_used,
-        dB_averaging=dB_averaging,
-        remove_small_objects=remove_small_objects,
-        min_km2_size=min_km2_size,
-        val_for_max_rad=val_for_max_rad,
-        max_rad_km=max_rad_km,
-        core=core,
-        nosfcecho=nosfcecho,
-        weakecho=weakecho,
-        bkgd_val=bkgd_val,
-        feat_val=feat_val,
-    )
+    _, _, feature_best = _feature_detection(ze, dx, dy, always_core_thres=always_core_thres, bkg_rad_km=bkg_rad_km,
+                                            use_cosine=use_cosine, max_diff=max_diff,
+                                            zero_diff_cos_val=zero_diff_cos_val, scalar_diff=scalar_diff,
+                                            use_addition=use_addition, calc_thres=calc_thres,
+                                            weak_echo_thres=weak_echo_thres, min_val_used=min_val_used,
+                                            dB_averaging=dB_averaging, remove_small_objects=remove_small_objects,
+                                            min_km2_size=min_km2_size, binary_close=binary_close,
+                                            val_for_max_rad=val_for_max_rad, max_rad_km=max_rad_km, core_val=core_val,
+                                            nosfcecho=nosfcecho, weakecho=weakecho, bkgd_val=bkgd_val,
+                                            feat_val=feat_val)
 
     # put data into a dictionary to be added as a field
     feature_dict = {
@@ -329,57 +317,27 @@ def feature_detection(
 
     # If estimation is True, run the algorithm on the field with offset subtracted and the field with the offset added
     if estimate_flag:
-        _, _, feature_under = _feature_detection(
-            ze - estimate_offset,
-            dx,
-            dy,
-            always_core_thres=always_core_thres,
-            bkg_rad_km=bkg_rad_km,
-            use_cosine=use_cosine,
-            max_diff=max_diff,
-            zero_diff_cos_val=zero_diff_cos_val,
-            scalar_diff=scalar_diff,
-            use_addition=use_addition,
-            calc_thres=calc_thres,
-            weak_echo_thres=weak_echo_thres,
-            min_val_used=min_val_used,
-            dB_averaging=dB_averaging,
-            remove_small_objects=remove_small_objects,
-            min_km2_size=min_km2_size,
-            val_for_max_rad=val_for_max_rad,
-            max_rad_km=max_rad_km,
-            core=core,
-            nosfcecho=nosfcecho,
-            weakecho=weakecho,
-            bkgd_val=bkgd_val,
-            feat_val=feat_val,
-        )
-
-        _, _, feature_over = _feature_detection(
-            ze + estimate_offset,
-            dx,
-            dy,
-            always_core_thres=always_core_thres,
-            bkg_rad_km=bkg_rad_km,
-            use_cosine=use_cosine,
-            max_diff=max_diff,
-            zero_diff_cos_val=zero_diff_cos_val,
-            scalar_diff=scalar_diff,
-            use_addition=use_addition,
-            calc_thres=calc_thres,
-            weak_echo_thres=weak_echo_thres,
-            min_val_used=min_val_used,
-            dB_averaging=dB_averaging,
-            remove_small_objects=remove_small_objects,
-            min_km2_size=min_km2_size,
-            val_for_max_rad=val_for_max_rad,
-            max_rad_km=max_rad_km,
-            core=core,
-            nosfcecho=nosfcecho,
-            weakecho=weakecho,
-            bkgd_val=bkgd_val,
-            feat_val=feat_val,
-        )
+        _, _, feature_under = _feature_detection(ze - estimate_offset, dx, dy, always_core_thres=always_core_thres,
+                                                 bkg_rad_km=bkg_rad_km, use_cosine=use_cosine, max_diff=max_diff,
+                                                 zero_diff_cos_val=zero_diff_cos_val, scalar_diff=scalar_diff,
+                                                 use_addition=use_addition, calc_thres=calc_thres,
+                                                 weak_echo_thres=weak_echo_thres, min_val_used=min_val_used,
+                                                 dB_averaging=dB_averaging, remove_small_objects=remove_small_objects,
+                                                 min_km2_size=min_km2_size, binary_close=binary_close,
+                                                 val_for_max_rad=val_for_max_rad, max_rad_km=max_rad_km, core_val=core_val,
+                                                 nosfcecho=nosfcecho, weakecho=weakecho, bkgd_val=bkgd_val,
+                                                 feat_val=feat_val)
+
+        _, _, feature_over = _feature_detection(ze + estimate_offset, dx, dy, always_core_thres=always_core_thres,
+                                                bkg_rad_km=bkg_rad_km, use_cosine=use_cosine, max_diff=max_diff,
+                                                zero_diff_cos_val=zero_diff_cos_val, scalar_diff=scalar_diff,
+                                                use_addition=use_addition, calc_thres=calc_thres,
+                                                weak_echo_thres=weak_echo_thres, min_val_used=min_val_used,
+                                                dB_averaging=dB_averaging, remove_small_objects=remove_small_objects,
+                                                min_km2_size=min_km2_size, binary_close=binary_close,
+                                                val_for_max_rad=val_for_max_rad, max_rad_km=max_rad_km, core_val=core_val,
+                                                nosfcecho=nosfcecho, weakecho=weakecho, bkgd_val=bkgd_val,
+                                                feat_val=feat_val)
 
         # save into dictionaries
         feature_dict["feature_under"] = {

From 6b23cc6b23807b2f380457691b0aeec09dc1903f Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Wed, 4 Oct 2023 17:01:10 -0400
Subject: [PATCH 04/21] Add author to functions

---
 pyart/retrieve/echo_class.py | 8 +++++---
 pyart/util/radar_utils.py    | 2 ++
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index ba125b24d9..192502b4e6 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -160,6 +160,8 @@ def feature_detection(
     etc.) described by Tomkins et al. (2023) and based on original convective-stratiform algorithms developed by
     Steiner et al. (1995), Yuter et al. (2005) and Yuter and Houze (1997) algorithm.
 
+    Author: Laura Tomkins (@lauratomkins)
+
     Parameters
     ----------
     grid : Grid
@@ -303,7 +305,7 @@ def feature_detection(
     # put data into a dictionary to be added as a field
     feature_dict = {
         "feature_detection": {
-            "data": feature_best,
+            "data": feature_best[None,:,:],
             "standard_name": "feature_detection",
             "long_name": "Feature Detection",
             "valid_min": 0,
@@ -341,7 +343,7 @@ def feature_detection(
 
         # save into dictionaries
         feature_dict["feature_under"] = {
-            "data": feature_under,
+            "data": feature_under[None,:,:],
             "standard_name": "feature_detection_under",
             "long_name": "Feature Detection Underestimate",
             "valid_min": 0,
@@ -353,7 +355,7 @@ def feature_detection(
         }
 
         feature_dict["feature_over"] = {
-            "data": feature_over,
+            "data": feature_over[None,:,:],
             "standard_name": "feature_detection_over",
             "long_name": "Feature Detection Overestimate",
             "valid_min": 0,
diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py
index a2eacafc5a..1a0246f1d7 100644
--- a/pyart/util/radar_utils.py
+++ b/pyart/util/radar_utils.py
@@ -499,6 +499,8 @@ def image_mute_radar(radar, field, mute_field, mute_threshold, field_threshold=N
     the correlation coefficient is less than a certain threshold to discern
     melting precipitation.
 
+    Author: Laura Tomkins (@lauratomkins)
+
     Parameters
     ----------
     radar : Radar

From f540d042c16069b02350fa6e57745058a851fc29 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Thu, 5 Oct 2023 15:58:15 -0400
Subject: [PATCH 05/21] Slight reformat for linting

---
 pyart/retrieve/_echo_class.py |  31 ++++++---
 pyart/retrieve/echo_class.py  | 122 +++++++++++++++++++++++-----------
 2 files changed, 107 insertions(+), 46 deletions(-)

diff --git a/pyart/retrieve/_echo_class.py b/pyart/retrieve/_echo_class.py
index 24daf851f8..38f8e32a37 100644
--- a/pyart/retrieve/_echo_class.py
+++ b/pyart/retrieve/_echo_class.py
@@ -378,13 +378,18 @@ def _feature_detection(
 
     # Get core array from cosine scheme, or scalar scheme
     if use_cosine:
-        core_array = core_cos_scheme(field, field_bkg, max_diff, zero_diff_cos_val, always_core_thres, core_val)
+        core_array = core_cos_scheme(
+            field, field_bkg, max_diff, zero_diff_cos_val, always_core_thres, core_val
+        )
     else:
-        core_array = core_scalar_scheme(field, field_bkg, scalar_diff, always_core_thres, core_val,
-                                        use_addition=use_addition)
+        core_array = core_scalar_scheme(
+            field, field_bkg, scalar_diff, always_core_thres, core_val, use_addition=use_addition
+        )
 
     # Assign radii based on background field
-    radius_array_km = assign_feature_radius_km(field_bkg, val_for_max_rad=val_for_max_rad, max_rad=max_rad_km)
+    radius_array_km = assign_feature_radius_km(
+        field_bkg, val_for_max_rad=val_for_max_rad, max_rad=max_rad_km
+    )
 
     # remove small objects in core array
     if remove_small_objects:
@@ -437,8 +442,18 @@ def _feature_detection(
 
     # Now do feature detection
     feature_array = np.zeros_like(field)
-    feature_array = classify_feature_array(field, feature_array, core_copy, nosfcecho, feat_val, bkgd_val, weakecho,
-                                           core_val, min_val_used, weak_echo_thres)
+    feature_array = classify_feature_array(
+        field,
+        feature_array,
+        core_copy,
+        nosfcecho,
+        feat_val,
+        bkgd_val,
+        weakecho,
+        core_val,
+        min_val_used,
+        weak_echo_thres
+    )
     # mask where field is masked
     feature_array = np.ma.masked_where(field.mask, feature_array)
 
@@ -674,9 +689,7 @@ def core_scalar_scheme(
     return core_array
 
 
-def create_radius_mask(
-    max_diameter, radius_km, x_spacing, y_spacing, center_mask_x
-):
+def create_radius_mask(max_diameter, radius_km, x_spacing, y_spacing, center_mask_x):
     """
     Creates a circular mask based on input diameter
 
diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index 192502b4e6..e46da1ca35 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -291,27 +291,43 @@ def feature_detection(
         ze = np.ma.copy(grid.fields[field]["data"][zslice, :, :])
 
     # run feature detection algorithm
-    _, _, feature_best = _feature_detection(ze, dx, dy, always_core_thres=always_core_thres, bkg_rad_km=bkg_rad_km,
-                                            use_cosine=use_cosine, max_diff=max_diff,
-                                            zero_diff_cos_val=zero_diff_cos_val, scalar_diff=scalar_diff,
-                                            use_addition=use_addition, calc_thres=calc_thres,
-                                            weak_echo_thres=weak_echo_thres, min_val_used=min_val_used,
-                                            dB_averaging=dB_averaging, remove_small_objects=remove_small_objects,
-                                            min_km2_size=min_km2_size, binary_close=binary_close,
-                                            val_for_max_rad=val_for_max_rad, max_rad_km=max_rad_km, core_val=core_val,
-                                            nosfcecho=nosfcecho, weakecho=weakecho, bkgd_val=bkgd_val,
-                                            feat_val=feat_val)
+    _, _, feature_best = _feature_detection(
+        ze,
+        dx,
+        dy,
+        always_core_thres=always_core_thres,
+        bkg_rad_km=bkg_rad_km,
+        use_cosine=use_cosine,
+        max_diff=max_diff,
+        zero_diff_cos_val=zero_diff_cos_val,
+        scalar_diff=scalar_diff,
+        use_addition=use_addition,
+        calc_thres=calc_thres,
+        weak_echo_thres=weak_echo_thres,
+        min_val_used=min_val_used,
+        dB_averaging=dB_averaging,
+        remove_small_objects=remove_small_objects,
+        min_km2_size=min_km2_size,
+        binary_close=binary_close,
+        val_for_max_rad=val_for_max_rad,
+        max_rad_km=max_rad_km,
+        core_val=core_val,
+        nosfcecho=nosfcecho,
+        weakecho=weakecho,
+        bkgd_val=bkgd_val,
+        feat_val=feat_val
+    )
 
     # put data into a dictionary to be added as a field
     feature_dict = {
         "feature_detection": {
-            "data": feature_best[None,:,:],
+            "data": feature_best[None, :, :],
             "standard_name": "feature_detection",
             "long_name": "Feature Detection",
             "valid_min": 0,
             "valid_max": 3,
             "comment_1": (
-                "{0} = No surface echo/Undefined, {1} = Background echo, {2} = Features, {3} = weak echo".format(
+                "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format(
                     nosfcecho, bkgd_val, feat_val, weakecho)
             ),
         }
@@ -319,49 +335,81 @@ def feature_detection(
 
     # If estimation is True, run the algorithm on the field with offset subtracted and the field with the offset added
     if estimate_flag:
-        _, _, feature_under = _feature_detection(ze - estimate_offset, dx, dy, always_core_thres=always_core_thres,
-                                                 bkg_rad_km=bkg_rad_km, use_cosine=use_cosine, max_diff=max_diff,
-                                                 zero_diff_cos_val=zero_diff_cos_val, scalar_diff=scalar_diff,
-                                                 use_addition=use_addition, calc_thres=calc_thres,
-                                                 weak_echo_thres=weak_echo_thres, min_val_used=min_val_used,
-                                                 dB_averaging=dB_averaging, remove_small_objects=remove_small_objects,
-                                                 min_km2_size=min_km2_size, binary_close=binary_close,
-                                                 val_for_max_rad=val_for_max_rad, max_rad_km=max_rad_km, core_val=core_val,
-                                                 nosfcecho=nosfcecho, weakecho=weakecho, bkgd_val=bkgd_val,
-                                                 feat_val=feat_val)
-
-        _, _, feature_over = _feature_detection(ze + estimate_offset, dx, dy, always_core_thres=always_core_thres,
-                                                bkg_rad_km=bkg_rad_km, use_cosine=use_cosine, max_diff=max_diff,
-                                                zero_diff_cos_val=zero_diff_cos_val, scalar_diff=scalar_diff,
-                                                use_addition=use_addition, calc_thres=calc_thres,
-                                                weak_echo_thres=weak_echo_thres, min_val_used=min_val_used,
-                                                dB_averaging=dB_averaging, remove_small_objects=remove_small_objects,
-                                                min_km2_size=min_km2_size, binary_close=binary_close,
-                                                val_for_max_rad=val_for_max_rad, max_rad_km=max_rad_km, core_val=core_val,
-                                                nosfcecho=nosfcecho, weakecho=weakecho, bkgd_val=bkgd_val,
-                                                feat_val=feat_val)
+        _, _, feature_under = _feature_detection(
+            ze - estimate_offset,
+            dx,
+            dy,
+            always_core_thres=always_core_thres,
+            bkg_rad_km=bkg_rad_km,
+            use_cosine=use_cosine,
+            max_diff=max_diff,
+            zero_diff_cos_val=zero_diff_cos_val,
+            scalar_diff=scalar_diff,
+            use_addition=use_addition,
+            calc_thres=calc_thres,
+            weak_echo_thres=weak_echo_thres,
+            min_val_used=min_val_used,
+            dB_averaging=dB_averaging,
+            remove_small_objects=remove_small_objects,
+            min_km2_size=min_km2_size,
+            binary_close=binary_close,
+            val_for_max_rad=val_for_max_rad,
+            max_rad_km=max_rad_km,
+            core_val=core_val,
+            nosfcecho=nosfcecho,
+            weakecho=weakecho,
+            bkgd_val=bkgd_val,
+            feat_val=feat_val
+        )
+
+        _, _, feature_over = _feature_detection(
+            ze + estimate_offset,
+            dx,
+            dy,
+            always_core_thres=always_core_thres,
+            bkg_rad_km=bkg_rad_km,
+            use_cosine=use_cosine,
+            max_diff=max_diff,
+            zero_diff_cos_val=zero_diff_cos_val,
+            scalar_diff=scalar_diff,
+            use_addition=use_addition,
+            calc_thres=calc_thres,
+            weak_echo_thres=weak_echo_thres,
+            min_val_used=min_val_used,
+            dB_averaging=dB_averaging,
+            remove_small_objects=remove_small_objects,
+            min_km2_size=min_km2_size,
+            binary_close=binary_close,
+            val_for_max_rad=val_for_max_rad,
+            max_rad_km=max_rad_km,
+            core_val=core_val,
+            nosfcecho=nosfcecho,
+            weakecho=weakecho,
+            bkgd_val=bkgd_val,
+            feat_val=feat_val
+        )
 
         # save into dictionaries
         feature_dict["feature_under"] = {
-            "data": feature_under[None,:,:],
+            "data": feature_under[None, :, :],
             "standard_name": "feature_detection_under",
             "long_name": "Feature Detection Underestimate",
             "valid_min": 0,
             "valid_max": 3,
             "comment_1": (
-                "{0} = No surface echo/Undefined, {1} = Background echo, {2} = Features, {3} = weak echo".format(
+                "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format(
                     nosfcecho, bkgd_val, feat_val, weakecho)
             ),
         }
 
         feature_dict["feature_over"] = {
-            "data": feature_over[None,:,:],
+            "data": feature_over[None, :, :],
             "standard_name": "feature_detection_over",
             "long_name": "Feature Detection Overestimate",
             "valid_min": 0,
             "valid_max": 3,
             "comment_1": (
-                "{0} = No surface echo/Undefined, {1} = Background echo, {2} = Features, {3} = weak echo".format(
+                "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format(
                     nosfcecho, bkgd_val, feat_val, weakecho)
             ),
         }

From f867ce40bd4308a2742076966f943bbc7bf49b2e Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Thu, 5 Oct 2023 16:01:22 -0400
Subject: [PATCH 06/21] More linting formatting

---
 pyart/retrieve/_echo_class.py |  9 +++++++--
 pyart/retrieve/echo_class.py  | 15 +++++++++------
 2 files changed, 16 insertions(+), 8 deletions(-)

diff --git a/pyart/retrieve/_echo_class.py b/pyart/retrieve/_echo_class.py
index 38f8e32a37..362d351eca 100644
--- a/pyart/retrieve/_echo_class.py
+++ b/pyart/retrieve/_echo_class.py
@@ -383,7 +383,12 @@ def _feature_detection(
         )
     else:
         core_array = core_scalar_scheme(
-            field, field_bkg, scalar_diff, always_core_thres, core_val, use_addition=use_addition
+            field,
+            field_bkg,
+            scalar_diff,
+            always_core_thres,
+            core_val,
+            use_addition=use_addition,
         )
 
     # Assign radii based on background field
@@ -452,7 +457,7 @@ def _feature_detection(
         weakecho,
         core_val,
         min_val_used,
-        weak_echo_thres
+        weak_echo_thres,
     )
     # mask where field is masked
     feature_array = np.ma.masked_where(field.mask, feature_array)
diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index e46da1ca35..afc67af9da 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -315,7 +315,7 @@ def feature_detection(
         nosfcecho=nosfcecho,
         weakecho=weakecho,
         bkgd_val=bkgd_val,
-        feat_val=feat_val
+        feat_val=feat_val,
     )
 
     # put data into a dictionary to be added as a field
@@ -328,7 +328,8 @@ def feature_detection(
             "valid_max": 3,
             "comment_1": (
                 "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format(
-                    nosfcecho, bkgd_val, feat_val, weakecho)
+                    nosfcecho, bkgd_val, feat_val, weakecho
+                )
             ),
         }
     }
@@ -359,7 +360,7 @@ def feature_detection(
             nosfcecho=nosfcecho,
             weakecho=weakecho,
             bkgd_val=bkgd_val,
-            feat_val=feat_val
+            feat_val=feat_val,
         )
 
         _, _, feature_over = _feature_detection(
@@ -386,7 +387,7 @@ def feature_detection(
             nosfcecho=nosfcecho,
             weakecho=weakecho,
             bkgd_val=bkgd_val,
-            feat_val=feat_val
+            feat_val=feat_val,
         )
 
         # save into dictionaries
@@ -398,7 +399,8 @@ def feature_detection(
             "valid_max": 3,
             "comment_1": (
                 "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format(
-                    nosfcecho, bkgd_val, feat_val, weakecho)
+                    nosfcecho, bkgd_val, feat_val, weakecho
+                )
             ),
         }
 
@@ -410,7 +412,8 @@ def feature_detection(
             "valid_max": 3,
             "comment_1": (
                 "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format(
-                    nosfcecho, bkgd_val, feat_val, weakecho)
+                    nosfcecho, bkgd_val, feat_val, weakecho
+                )
             ),
         }
 

From a778ee7592f910525dbe6eb6d71d7c2c678aedf3 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Thu, 9 Nov 2023 11:37:45 -0500
Subject: [PATCH 07/21] Add alternate option for getting under/over estimate
 field

---
 pyart/retrieve/echo_class.py | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index afc67af9da..dd524eaadd 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -154,6 +154,8 @@ def feature_detection(
     field=None,
     estimate_flag=True,
     estimate_offset=5,
+    overest_field=None,
+    underest_field=None,
 ):
     """
     This function can be used to detect features in a field (e.g. reflectivity, rain rate, snow rate,
@@ -228,6 +230,10 @@ def feature_detection(
         Default is True (recommended)
     estimate_offset : float, optional
         Value used to offset the field values by for the over/underestimation application. Default value is 5 dBZ.
+    overest_field : str, optional
+        Name of field in grid object used to calculate the overestimate if estimate_flag is True.
+    underest_field : str, optional
+        Name of field in grid object used to calculate the underestimate if estimate_flag is True.
 
     Returns
     -------
@@ -336,8 +342,16 @@ def feature_detection(
 
     # If estimation is True, run the algorithm on the field with offset subtracted and the field with the offset added
     if estimate_flag:
+
+        # get underestimate field
+        if underest_field is None:
+            under_field = ze - estimate_offset
+        elif underest_field is not None:
+            under_field = np.ma.copy(grid.fields[underest_field]["data"][0, :, :])
+
+        # run algorithm to get feature detection underestimate
         _, _, feature_under = _feature_detection(
-            ze - estimate_offset,
+            under_field,
             dx,
             dy,
             always_core_thres=always_core_thres,
@@ -363,8 +377,15 @@ def feature_detection(
             feat_val=feat_val,
         )
 
+        # get overestimate field
+        if overest_field is None:
+            over_field = ze + estimate_offset
+        elif overest_field is not None:
+            over_field = np.ma.copy(grid.fields[overest_field]["data"][0, :, :])
+
+        # run algorithm to get feature detection underestimate
         _, _, feature_over = _feature_detection(
-            ze + estimate_offset,
+            over_field,
             dx,
             dy,
             always_core_thres=always_core_thres,

From 337d8f0b57ac76511d7d64c39fda199871019994 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Thu, 9 Nov 2023 11:37:54 -0500
Subject: [PATCH 08/21] Create plot_feature_detection.py

---
 examples/retrieve/plot_feature_detection.py | 591 ++++++++++++++++++++
 1 file changed, 591 insertions(+)
 create mode 100644 examples/retrieve/plot_feature_detection.py

diff --git a/examples/retrieve/plot_feature_detection.py b/examples/retrieve/plot_feature_detection.py
new file mode 100644
index 0000000000..c7fb9ec29b
--- /dev/null
+++ b/examples/retrieve/plot_feature_detection.py
@@ -0,0 +1,591 @@
+"""
+=======================================
+Feature Detection classification
+=======================================
+This example shows how to use the feature detection algorithm. We show how to use the algorithm to detect convective
+and stratiform regions in warm season events and how the algorithm can be used to detect both weak and strong
+features in cool-season events.
+"""
+
+print(__doc__)
+
+# Author: Laura Tomkins (lmtomkin@ncsu.edu)
+# License: BSD 3 clause
+
+
+import cartopy.crs as ccrs
+import matplotlib.colors as mcolors
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pyart
+from open_radar_data import DATASETS
+
+######################################
+# How the algorithm works
+# ----------
+# The feature detection algorithm works by identifying features that exceed the background value by an amount that
+# varies with the background value. The algorithm is heavily customizable and is designed to work with a variety of
+# datasets. Here, we show several examples of how to use the algorithm with different types of radar data.
+#
+# See Steiner et al. (1995), Yuter and Houze (1997), and Yuter et al. (2005) for full details on the algorithm. A
+# manuscript (Tomkins et al. 2024) is in prep to describe feature detection in cool-season events (part 2).
+
+######################################
+# Part 1: Warm-season convective-stratiform classification
+# ----------
+# **Classification of summer convective example**
+#
+# Our first example classifies echo from a summer convective event.
+
+# read in file
+filename = pyart.testing.get_test_data("swx_20120520_0641.nc")
+radar = pyart.io.read(filename)
+
+# extract the lowest sweep
+radar = radar.extract_sweeps([0])
+
+# interpolate to grid
+grid = pyart.map.grid_from_radars(
+    (radar,),
+    grid_shape=(1, 201, 201),
+    grid_limits=((0, 10000), (-50000.0, 50000.0), (-50000.0, 50000.0)),
+    fields=["reflectivity_horizontal"],
+)
+
+# get dx dy
+dx = grid.x["data"][1] - grid.x["data"][0]
+dy = grid.y["data"][1] - grid.y["data"][0]
+
+# feature detection
+feature_dict = pyart.retrieve.feature_detection(
+    grid,
+    dx,
+    dy,
+    field="reflectivity_horizontal",
+    always_core_thres=40,
+    bkg_rad_km=20,
+    use_cosine=True,
+    max_diff=5,
+    zero_diff_cos_val=55,
+    weak_echo_thres=10,
+    max_rad_km=2,
+)
+
+# add to grid object
+# mask zero values (no surface echo)
+feature_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0)
+# mask 3 values (weak echo)
+feature_masked = np.ma.masked_equal(feature_masked, 3)
+# add dimension to array to add to grid object
+feature_dict["feature_detection"]["data"] = feature_masked
+# add field
+grid.add_field("feature_detection", feature_dict["feature_detection"], replace_existing=True)
+
+# create plot using GridMapDisplay
+# plot variables
+display = pyart.graph.GridMapDisplay(grid)
+magma_r_cmap = plt.get_cmap("magma_r")
+ref_cmap = mcolors.LinearSegmentedColormap.from_list(
+    "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))
+)
+projection = ccrs.AlbersEqualArea(
+    central_latitude=radar.latitude["data"][0],
+    central_longitude=radar.longitude["data"][0],
+)
+
+# plot
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(1, 2, 1, projection=projection)
+display.plot_grid(
+    "reflectivity_horizontal",
+    vmin=5,
+    vmax=45,
+    cmap=ref_cmap,
+    transform=ccrs.PlateCarree(),
+    ax=ax1,
+)
+ax2 = plt.subplot(1, 2, 2, projection=projection)
+display.plot_grid(
+    "feature_detection",
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
+    ax=ax2,
+    transform=ccrs.PlateCarree(),
+    ticks=[1 / 3, 1, 5 / 3],
+    ticklabs=["", "Stratiform", "Convective"],
+)
+plt.show()
+
+
+######################################
+# In addition to the default feature detection classification, the function also returns an underestimate
+# (``feature_under``) and an overestimate (``feature_over``) to take into consideration the uncertainty when choosing
+# classification parameters. The under and overestimate use the same parameters, but vary the input field by a
+# certain value (default is 5 dBZ, can be changed with ``estimate_offset``). The estimation can be turned off (
+# ``estimate_flag=False``), but we recommend keeping it turned on.
+
+# mask weak echo and no surface echo
+feature_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0)
+feature_masked = np.ma.masked_equal(feature_masked, 3)
+feature_dict["feature_detection"]["data"] = feature_masked
+# underest.
+feature_masked = np.ma.masked_equal(feature_dict["feature_under"]["data"], 0)
+feature_masked = np.ma.masked_equal(feature_masked, 3)
+feature_dict["feature_under"]["data"] = feature_masked
+# overest.
+feature_masked = np.ma.masked_equal(feature_dict["feature_over"]["data"], 0)
+feature_masked = np.ma.masked_equal(feature_masked, 3)
+feature_dict["feature_over"]["data"] = feature_masked
+
+# Plot each estimation
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(131)
+ax1.pcolormesh(
+    feature_dict["feature_detection"]["data"][0, :, :],
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
+)
+ax1.set_title("Best estimate")
+ax1.set_aspect("equal")
+ax2 = plt.subplot(132)
+ax2.pcolormesh(
+    feature_dict["feature_under"]["data"][0, :, :],
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3)
+)
+ax2.set_title("Underestimate")
+ax2.set_aspect("equal")
+ax3 = plt.subplot(133)
+ax3.pcolormesh(
+    feature_dict["feature_over"]["data"][0, :, :],
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3)
+)
+ax3.set_title("Overestimate")
+ax3.set_aspect("equal")
+plt.show()
+
+######################################
+# **Tropical example**
+#
+# Let's get a NEXRAD file from Hurricane Ian
+
+# Read in file
+nexrad_file = "s3://noaa-nexrad-level2/2022/09/28/KTBW/KTBW20220928_190142_V06"
+radar = pyart.io.read_nexrad_archive(nexrad_file)
+
+# extract the lowest sweep
+radar = radar.extract_sweeps([0])
+
+# interpolate to grid
+grid = pyart.map.grid_from_radars(
+    (radar,),
+    grid_shape=(1, 201, 201),
+    grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)),
+    fields=["reflectivity"],
+)
+
+# get dx dy
+dx = grid.x["data"][1] - grid.x["data"][0]
+dy = grid.y["data"][1] - grid.y["data"][0]
+
+# feature detection
+feature_dict = pyart.retrieve.feature_detection(
+    grid,
+    dx,
+    dy,
+    field="reflectivity",
+    always_core_thres=40,
+    bkg_rad_km=20,
+    use_cosine=True,
+    max_diff=3,
+    zero_diff_cos_val=55,
+    weak_echo_thres=5,
+    max_rad_km=2,
+    estimate_flag=False,
+)
+
+# add to grid object
+# mask zero values (no surface echo)
+feature_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0)
+# mask 3 values (weak echo)
+feature_masked = np.ma.masked_equal(feature_masked, 3)
+# add dimension to array to add to grid object
+feature_dict["feature_detection"]["data"] = feature_masked
+# add field
+grid.add_field("feature_detection", feature_dict["feature_detection"], replace_existing=True)
+
+# create plot using GridMapDisplay
+# plot variables
+display = pyart.graph.GridMapDisplay(grid)
+magma_r_cmap = plt.get_cmap("magma_r")
+ref_cmap = mcolors.LinearSegmentedColormap.from_list(
+    "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))
+)
+projection = ccrs.AlbersEqualArea(
+    central_latitude=radar.latitude["data"][0],
+    central_longitude=radar.longitude["data"][0],
+)
+# plot
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(1, 2, 1, projection=projection)
+display.plot_grid(
+    "reflectivity",
+    vmin=5,
+    vmax=45,
+    cmap=ref_cmap,
+    transform=ccrs.PlateCarree(),
+    ax=ax1,
+    axislabels_flag=False,
+)
+ax2 = plt.subplot(1, 2, 2, projection=projection)
+display.plot_grid(
+    "feature_detection",
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
+    axislabels_flag=False,
+    transform=ccrs.PlateCarree(),
+    ticks=[1 / 3, 1, 5 / 3],
+    ticklabs=["", "Stratiform", "Convective"],
+    ax=ax2,
+)
+plt.show()
+
+
+######################################
+# **Comparison with original C++ algorithm**
+#
+# To ensure that our algorithm here matches the original algorithm developed in C++, we compare the results from an
+# example during the KWAJEX field campaign. The file contains the original convective-stratiform classification from
+# the C++ algorithm. We use the algorithm implemented in Py-ART with the same settings to produce identical results.
+
+# read in file
+filename = DATASETS.fetch("convsf.19990811.221202.cdf")
+grid = pyart.io.read_grid(filename)
+
+# colormap
+convsf_cmap = mcolors.LinearSegmentedColormap.from_list('convsf', [(33 / 255, 145 / 255, 140 / 255),
+                                                                    (253 / 255, 231 / 255, 37 / 255),
+                                                                    (210 / 255, 180 / 255, 140 / 255)], N=3)
+
+# get original data from file processed in C++
+x = grid.x['data']
+y = grid.y['data']
+ref = grid.fields['maxdz']['data'][0,:,:]
+csb = grid.fields['convsf']['data'][0,:,:]
+csb_masked = np.ma.masked_less_equal(csb, 0)
+csl = grid.fields['convsf_lo']['data'][0,:,:]
+csl_masked = np.ma.masked_less_equal(csl, 0)
+csh = grid.fields['convsf_hi']['data'][0,:,:]
+csh_masked = np.ma.masked_less_equal(csh, 0)
+
+# plot
+fig, axs = plt.subplots(2,2, figsize=(12,10))
+# reflectivity
+rpm = axs[0,0].pcolormesh(x/1000, y/1000, ref, vmin=0, vmax=45, cmap='magma_r')
+axs[0,0].set_aspect('equal')
+axs[0,0].set_title('Reflectivity [dBZ]')
+plt.colorbar(rpm, ax=axs[0,0])
+# convsf
+csbpm = axs[0,1].pcolormesh(x/1000, y/1000, csb_masked, vmin=1, vmax=3, cmap=convsf_cmap)
+axs[0,1].set_aspect('equal')
+axs[0,1].set_title('convsf')
+cb = plt.colorbar(csbpm, ax=axs[0,1], ticks=[4/3, 2, 8/3])
+cb.ax.set_yticklabels(['Stratiform', 'Convective', 'Weak Echo'])
+# convsf lo
+cslpm = axs[1,0].pcolormesh(x/1000, y/1000, csl_masked, vmin=1, vmax=3, cmap=convsf_cmap)
+axs[1,0].set_aspect('equal')
+axs[1,0].set_title('convsf lo')
+cb = plt.colorbar(cslpm, ax=axs[1,0], ticks=[])
+# convsf hi
+cshpm = axs[1,1].pcolormesh(x/1000, y/1000, csh_masked, vmin=1, vmax=3, cmap=convsf_cmap)
+axs[1,1].set_aspect('equal')
+axs[1,1].set_title('convsf hi')
+cb = plt.colorbar(cshpm, ax=axs[1,1], ticks=[4/3, 2, 8/3])
+cb.ax.set_yticklabels(['Stratiform', 'Convective', 'Weak Echo'])
+plt.show()
+
+# now let's compare with the Python algorithm
+convsf_dict = pyart.retrieve.feature_detection(grid, dx=x[1]-x[0], dy=y[1]-y[0], level_m=None, always_core_thres=40,
+                                              bkg_rad_km=11, use_cosine=True, max_diff=8, zero_diff_cos_val=55,
+                                              scalar_diff=None, use_addition=False, calc_thres=0,
+                                              weak_echo_thres=15, min_val_used=0, dB_averaging=True,
+                                              remove_small_objects=False, min_km2_size=None, val_for_max_rad=30,
+                                              max_rad_km=5.0, core_val=3, nosfcecho=0, weakecho=3, bkgd_val=1,
+                                              feat_val=2, field='maxdz', estimate_flag=True, estimate_offset=5)
+
+# new data
+csb_lt = convsf_dict['feature_detection']['data'][0, :, :]
+csb_lt_masked = np.ma.masked_less_equal(csb_lt, 0)
+csu_lt = convsf_dict['feature_under']['data'][0, :, :]
+csu_lt_masked = np.ma.masked_less_equal(csu_lt, 0)
+cso_lt = convsf_dict['feature_over']['data'][0, :, :]
+cso_lt_masked = np.ma.masked_less_equal(cso_lt, 0)
+
+fig, axs = plt.subplots(2,2, figsize=(12,10))
+# reflectivity
+rpm = axs[0,0].pcolormesh(x/1000, y/1000, ref, vmin=0, vmax=45, cmap='magma_r')
+axs[0,0].set_aspect('equal')
+axs[0,0].set_title('Reflectivity [dBZ]')
+plt.colorbar(rpm, ax=axs[0,0])
+# convsf best
+csbpm = axs[0,1].pcolormesh(x/1000, y/1000, csb_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap)
+axs[0,1].set_aspect('equal')
+axs[0,1].set_title('convsf')
+cb = plt.colorbar(csbpm, ax=axs[0,1], ticks=[4/3, 2, 8/3])
+cb.ax.set_yticklabels(['Stratiform', 'Convective', 'Weak Echo'])
+# convsf under
+csupm = axs[1,0].pcolormesh(x/1000, y/1000, csu_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap)
+axs[1,0].set_aspect('equal')
+axs[1,0].set_title('convsf under')
+cb = plt.colorbar(csupm, ax=axs[1,0], ticks=[])
+# convsf over
+csopm = axs[1,1].pcolormesh(x/1000, y/1000, cso_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap)
+axs[1,1].set_aspect('equal')
+axs[1,1].set_title('convsf over')
+cb = plt.colorbar(csopm, ax=axs[1,1], ticks=[4/3, 2, 8/3])
+cb.ax.set_yticklabels(['Stratiform', 'Convective', 'Weak Echo'])
+plt.show()
+
+######################################
+# Part 2: Cool-season feature detection
+# ----------
+######################################
+# **Winter storm example**
+#
+# In this example, we will show how to algorithm can be used to detect features (snow bands) in winter storms. Here, we
+# rescale the reflectivity to snow rate (Rasumussen et al. 2003) in order to better represent the snow field. Note
+# that we are not using the absolute values of the snow rate, we are only using the relative anomalies relative to
+# the background average. We recommend using a rescaled reflectivity for this algorithm, but note that you need to
+# change ``dB_averaging = False`` when using a rescaled field (set ``dB_averaging`` to True for reflectivity fields in
+# dBZ units). Also note that since we are now using a snow rate field with different values, most of our parameters
+# have changed as well.
+
+# Read in file
+nexrad_file = "s3://noaa-nexrad-level2/2021/02/07/KOKX/KOKX20210207_161413_V06"
+radar = pyart.io.read_nexrad_archive(nexrad_file)
+
+# extract the lowest sweep
+radar = radar.extract_sweeps([0])
+
+# interpolate to grid
+grid = pyart.map.grid_from_radars(
+    (radar,),
+    grid_shape=(1, 201, 201),
+    grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)),
+    fields=["reflectivity", "cross_correlation_ratio"],
+)
+
+# image mute grid object
+grid = pyart.util.image_mute_radar(
+    grid, "reflectivity", "cross_correlation_ratio", 0.97, 20
+)
+
+# convect non-muted reflectivity to snow rate
+ref = grid.fields["reflectivity"]["data"][0, :, :]
+ref = np.ma.masked_invalid(ref)
+
+ref_linear = 10 ** (ref / 10)  # mm6/m3
+snow_rate = (ref_linear / 57.3) ** (1 / 1.67)  # mm/hr
+
+# add to grid
+snow_rate_dict = {
+    "data": snow_rate[None, :, :],
+    "standard_name": "snow_rate",
+    "long_name": "Snow rate converted from linear reflectivity",
+    "units": "mm/hr",
+    "valid_min": 0,
+    "valid_max": 40500,
+}
+grid.add_field("snow_rate", snow_rate_dict, replace_existing=True)
+
+# get dx dy
+dx = grid.x["data"][1] - grid.x["data"][0]
+dy = grid.y["data"][1] - grid.y["data"][0]
+
+# feature detection
+feature_dict = pyart.retrieve.feature_detection(
+    grid,
+    dx,
+    dy,
+    field="snow_rate",
+    dB_averaging=False,
+    always_core_thres=4,
+    bkg_rad_km=40,
+    use_cosine=True,
+    max_diff=1.5,
+    zero_diff_cos_val=5,
+    weak_echo_thres=0,
+    min_val_used=0,
+    max_rad_km=1,
+    estimate_flag=False,
+)
+
+# add to grid object
+# mask zero values (no surface echo)
+feature_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0)
+# mask 3 values (weak echo)
+feature_masked = np.ma.masked_equal(feature_masked, 3)
+# add dimension to array to add to grid object
+feature_dict["feature_detection"]["data"] = feature_masked
+# add field
+grid.add_field("feature_detection", feature_dict["feature_detection"], replace_existing=True)
+
+# create plot using GridMapDisplay
+# plot variables
+display = pyart.graph.GridMapDisplay(grid)
+
+projection = ccrs.AlbersEqualArea(
+    central_latitude=radar.latitude["data"][0],
+    central_longitude=radar.longitude["data"][0],
+)
+# plot
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(1, 2, 1, projection=projection)
+display.plot_grid(
+    "snow_rate",
+    vmin=0,
+    vmax=10,
+    cmap=plt.get_cmap("viridis"),
+    transform=ccrs.PlateCarree(),
+    ax=ax1,
+    axislabels_flag=False,
+)
+ax2 = plt.subplot(1, 2, 2, projection=projection)
+display.plot_grid(
+    "feature_detection",
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
+    axislabels_flag=False,
+    transform=ccrs.PlateCarree(),
+    ticks=[1 / 3, 1, 5 / 3],
+    ticklabs=["", "Background", "Feature"],
+    ax=ax2,
+)
+plt.show()
+
+######################################
+# **Strong and Faint features**
+#
+# We have developed a new technique which runs the algorithm twice to find features that are very distinct from the
+# background (strong features) and features that are not very distinct from the background (faint features).
+
+# run the algorithm again, but to find objects that are also faint
+# feature detection
+feature_dict = pyart.retrieve.feature_detection(
+    grid,
+    dx,
+    dy,
+    field="snow_rate",
+    dB_averaging=False,
+    always_core_thres=4,
+    bkg_rad_km=40,
+    use_cosine=False,
+    scalar_diff=1.5,
+    use_addition=False,
+    weak_echo_thres=0,
+    min_val_used=0,
+    max_rad_km=1,
+    estimate_flag=False,
+)
+
+# separate strong vs. faint
+# mask zero values (no surface echo)
+scalar_features_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0)
+# mask 3 values (weak echo)
+scalar_features_masked = np.ma.masked_equal(scalar_features_masked, 3)
+# get cosine features
+cosine_features_masked = grid.fields['feature_detection']['data']
+# isolate faint features only
+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)
+
+# add to grid object
+new_dict = {
+    "features_faint_strong" : {
+        "data": features_faint_strong_masked,
+        "standard_name": "features_faint_strong",
+        "long name": "Faint and Strong Features",
+        "valid_min": 0,
+        "valid_max": 10,
+        "comment_1": "3: Faint features, 2: Strong features, 1: background"
+    }
+}
+# add field
+grid.add_field("features_faint_strong", new_dict["features_faint_strong"], replace_existing=True)
+
+# create plot using GridMapDisplay
+# plot variables
+display = pyart.graph.GridMapDisplay(grid)
+
+projection = ccrs.AlbersEqualArea(
+    central_latitude=radar.latitude["data"][0],
+    central_longitude=radar.longitude["data"][0],
+)
+
+faint_strong_cmap = mcolors.LinearSegmentedColormap.from_list('faint_strong',
+                                                [(32 / 255, 144 / 255, 140 / 255), (254 / 255, 238 / 255, 93 / 255),
+                                                 (254 / 255, 152 / 255, 93 / 255)], N=3)
+
+# plot
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(1, 2, 1, projection=projection)
+display.plot_grid(
+    "snow_rate",
+    vmin=0,
+    vmax=10,
+    cmap=plt.get_cmap("viridis"),
+    transform=ccrs.PlateCarree(),
+    ax=ax1,
+    axislabels_flag=False,
+)
+ax2 = plt.subplot(1, 2, 2, projection=projection)
+display.plot_grid(
+    "features_faint_strong",
+    vmin=1,
+    vmax=3,
+    cmap=faint_strong_cmap,
+    axislabels_flag=False,
+    transform=ccrs.PlateCarree(),
+    ticks=[1.33, 2, 2.66],
+    ticklabs=["Background", "Strong Feature", "Faint Feature"],
+    ax=ax2,
+)
+plt.show()
+######################################
+# Summary of recommendations and best practices
+# ----------
+# * Tune your parameters to your specific purpose
+# * Use a rescaled field if possible (i.e. linear reflectivity, rain or snow rate)
+# * Keep ``estimate_flag=True`` to see uncertainty in classification
+#
+# References
+# ----------
+# Steiner, M. R., R. A. Houze Jr., and S. E. Yuter, 1995: Climatological
+# Characterization of Three-Dimensional Storm Structure from Operational
+# Radar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007.
+# https://doi.org/10.1175/1520-0450(1995)034<1978:CCOTDS>2.0.CO;2.
+#
+# Yuter, S. E., and R. A. Houze, Jr., 1997: Measurements of raindrop size
+# distributions over the Pacific warm pool and implications for Z-R relations.
+# J. Appl. Meteor., 36, 847-867.
+# https://doi.org/10.1175/1520-0450(1997)036%3C0847:MORSDO%3E2.0.CO;2
+#
+# Yuter, S. E., R. A. Houze, Jr., E. A. Smith, T. T. Wilheit, and E. Zipser,
+# 2005: Physical characterization of tropical oceanic convection observed in
+# KWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1
+#
+# Rasmussen, R., M. Dixon, S. Vasiloff, F. Hage, S. Knight, J. Vivekanandan,
+# and M. Xu, 2003: Snow Nowcasting Using a Real-Time Correlation of Radar
+# Reflectivity with Snow Gauge Accumulation. J. Appl. Meteorol. Climatol., 42, 20–36.
+# https://doi.org/10.1175/1520-0450(2003)042%3C0020:SNUART%3E2.0.CO;2

From 2bff81ac0d8c33b3a51c4d6485d43994d52956ad Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Thu, 9 Nov 2023 11:37:57 -0500
Subject: [PATCH 09/21] Delete plot_convective_stratiform.py

---
 .../retrieve/plot_convective_stratiform.py    | 467 ------------------
 1 file changed, 467 deletions(-)
 delete mode 100644 examples/retrieve/plot_convective_stratiform.py

diff --git a/examples/retrieve/plot_convective_stratiform.py b/examples/retrieve/plot_convective_stratiform.py
deleted file mode 100644
index da0f22cc9f..0000000000
--- a/examples/retrieve/plot_convective_stratiform.py
+++ /dev/null
@@ -1,467 +0,0 @@
-"""
-=======================================
-Convective-Stratiform classification
-=======================================
-This example shows how to use the updated convective stratiform classifcation algorithm. We show 3 examples,
-a summer convective example, an example from Hurricane Ian, and an example from a winter storm.
-"""
-
-print(__doc__)
-
-# Author: Laura Tomkins (lmtomkin@ncsu.edu)
-# License: BSD 3 clause
-
-
-import cartopy.crs as ccrs
-import matplotlib.colors as mcolors
-import matplotlib.pyplot as plt
-import numpy as np
-
-import pyart
-
-######################################
-# How the algorithm works
-# ----------
-# This first section describes how the convective-stratiform algorithm works (see references for full details). This
-# algorithm is a feature detection algorithm and classifies fields as "convective" or "stratiform". The algorithm is
-# designed to detect features in a reflectivity field but can also detect features in fields such as rain rate or
-# snow rate. In this section we describe the steps of the convective stratiform algorithm and the variables used in
-# the function.
-# The first step of the algorithm calculates a background average of the field with a circular footprint using a radius
-# provided by ``bkg_rad_km``. A larger radius will yield a smoother field. The radius needs to be at least double the
-# grid spacing, but we recommend at least three times the grid spacing. If using reflectivity, ``dB_averaging`` should be set
-# to True to convert reflectivity to linear Z before averaging, False for rescaled fields such as rain or snow rate.
-# ``calc_thres`` determines the minimum fraction of a circle that is considered in the background average calculation
-# (default is 0.75, so the points along the edges where there is less than 75% of a full circle of data,
-# the algorithm is not run).
-# Once the background average has been calculated, the original field is compared to the background average.  In
-# order for points to be considered "convective cores" they must exceed the background value by a certain value or
-# simply be greater than the ``always_core_thres``. This value is determined by either a cosine scheme, or a scalar
-# value (i.e. the reflectivity value must be X times the background value (multiplier; ``use_addition=False``),
-# or  X greater than the background value (``use_addition=True``) where X is the ``scalar_diff``).
-# ``use_cosine`` determines if a cosine scheme or scalar scheme is to be used. If ``use_cosine`` is True,
-# then the ``max_diff`` and ``zero_diff_cos_val`` come into use. These values define the cosine scheme that is used  to
-# determine the minimum difference between the background and reflectivity value in order for a core to be
-# identified. ``max_diff`` is the maximum difference between the field and the background for a core to be identified,
-# or where the cosine function crosses the y-axis. The ``zero_diff_cos_val`` is where the difference between the  field
-# and the background is zero, or where the cosine function crosses the x-axis. Note, if
-# ``always_core_thres`` < ``zero_diff_cos_val``, ``zero_diff_cos_val`` only helps define the shape of the cosine curve and
-# all values greater than ``always_core_thres`` will be considered a convective core. If
-# ``always_core_thres`` > ``zero_diff_cos_val`` then all values greater than ``zero_diff_cos_val`` will be considered a
-# convective core. We plot some examples of the schemes below:
-
-######################################
-# Example of the cosine scheme:
-pyart.graph.plot_convstrat_scheme(
-    always_core_thres=30, use_cosine=True, max_diff=5, zero_diff_cos_val=45
-)
-
-######################################
-# when zero_diff_cos_val is greater than always_core_thres, the difference becomes zero at the zero_diff_cos_val
-pyart.graph.plot_convstrat_scheme(
-    always_core_thres=55, use_cosine=True, max_diff=5, zero_diff_cos_val=45
-)
-
-######################################
-# alternatively, we can use a simpler scalar difference instead of a cosine scheme
-pyart.graph.plot_convstrat_scheme(
-    always_core_thres=40,
-    use_cosine=False,
-    max_diff=None,
-    zero_diff_cos_val=None,
-    use_addition=True,
-    scalar_diff=2,
-)
-
-######################################
-# if you are interested in picking up weak features, you can also use the scalar difference as a multiplier instead,
-# so very weak features do not have to be that different from the background to be classified as convective.
-pyart.graph.plot_convstrat_scheme(
-    always_core_thres=40,
-    use_cosine=False,
-    max_diff=None,
-    zero_diff_cos_val=None,
-    use_addition=False,
-    scalar_diff=2,
-)
-
-######################################
-# Once the cores are identified, there is an option to remove speckles (``remove_small_objects``) smaller than a  given
-# size (``min_km2_size``).
-# After the convective cores are identified, We then incorporate convective radii using
-# ``val_for_max_conv_rad`` and ``max_conv_rad_km``. The convective radii act as a dilation and are used to classify
-# additional points around the cores as convective that may not have been identified previously.  The
-# ``val_for_max_conv_rad`` is the value where the maximum convective radius is applied and the ``max_conv_rad_km`` is the
-# maximum convective radius. Values less than the ``val_for_max_conv_rad`` are assigned a convective radius using a step
-# function.
-# Finally, the points are classified as NOSFCECHO (threshold set with ``min_dBZ_used``; 0), WEAKECHO (threshold set with
-# ``weak_echo_thres``; 3), SF (stratiform; 1), CONV (convective; 2).
-
-######################################
-# Examples
-# ----------
-# **Classification of summer convective example**
-#
-# Our first example classifies echo from a summer convective event. We use a cosine scheme to classify the convective
-# points.
-
-# Now let's do a classification with our parameters
-# read in file
-filename = pyart.testing.get_test_data("swx_20120520_0641.nc")
-radar = pyart.io.read(filename)
-
-# extract the lowest sweep
-radar = radar.extract_sweeps([0])
-
-# interpolate to grid
-grid = pyart.map.grid_from_radars(
-    (radar,),
-    grid_shape=(1, 201, 201),
-    grid_limits=((0, 10000), (-50000.0, 50000.0), (-50000.0, 50000.0)),
-    fields=["reflectivity_horizontal"],
-)
-
-# get dx dy
-dx = grid.x["data"][1] - grid.x["data"][0]
-dy = grid.y["data"][1] - grid.y["data"][0]
-
-# convective stratiform classification
-convsf_dict = pyart.retrieve.conv_strat_yuter(
-    grid,
-    dx,
-    dy,
-    refl_field="reflectivity_horizontal",
-    always_core_thres=40,
-    bkg_rad_km=20,
-    use_cosine=True,
-    max_diff=5,
-    zero_diff_cos_val=55,
-    weak_echo_thres=10,
-    max_conv_rad_km=2,
-)
-
-# add to grid object
-# mask zero values (no surface echo)
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
-# mask 3 values (weak echo)
-convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-# add dimension to array to add to grid object
-convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
-# add field
-grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
-
-# create plot using GridMapDisplay
-# plot variables
-display = pyart.graph.GridMapDisplay(grid)
-magma_r_cmap = plt.get_cmap("magma_r")
-ref_cmap = mcolors.LinearSegmentedColormap.from_list(
-    "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))
-)
-projection = ccrs.AlbersEqualArea(
-    central_latitude=radar.latitude["data"][0],
-    central_longitude=radar.longitude["data"][0],
-)
-
-# plot
-plt.figure(figsize=(10, 4))
-ax1 = plt.subplot(1, 2, 1, projection=projection)
-display.plot_grid(
-    "reflectivity_horizontal",
-    vmin=5,
-    vmax=45,
-    cmap=ref_cmap,
-    transform=ccrs.PlateCarree(),
-    ax=ax1,
-)
-ax2 = plt.subplot(1, 2, 2, projection=projection)
-display.plot_grid(
-    "convsf",
-    vmin=0,
-    vmax=2,
-    cmap=plt.get_cmap("viridis", 3),
-    ax=ax2,
-    transform=ccrs.PlateCarree(),
-    ticks=[1 / 3, 1, 5 / 3],
-    ticklabs=["", "Stratiform", "Convective"],
-)
-plt.show()
-
-
-######################################
-# In addition to the default convective-stratiform classification, the function also returns an underestimate
-# (convsf_under) and an overestimate (convsf_over) to take into consideration the uncertainty when choosing
-# classification parameters. The under and overestimate use the same parameters, but vary the input field by a
-# certain value (default is 5 dBZ, can be changed with ``estimate_offset``). The estimation can be turned off (
-# ``estimate_flag=False``), but we recommend keeping it turned on.
-
-# mask weak echo and no surface echo
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
-convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-convsf_dict["convsf"]["data"] = convsf_masked
-# underest.
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf_under"]["data"], 0)
-convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-convsf_dict["convsf_under"]["data"] = convsf_masked
-# overest.
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf_over"]["data"], 0)
-convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-convsf_dict["convsf_over"]["data"] = convsf_masked
-
-# Plot each estimation
-plt.figure(figsize=(10, 4))
-ax1 = plt.subplot(131)
-ax1.pcolormesh(
-    convsf_dict["convsf"]["data"][0, :, :],
-    vmin=0,
-    vmax=2,
-    cmap=plt.get_cmap("viridis", 3),
-)
-ax1.set_title("Best estimate")
-ax1.set_aspect("equal")
-ax2 = plt.subplot(132)
-ax2.pcolormesh(
-    convsf_dict["convsf_under"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
-)
-ax2.set_title("Underestimate")
-ax2.set_aspect("equal")
-ax3 = plt.subplot(133)
-ax3.pcolormesh(
-    convsf_dict["convsf_over"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
-)
-ax3.set_title("Overestimate")
-ax3.set_aspect("equal")
-plt.show()
-
-######################################
-# **Tropical example**
-#
-# Let's get a NEXRAD file from Hurricane Ian
-
-# Read in file
-nexrad_file = "s3://noaa-nexrad-level2/2022/09/28/KTBW/KTBW20220928_190142_V06"
-radar = pyart.io.read_nexrad_archive(nexrad_file)
-
-# extract the lowest sweep
-radar = radar.extract_sweeps([0])
-
-# interpolate to grid
-grid = pyart.map.grid_from_radars(
-    (radar,),
-    grid_shape=(1, 201, 201),
-    grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)),
-    fields=["reflectivity"],
-)
-
-# get dx dy
-dx = grid.x["data"][1] - grid.x["data"][0]
-dy = grid.y["data"][1] - grid.y["data"][0]
-
-# convective stratiform classification
-convsf_dict = pyart.retrieve.conv_strat_yuter(
-    grid,
-    dx,
-    dy,
-    refl_field="reflectivity",
-    always_core_thres=40,
-    bkg_rad_km=20,
-    use_cosine=True,
-    max_diff=3,
-    zero_diff_cos_val=55,
-    weak_echo_thres=5,
-    max_conv_rad_km=2,
-    estimate_flag=False,
-)
-
-# add to grid object
-# mask zero values (no surface echo)
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
-# mask 3 values (weak echo)
-convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-# add dimension to array to add to grid object
-convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
-# add field
-grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
-
-# create plot using GridMapDisplay
-# plot variables
-display = pyart.graph.GridMapDisplay(grid)
-magma_r_cmap = plt.get_cmap("magma_r")
-ref_cmap = mcolors.LinearSegmentedColormap.from_list(
-    "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))
-)
-projection = ccrs.AlbersEqualArea(
-    central_latitude=radar.latitude["data"][0],
-    central_longitude=radar.longitude["data"][0],
-)
-# plot
-plt.figure(figsize=(10, 4))
-ax1 = plt.subplot(1, 2, 1, projection=projection)
-display.plot_grid(
-    "reflectivity",
-    vmin=5,
-    vmax=45,
-    cmap=ref_cmap,
-    transform=ccrs.PlateCarree(),
-    ax=ax1,
-    axislabels_flag=False,
-)
-ax2 = plt.subplot(1, 2, 2, projection=projection)
-display.plot_grid(
-    "convsf",
-    vmin=0,
-    vmax=2,
-    cmap=plt.get_cmap("viridis", 3),
-    axislabels_flag=False,
-    transform=ccrs.PlateCarree(),
-    ticks=[1 / 3, 1, 5 / 3],
-    ticklabs=["", "Stratiform", "Convective"],
-    ax=ax2,
-)
-plt.show()
-
-######################################
-# **Winter storm example with image muting**
-#
-# Here is a final example of the convective stratiform classification using an example from a winter storm. Before
-# doing the classification, we image mute the reflectivity to remove regions with melting or mixed precipitation. We
-# then rescale the reflectivity to snow rate (Rasumussen et al. 2003). We recommend using a rescaled reflectivity
-# to do the classification, but if you do make sure to changed dB_averaging to False because this parameter is  used
-# to convert reflectivity to a linear value before averaging (set dB_averaging to True for reflectivity fields in
-# dBZ units).
-# In this example, note how we change some of the other parameters since we are classifying snow rate instead of
-# reflecitivity.
-
-# Read in file
-nexrad_file = "s3://noaa-nexrad-level2/2021/02/07/KOKX/KOKX20210207_161413_V06"
-radar = pyart.io.read_nexrad_archive(nexrad_file)
-
-# extract the lowest sweep
-radar = radar.extract_sweeps([0])
-
-# interpolate to grid
-grid = pyart.map.grid_from_radars(
-    (radar,),
-    grid_shape=(1, 201, 201),
-    grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)),
-    fields=["reflectivity", "cross_correlation_ratio"],
-)
-
-# image mute grid object
-grid = pyart.util.image_mute_radar(
-    grid, "reflectivity", "cross_correlation_ratio", 0.97, 20
-)
-
-# convect non-muted reflectivity to snow rate
-nonmuted_ref = grid.fields["nonmuted_reflectivity"]["data"][0, :, :]
-nonmuted_ref = np.ma.masked_invalid(nonmuted_ref)
-
-nonmuted_ref_linear = 10 ** (nonmuted_ref / 10)  # mm6/m3
-snow_rate = (nonmuted_ref_linear / 57.3) ** (1 / 1.67)  #
-
-# add to grid
-snow_rate_dict = {
-    "data": snow_rate[None, :, :],
-    "standard_name": "snow_rate",
-    "long_name": "Snow rate converted from linear reflectivity",
-    "units": "mm/hr",
-    "valid_min": 0,
-    "valid_max": 40500,
-}
-grid.add_field("snow_rate", snow_rate_dict, replace_existing=True)
-
-# get dx dy
-dx = grid.x["data"][1] - grid.x["data"][0]
-dy = grid.y["data"][1] - grid.y["data"][0]
-
-# convective stratiform classification
-convsf_dict = pyart.retrieve.conv_strat_yuter(
-    grid,
-    dx,
-    dy,
-    refl_field="snow_rate",
-    dB_averaging=False,
-    always_core_thres=4,
-    bkg_rad_km=40,
-    use_cosine=True,
-    max_diff=1.5,
-    zero_diff_cos_val=5,
-    weak_echo_thres=0,
-    min_dBZ_used=0,
-    max_conv_rad_km=1,
-    estimate_flag=False,
-)
-
-# add to grid object
-# mask zero values (no surface echo)
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
-# mask 3 values (weak echo)
-convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-# add dimension to array to add to grid object
-convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
-# add field
-grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
-
-# create plot using GridMapDisplay
-# plot variables
-display = pyart.graph.GridMapDisplay(grid)
-magma_r_cmap = plt.get_cmap("magma_r")
-ref_cmap = mcolors.LinearSegmentedColormap.from_list(
-    "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))
-)
-projection = ccrs.AlbersEqualArea(
-    central_latitude=radar.latitude["data"][0],
-    central_longitude=radar.longitude["data"][0],
-)
-# plot
-plt.figure(figsize=(10, 4))
-ax1 = plt.subplot(1, 2, 1, projection=projection)
-display.plot_grid(
-    "snow_rate",
-    vmin=0,
-    vmax=10,
-    cmap=plt.get_cmap("viridis"),
-    transform=ccrs.PlateCarree(),
-    ax=ax1,
-    axislabels_flag=False,
-)
-ax2 = plt.subplot(1, 2, 2, projection=projection)
-display.plot_grid(
-    "convsf",
-    vmin=0,
-    vmax=2,
-    cmap=plt.get_cmap("viridis", 3),
-    axislabels_flag=False,
-    transform=ccrs.PlateCarree(),
-    ticks=[1 / 3, 1, 5 / 3],
-    ticklabs=["", "Stratiform", "Convective"],
-    ax=ax2,
-)
-plt.show()
-
-######################################
-# Summary of recommendations and best practices
-# ----------
-# * Tune your parameters to your specific purpose
-# * Use a rescaled field if possible (i.e. linear reflectivity, rain or snow rate)
-# * Keep ``estimate_flag=True`` to see uncertainty in classification
-#
-# References
-# ----------
-# Steiner, M. R., R. A. Houze Jr., and S. E. Yuter, 1995: Climatological
-# Characterization of Three-Dimensional Storm Structure from Operational
-# Radar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007.
-# https://doi.org/10.1175/1520-0450(1995)034<1978:CCOTDS>2.0.CO;2.
-#
-# Yuter, S. E., and R. A. Houze, Jr., 1997: Measurements of raindrop size
-# distributions over the Pacific warm pool and implications for Z-R relations.
-# J. Appl. Meteor., 36, 847-867.
-# https://doi.org/10.1175/1520-0450(1997)036%3C0847:MORSDO%3E2.0.CO;2
-#
-# Yuter, S. E., R. A. Houze, Jr., E. A. Smith, T. T. Wilheit, and E. Zipser,
-# 2005: Physical characterization of tropical oceanic convection observed in
-# KWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1
-#
-# Rasmussen, R., M. Dixon, S. Vasiloff, F. Hage, S. Knight, J. Vivekanandan,
-# and M. Xu, 2003: Snow Nowcasting Using a Real-Time Correlation of Radar
-# Reflectivity with Snow Gauge Accumulation. J. Appl. Meteorol. Climatol., 42, 20–36.
-# https://doi.org/10.1175/1520-0450(2003)042%3C0020:SNUART%3E2.0.CO;2

From 11c6c40754c0ac37d9eba9bc3550de224519f708 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Mon, 13 Nov 2023 14:39:01 -0500
Subject: [PATCH 10/21] Update echo_class.py

---
 pyart/retrieve/echo_class.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index dd524eaadd..a52f1a2063 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -255,8 +255,8 @@ def feature_detection(
     2005: Physical characterization of tropical oceanic convection observed in
     KWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1
 
-    Tomkins, L. M., S. E. Yuter, and M. A. Miller, 2023: Objective identification
-    of faint and strong snow bands in radar observations of winter storms. in prep.
+    Tomkins, L. M., S. E. Yuter, and M. A. Miller, 2024: Objective identification
+    of faint and strong features in radar observations of winter storms. in prep.
 
     """
 
@@ -342,7 +342,6 @@ def feature_detection(
 
     # If estimation is True, run the algorithm on the field with offset subtracted and the field with the offset added
     if estimate_flag:
-
         # get underestimate field
         if underest_field is None:
             under_field = ze - estimate_offset

From 0562ad3b40a374d9c2823c77e9fab6d0fc2d099d Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Mon, 13 Nov 2023 15:10:11 -0500
Subject: [PATCH 11/21] Update plot_feature_detection.py

---
 examples/retrieve/plot_feature_detection.py | 381 +++++++++++++++-----
 1 file changed, 291 insertions(+), 90 deletions(-)

diff --git a/examples/retrieve/plot_feature_detection.py b/examples/retrieve/plot_feature_detection.py
index c7fb9ec29b..f4489302b8 100644
--- a/examples/retrieve/plot_feature_detection.py
+++ b/examples/retrieve/plot_feature_detection.py
@@ -17,9 +17,9 @@
 import matplotlib.colors as mcolors
 import matplotlib.pyplot as plt
 import numpy as np
+from open_radar_data import DATASETS
 
 import pyart
-from open_radar_data import DATASETS
 
 ######################################
 # How the algorithm works
@@ -80,7 +80,9 @@
 # add dimension to array to add to grid object
 feature_dict["feature_detection"]["data"] = feature_masked
 # add field
-grid.add_field("feature_detection", feature_dict["feature_detection"], replace_existing=True)
+grid.add_field(
+    "feature_detection", feature_dict["feature_detection"], replace_existing=True
+)
 
 # create plot using GridMapDisplay
 # plot variables
@@ -118,7 +120,6 @@
 )
 plt.show()
 
-
 ######################################
 # In addition to the default feature detection classification, the function also returns an underestimate
 # (``feature_under``) and an overestimate (``feature_over``) to take into consideration the uncertainty when choosing
@@ -155,7 +156,7 @@
     feature_dict["feature_under"]["data"][0, :, :],
     vmin=0,
     vmax=2,
-    cmap=plt.get_cmap("viridis", 3)
+    cmap=plt.get_cmap("viridis", 3),
 )
 ax2.set_title("Underestimate")
 ax2.set_aspect("equal")
@@ -164,7 +165,7 @@
     feature_dict["feature_over"]["data"][0, :, :],
     vmin=0,
     vmax=2,
-    cmap=plt.get_cmap("viridis", 3)
+    cmap=plt.get_cmap("viridis", 3),
 )
 ax3.set_title("Overestimate")
 ax3.set_aspect("equal")
@@ -218,7 +219,9 @@
 # add dimension to array to add to grid object
 feature_dict["feature_detection"]["data"] = feature_masked
 # add field
-grid.add_field("feature_detection", feature_dict["feature_detection"], replace_existing=True)
+grid.add_field(
+    "feature_detection", feature_dict["feature_detection"], replace_existing=True
+)
 
 # create plot using GridMapDisplay
 # plot variables
@@ -257,7 +260,6 @@
 )
 plt.show()
 
-
 ######################################
 # **Comparison with original C++ algorithm**
 #
@@ -270,87 +272,127 @@
 grid = pyart.io.read_grid(filename)
 
 # colormap
-convsf_cmap = mcolors.LinearSegmentedColormap.from_list('convsf', [(33 / 255, 145 / 255, 140 / 255),
-                                                                    (253 / 255, 231 / 255, 37 / 255),
-                                                                    (210 / 255, 180 / 255, 140 / 255)], N=3)
+convsf_cmap = mcolors.LinearSegmentedColormap.from_list(
+    "convsf",
+    [
+        (33 / 255, 145 / 255, 140 / 255),
+        (253 / 255, 231 / 255, 37 / 255),
+        (210 / 255, 180 / 255, 140 / 255)
+    ],
+    N=3,
+)
 
 # get original data from file processed in C++
-x = grid.x['data']
-y = grid.y['data']
-ref = grid.fields['maxdz']['data'][0,:,:]
-csb = grid.fields['convsf']['data'][0,:,:]
+x = grid.x["data"]
+y = grid.y["data"]
+ref = grid.fields["maxdz"]["data"][0, :, :]
+csb = grid.fields["convsf"]["data"][0, :, :]
 csb_masked = np.ma.masked_less_equal(csb, 0)
-csl = grid.fields['convsf_lo']['data'][0,:,:]
+csl = grid.fields["convsf_lo"]["data"][0, :, :]
 csl_masked = np.ma.masked_less_equal(csl, 0)
-csh = grid.fields['convsf_hi']['data'][0,:,:]
+csh = grid.fields["convsf_hi"]["data"][0, :, :]
 csh_masked = np.ma.masked_less_equal(csh, 0)
 
 # plot
-fig, axs = plt.subplots(2,2, figsize=(12,10))
+fig, axs = plt.subplots(2, 2, figsize=(12, 10))
 # reflectivity
-rpm = axs[0,0].pcolormesh(x/1000, y/1000, ref, vmin=0, vmax=45, cmap='magma_r')
-axs[0,0].set_aspect('equal')
-axs[0,0].set_title('Reflectivity [dBZ]')
-plt.colorbar(rpm, ax=axs[0,0])
+rpm = axs[0, 0].pcolormesh(x / 1000, y / 1000, ref, vmin=0, vmax=45, cmap="magma_r")
+axs[0, 0].set_aspect("equal")
+axs[0, 0].set_title("Reflectivity [dBZ]")
+plt.colorbar(rpm, ax=axs[0, 0])
 # convsf
-csbpm = axs[0,1].pcolormesh(x/1000, y/1000, csb_masked, vmin=1, vmax=3, cmap=convsf_cmap)
-axs[0,1].set_aspect('equal')
-axs[0,1].set_title('convsf')
-cb = plt.colorbar(csbpm, ax=axs[0,1], ticks=[4/3, 2, 8/3])
-cb.ax.set_yticklabels(['Stratiform', 'Convective', 'Weak Echo'])
+csbpm = axs[0, 1].pcolormesh(
+    x / 1000, y / 1000, csb_masked, vmin=1, vmax=3, cmap=convsf_cmap
+)
+axs[0, 1].set_aspect("equal")
+axs[0, 1].set_title("convsf")
+cb = plt.colorbar(csbpm, ax=axs[0, 1], ticks=[4 / 3, 2, 8 / 3])
+cb.ax.set_yticklabels(["Stratiform", "Convective", "Weak Echo"])
 # convsf lo
-cslpm = axs[1,0].pcolormesh(x/1000, y/1000, csl_masked, vmin=1, vmax=3, cmap=convsf_cmap)
-axs[1,0].set_aspect('equal')
-axs[1,0].set_title('convsf lo')
-cb = plt.colorbar(cslpm, ax=axs[1,0], ticks=[])
+cslpm = axs[1, 0].pcolormesh(
+    x / 1000, y / 1000, csl_masked, vmin=1, vmax=3, cmap=convsf_cmap
+)
+axs[1, 0].set_aspect("equal")
+axs[1, 0].set_title("convsf lo")
+cb = plt.colorbar(cslpm, ax=axs[1, 0], ticks=[])
 # convsf hi
-cshpm = axs[1,1].pcolormesh(x/1000, y/1000, csh_masked, vmin=1, vmax=3, cmap=convsf_cmap)
-axs[1,1].set_aspect('equal')
-axs[1,1].set_title('convsf hi')
-cb = plt.colorbar(cshpm, ax=axs[1,1], ticks=[4/3, 2, 8/3])
-cb.ax.set_yticklabels(['Stratiform', 'Convective', 'Weak Echo'])
+cshpm = axs[1, 1].pcolormesh(
+    x / 1000, y / 1000, csh_masked, vmin=1, vmax=3, cmap=convsf_cmap
+)
+axs[1, 1].set_aspect("equal")
+axs[1, 1].set_title("convsf hi")
+cb = plt.colorbar(cshpm, ax=axs[1, 1], ticks=[4 / 3, 2, 8 / 3])
+cb.ax.set_yticklabels(["Stratiform", "Convective", "Weak Echo"])
 plt.show()
 
 # now let's compare with the Python algorithm
-convsf_dict = pyart.retrieve.feature_detection(grid, dx=x[1]-x[0], dy=y[1]-y[0], level_m=None, always_core_thres=40,
-                                              bkg_rad_km=11, use_cosine=True, max_diff=8, zero_diff_cos_val=55,
-                                              scalar_diff=None, use_addition=False, calc_thres=0,
-                                              weak_echo_thres=15, min_val_used=0, dB_averaging=True,
-                                              remove_small_objects=False, min_km2_size=None, val_for_max_rad=30,
-                                              max_rad_km=5.0, core_val=3, nosfcecho=0, weakecho=3, bkgd_val=1,
-                                              feat_val=2, field='maxdz', estimate_flag=True, estimate_offset=5)
+convsf_dict = pyart.retrieve.feature_detection(
+    grid,
+    dx=x[1] - x[0],
+    dy=y[1] - y[0],
+    level_m=None,
+    always_core_thres=40,
+    bkg_rad_km=11,
+    use_cosine=True,
+    max_diff=8,
+    zero_diff_cos_val=55,
+    scalar_diff=None,
+    use_addition=False,
+    calc_thres=0,
+    weak_echo_thres=15,
+    min_val_used=0,
+    dB_averaging=True,
+    remove_small_objects=False,
+    min_km2_size=None,
+    val_for_max_rad=30,
+    max_rad_km=5.0,
+    core_val=3,
+    nosfcecho=0,
+    weakecho=3,
+    bkgd_val=1,
+    feat_val=2,
+    field='maxdz',
+    estimate_flag=True,
+    estimate_offset=5,
+)
 
 # new data
-csb_lt = convsf_dict['feature_detection']['data'][0, :, :]
+csb_lt = convsf_dict["feature_detection"]["data"][0, :, :]
 csb_lt_masked = np.ma.masked_less_equal(csb_lt, 0)
-csu_lt = convsf_dict['feature_under']['data'][0, :, :]
+csu_lt = convsf_dict["feature_under"]["data"][0, :, :]
 csu_lt_masked = np.ma.masked_less_equal(csu_lt, 0)
-cso_lt = convsf_dict['feature_over']['data'][0, :, :]
+cso_lt = convsf_dict["feature_over"]["data"][0, :, :]
 cso_lt_masked = np.ma.masked_less_equal(cso_lt, 0)
 
-fig, axs = plt.subplots(2,2, figsize=(12,10))
+fig, axs = plt.subplots(2, 2, figsize=(12, 10))
 # reflectivity
-rpm = axs[0,0].pcolormesh(x/1000, y/1000, ref, vmin=0, vmax=45, cmap='magma_r')
-axs[0,0].set_aspect('equal')
-axs[0,0].set_title('Reflectivity [dBZ]')
-plt.colorbar(rpm, ax=axs[0,0])
+rpm = axs[0, 0].pcolormesh(x / 1000, y / 1000, ref, vmin=0, vmax=45, cmap="magma_r")
+axs[0, 0].set_aspect("equal")
+axs[0, 0].set_title("Reflectivity [dBZ]")
+plt.colorbar(rpm, ax=axs[0, 0])
 # convsf best
-csbpm = axs[0,1].pcolormesh(x/1000, y/1000, csb_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap)
-axs[0,1].set_aspect('equal')
-axs[0,1].set_title('convsf')
-cb = plt.colorbar(csbpm, ax=axs[0,1], ticks=[4/3, 2, 8/3])
-cb.ax.set_yticklabels(['Stratiform', 'Convective', 'Weak Echo'])
+csbpm = axs[0, 1].pcolormesh(
+    x / 1000, y / 1000, csb_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap
+)
+axs[0, 1].set_aspect("equal")
+axs[0, 1].set_title("convsf")
+cb = plt.colorbar(csbpm, ax=axs[0, 1], ticks=[4 / 3, 2, 8 / 3])
+cb.ax.set_yticklabels(["Stratiform", "Convective", "Weak Echo"])
 # convsf under
-csupm = axs[1,0].pcolormesh(x/1000, y/1000, csu_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap)
-axs[1,0].set_aspect('equal')
-axs[1,0].set_title('convsf under')
-cb = plt.colorbar(csupm, ax=axs[1,0], ticks=[])
+csupm = axs[1, 0].pcolormesh(
+    x / 1000, y / 1000, csu_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap
+)
+axs[1, 0].set_aspect("equal")
+axs[1, 0].set_title("convsf under")
+cb = plt.colorbar(csupm, ax=axs[1, 0], ticks=[])
 # convsf over
-csopm = axs[1,1].pcolormesh(x/1000, y/1000, cso_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap)
-axs[1,1].set_aspect('equal')
-axs[1,1].set_title('convsf over')
-cb = plt.colorbar(csopm, ax=axs[1,1], ticks=[4/3, 2, 8/3])
-cb.ax.set_yticklabels(['Stratiform', 'Convective', 'Weak Echo'])
+csopm = axs[1, 1].pcolormesh(
+    x / 1000, y / 1000, cso_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap
+)
+axs[1, 1].set_aspect("equal")
+axs[1, 1].set_title("convsf over")
+cb = plt.colorbar(csopm, ax=axs[1, 1], ticks=[4 / 3, 2, 8 / 3])
+cb.ax.set_yticklabels(["Stratiform", "Convective", "Weak Echo"])
 plt.show()
 
 ######################################
@@ -387,23 +429,8 @@
     grid, "reflectivity", "cross_correlation_ratio", 0.97, 20
 )
 
-# convect non-muted reflectivity to snow rate
-ref = grid.fields["reflectivity"]["data"][0, :, :]
-ref = np.ma.masked_invalid(ref)
-
-ref_linear = 10 ** (ref / 10)  # mm6/m3
-snow_rate = (ref_linear / 57.3) ** (1 / 1.67)  # mm/hr
-
-# add to grid
-snow_rate_dict = {
-    "data": snow_rate[None, :, :],
-    "standard_name": "snow_rate",
-    "long_name": "Snow rate converted from linear reflectivity",
-    "units": "mm/hr",
-    "valid_min": 0,
-    "valid_max": 40500,
-}
-grid.add_field("snow_rate", snow_rate_dict, replace_existing=True)
+# rescale reflectivity to 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]
@@ -435,7 +462,9 @@
 # add dimension to array to add to grid object
 feature_dict["feature_detection"]["data"] = feature_masked
 # add field
-grid.add_field("feature_detection", feature_dict["feature_detection"], replace_existing=True)
+grid.add_field(
+    "feature_detection", feature_dict["feature_detection"], replace_existing=True
+)
 
 # create plot using GridMapDisplay
 # plot variables
@@ -471,6 +500,104 @@
 )
 plt.show()
 
+######################################
+# **Under and over estimate in snow**
+#
+# Here we show how to bound the snow rate with an under and over estimate
+
+# get reflectivity field and copy
+ref_field = grid.fields["reflectivity"]
+ref_field_over = ref_field.copy()
+ref_field_under = ref_field.copy()
+
+# offset original field by +/- 2dB
+ref_field_over["data"] = ref_field["data"] + 2
+ref_field_under["data"] = ref_field["data"] - 2
+
+# adjust dictionary
+ref_field_over["standard_name"] = ref_field["standard_name"] + "_overest"
+ref_field_over["long_name"] = ref_field["long_name"] + " overestimate"
+ref_field_under["standard_name"] = ref_field["standard_name"] + "_underest"
+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
+)
+
+# convert to snow rate
+grid = pyart.retrieve.qpe.ZtoR(
+    grid, ref_field="reflectivity_over", a=57.3, b=1.67, save_name="snow_rate_over"
+)
+grid = pyart.retrieve.qpe.ZtoR(
+    grid, ref_field="reflectivity_under", a=57.3, b=1.67, save_name="snow_rate_under"
+)
+
+# now do feature detection with under and over estimate fields
+feature_estimate_dict = pyart.retrieve.feature_detection(
+    grid,
+    dx,
+    dy,
+    field="snow_rate",
+    dB_averaging=False,
+    always_core_thres=4,
+    bkg_rad_km=40,
+    use_cosine=True,
+    max_diff=1.5,
+    zero_diff_cos_val=5,
+    weak_echo_thres=0,
+    min_val_used=0,
+    max_rad_km=1,
+    estimate_flag=True,
+    overest_field="snow_rate_over",
+    underest_field="snow_rate_under",
+)
+
+# x and y data for plotting
+x = grid.x["data"]
+y = grid.y["data"]
+
+# now let's plot
+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,
+)
+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),
+)
+axs[0, 1].set_aspect("equal")
+axs[0, 1].set_title("Feature detection (best)")
+cb = plt.colorbar(bpm, ax=axs[0, 1], ticks=[3 / 2, 5 / 2])
+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),
+)
+axs[1, 0].set_aspect("equal")
+axs[1, 0].set_title("Feature detection (underestimate)")
+cb = plt.colorbar(upm, ax=axs[1, 0], ticks=[3 / 2, 5 / 2])
+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),
+)
+axs[1, 1].set_aspect("equal")
+axs[1, 1].set_title("Feature detection (overestimate)")
+cb = plt.colorbar(opm, ax=axs[1, 1], ticks=[3 / 2, 5 / 2])
+cb.ax.set_yticklabels(["Background", "Feature"])
+plt.show()
+
 ######################################
 # **Strong and Faint features**
 #
@@ -498,11 +625,13 @@
 
 # separate strong vs. faint
 # mask zero values (no surface echo)
-scalar_features_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0)
+scalar_features_masked = np.ma.masked_equal(
+    feature_dict["feature_detection"]["data"], 0
+)
 # mask 3 values (weak echo)
 scalar_features_masked = np.ma.masked_equal(scalar_features_masked, 3)
 # get cosine features
-cosine_features_masked = grid.fields['feature_detection']['data']
+cosine_features_masked = grid.fields["feature_detection"]["data"]
 # isolate faint features only
 faint_features = scalar_features_masked - cosine_features_masked
 faint_features_masked = np.ma.masked_less_equal(faint_features, 0)
@@ -512,17 +641,19 @@
 
 # add to grid object
 new_dict = {
-    "features_faint_strong" : {
+    "features_faint_strong": {
         "data": features_faint_strong_masked,
         "standard_name": "features_faint_strong",
         "long name": "Faint and Strong Features",
         "valid_min": 0,
         "valid_max": 10,
-        "comment_1": "3: Faint features, 2: Strong features, 1: background"
+        "comment_1": "3: Faint features, 2: Strong features, 1: background",
     }
 }
 # add field
-grid.add_field("features_faint_strong", new_dict["features_faint_strong"], replace_existing=True)
+grid.add_field(
+    "features_faint_strong", new_dict["features_faint_strong"], replace_existing=True
+)
 
 # create plot using GridMapDisplay
 # plot variables
@@ -533,9 +664,15 @@
     central_longitude=radar.longitude["data"][0],
 )
 
-faint_strong_cmap = mcolors.LinearSegmentedColormap.from_list('faint_strong',
-                                                [(32 / 255, 144 / 255, 140 / 255), (254 / 255, 238 / 255, 93 / 255),
-                                                 (254 / 255, 152 / 255, 93 / 255)], N=3)
+faint_strong_cmap = mcolors.LinearSegmentedColormap.from_list(
+    "faint_strong",
+    [
+        (32 / 255, 144 / 255, 140 / 255),
+        (254 / 255, 238 / 255, 93 / 255),
+        (254 / 255, 152 / 255, 93 / 255)
+    ],
+    N=3,
+)
 
 # plot
 plt.figure(figsize=(10, 4))
@@ -561,6 +698,65 @@
     ticklabs=["Background", "Strong Feature", "Faint Feature"],
     ax=ax2,
 )
+plt.show()
+
+######################################
+# **Image muted Strong and Faint features**
+#
+# 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)
+    muted_field = np.ma.masked_where(muted_ref.mask, field)
+
+    return nonmuted_field, muted_field
+
+# get fields
+snow_rate = grid.fields["snow_rate"]["data"][0, :, :]
+muted_ref = grid.fields["muted_reflectivity"]["data"][0, :, :]
+features = grid.fields["features_faint_strong"]["data"][0, :, :]
+x = grid.x["data"]
+y = grid.y["data"]
+
+# mute
+snow_rate_nonmuted, snow_rate_muted = quick_image_mute(snow_rate, muted_ref)
+features_nonmuted, features_muted = quick_image_mute(features, muted_ref)
+
+# muted colormap
+grays_cmap = plt.get_cmap("gray_r")
+muted_cmap = mcolors.LinearSegmentedColormap.from_list(
+    "muted_cmap", grays_cmap(np.linspace(0.2, 0.8, grays_cmap.N)), 3
+)
+
+
+# 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
+)
+srpmm = axs[0].pcolormesh(
+    x / 1000, y / 1000, snow_rate_muted, vmin=0, vmax=10, cmap=muted_cmap
+)
+axs[0].set_aspect("equal")
+axs[0].set_title("Snow rate [mm hr$^{-1}$]")
+plt.colorbar(srpm, ax=axs[0])
+
+# features
+csbpm = axs[1].pcolormesh(
+    x / 1000, y / 1000, features_nonmuted, vmin=1, vmax=3, cmap=faint_strong_cmap
+)
+csbpmm = axs[1].pcolormesh(
+    x / 1000, y / 1000, features_muted, vmin=1, vmax=3, cmap=muted_cmap
+)
+
+axs[1].set_aspect("equal")
+axs[1].set_title("Feature detection")
+cb = plt.colorbar(csbpm, ax=axs[1], ticks=[1.33, 2, 2.66])
+cb.ax.set_yticklabels(["Background", "Strong Feature", "Faint Feature"])
+
 plt.show()
 ######################################
 # Summary of recommendations and best practices
@@ -589,3 +785,8 @@
 # and M. Xu, 2003: Snow Nowcasting Using a Real-Time Correlation of Radar
 # Reflectivity with Snow Gauge Accumulation. J. Appl. Meteorol. Climatol., 42, 20–36.
 # https://doi.org/10.1175/1520-0450(2003)042%3C0020:SNUART%3E2.0.CO;2
+#
+# Tomkins, L. M., Yuter, S. E., Miller, M. A., and Allen, L. R., 2022:
+# Image muting of mixed precipitation to improve identification of regions
+# of heavy snow in radar data. Atmos. Meas. Tech., 15, 5515–5525,
+# https://doi.org/10.5194/amt-15-5515-2022

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 12/21] 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
 )

From 68987521ea9cbc19462f42474894e0dbb39eb5b2 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Mon, 13 Nov 2023 15:33:05 -0500
Subject: [PATCH 13/21] Update test_echo_class.py

---
 tests/retrieve/test_echo_class.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/tests/retrieve/test_echo_class.py b/tests/retrieve/test_echo_class.py
index ff11f74dd3..288ca32dd9 100644
--- a/tests/retrieve/test_echo_class.py
+++ b/tests/retrieve/test_echo_class.py
@@ -67,15 +67,15 @@ def test_steiner_conv_strat_modify_area(area_relation):
     assert eclass["data"].max() == 2
 
 
-def test_conv_strat__yuter_default():
+def test_feature_detection_default():
     grid = pyart.testing.make_storm_grid()
-    dict = pyart.retrieve.conv_strat_yuter(grid, bkg_rad_km=50)
+    dict = pyart.retrieve.feature_detection(grid, bkg_rad_km=50)
 
-    assert "convsf" in dict.keys()
-    assert "convsf_under" in dict.keys()
-    assert "convsf_over" in dict.keys()
+    assert "feature_detection" in dict.keys()
+    assert "feature_under" in dict.keys()
+    assert "feature_over" in dict.keys()
     assert np.all(
-        dict["convsf"]["data"][25]
+        dict["feature_detection"]["data"][0, 25, :]
         == np.array(
             [
                 0,
@@ -123,13 +123,13 @@ def test_conv_strat__yuter_default():
     )
 
 
-def test_conv_strat_yuter_noest():
+def test_feature_detection_noest():
     grid = pyart.testing.make_storm_grid()
-    dict = pyart.retrieve.conv_strat_yuter(grid, bkg_rad_km=50, estimate_flag=False)
+    dict = pyart.retrieve.feature_detection(grid, bkg_rad_km=50, estimate_flag=False)
 
-    assert "convsf" in dict.keys()
-    assert "convsf_under" not in dict.keys()
-    assert "convsf_over" not in dict.keys()
+    assert "feature_detection" in dict.keys()
+    assert "feature_under" not in dict.keys()
+    assert "feature_over" not in dict.keys()
 
 
 def test_hydroclass_semisupervised():

From ff37542c2c5520b7bfa64c6583b895e06c0fef1f Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Tue, 14 Nov 2023 09:53:33 -0500
Subject: [PATCH 14/21] Keep original function and add warning

---
 pyart/retrieve/echo_class.py | 162 +++++++++++++++++++++++++++++++++++
 1 file changed, 162 insertions(+)

diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index a52f1a2063..c312cdeeff 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -125,6 +125,168 @@ def steiner_conv_strat(
     }
 
 
+def conv_strat_yuter(
+        grid,
+        dx=None,
+        dy=None,
+        level_m=None,
+        always_core_thres=42,
+        bkg_rad_km=11,
+        use_cosine=True,
+        max_diff=5,
+        zero_diff_cos_val=55,
+        scalar_diff=1.5,
+        use_addition=True,
+        calc_thres=0.75,
+        weak_echo_thres=5.0,
+        min_dBZ_used=5.0,
+        dB_averaging=True,
+        remove_small_objects=True,
+        min_km2_size=10,
+        val_for_max_conv_rad=30,
+        max_conv_rad_km=5.0,
+        cs_core=3,
+        nosfcecho=0,
+        weakecho=3,
+        sf=1,
+        conv=2,
+        refl_field=None,
+        estimate_flag=True,
+        estimate_offset=5,
+):
+    """
+    Partition reflectivity into convective-stratiform using the Yuter et al. (2005)
+    and Yuter and Houze (1997) algorithm.
+
+    Parameters
+    ----------
+    grid : Grid
+        Grid containing reflectivity field to partition.
+    dx, dy : float, optional
+        The x- and y-dimension resolutions in meters, respectively. If None
+        the resolution is determined from the first two axes values parsed from grid object.
+    level_m : float, optional
+        Desired height in meters to classify with convective stratiform algorithm.
+    always_core_thres : float, optional
+        Threshold for points that are always convective. All values above the threshold are classifed as convective
+        See Yuter et al. (2005) for more detail.
+    bkg_rad_km : float, optional
+        Radius to compute background reflectivity in kilometers. Default is 11 km. Recommended to be at least 3 x
+        grid spacing
+    use_cosine : bool, optional
+        Boolean used to determine if a cosine scheme (see Yuter and Houze (1997)) should be used for identifying
+        convective cores (True) or if a simpler scalar scheme (False) should be used.
+    max_diff : float, optional
+        Maximum difference between background average and reflectivity in order to be classified as convective.
+        "a" value in Eqn. B1 in Yuter and Houze (1997)
+    zero_diff_cos_val : float, optional
+        Value where difference between background average and reflectivity is zero in the cosine function
+        "b" value in Eqn. B1 in Yuter and Houze (1997)
+    scalar_diff : float, optional
+        If using a scalar difference scheme, this value is the multiplier or addition to the background average
+    use_addition : bool, optional
+        Determines if a multiplier (False) or addition (True) in the scalar difference scheme should be used
+    calc_thres : float, optional
+        Minimum percentage of points needed to be considered in background average calculation
+    weak_echo_thres : float, optional
+        Threshold for determining weak echo. All values below this threshold will be considered weak echo
+    min_dBZ_used : float, optional
+        Minimum dBZ value used for classification. All values below this threshold will be considered no surface echo
+        See Yuter and Houze (1997) and Yuter et al. (2005) for more detail.
+    dB_averaging : bool, optional
+        True if using dBZ reflectivity values that need to be converted to linear Z before averaging. False for
+        other non-dBZ values (i.e. snow rate)
+    remove_small_objects : bool, optional
+        Determines if small objects should be removed from convective core array. Default is True.
+    min_km2_size : float, optional
+        Minimum size of convective cores to be considered. Cores less than this size will be removed. Default is 10
+        km^2.
+    val_for_max_conv_rad : float, optional
+        dBZ for maximum convective radius. Convective cores with values above this will have the maximum convective
+        radius
+    max_conv_rad_km : float, optional
+        Maximum radius around convective cores to classify as convective. Default is 5 km
+    cs_core : int, optional
+        Value for points classified as convective cores
+    nosfcecho : int, optional
+        Value for points classified as no surface echo, based on min_dBZ_used
+    weakecho : int, optional
+        Value for points classified as weak echo, based on weak_echo_thres
+    sf : int, optional
+        Value for points classified as stratiform
+    conv : int, optional
+        Value for points classified as convective
+    refl_field : str, optional
+        Field in grid to use as the reflectivity during partitioning. None will use the default reflectivity
+        field name from the Py-ART configuration file.
+    estimate_flag : bool, optional
+        Determines if over/underestimation should be applied. If true, the algorithm will also be run on the same field
+        wih the estimate_offset added and the same field with the estimate_offset subtracted.
+        Default is True (recommended)
+    estimate_offset : float, optional
+        Value used to offset the reflectivity values by for the over/underestimation application. Default value is 5
+        dBZ.
+
+    Returns
+    -------
+    convsf_dict : dict
+        Convective-stratiform classification dictionary.
+
+    References
+    ----------
+    Yuter, S. E., and R. A. Houze, Jr., 1997: Measurements of raindrop size
+    distributions over the Pacific warm pool and implications for Z-R relations.
+    J. Appl. Meteor., 36, 847-867.
+    https://doi.org/10.1175/1520-0450(1997)036%3C0847:MORSDO%3E2.0.CO;2
+
+    Yuter, S. E., R. A. Houze, Jr., E. A. Smith, T. T. Wilheit, and E. Zipser,
+    2005: Physical characterization of tropical oceanic convection observed in
+    KWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1
+
+    """
+
+    warn(
+        "This function will be deprecated in Py-ART 2.0."
+        " Please use feature_detection function.",
+        DeprecationWarning,
+    )
+
+    feature_dict = feature_detection(
+        grid,
+        dx=dx,
+        dy=dy,
+        level_m=level_m,
+        always_core_thres=always_core_thres,
+        bkg_rad_km=bkg_rad_km,
+        use_cosine=use_cosine,
+        max_diff=max_diff,
+        zero_diff_cos_val=zero_diff_cos_val,
+        scalar_diff=scalar_diff,
+        use_addition=use_addition,
+        calc_thres=calc_thres,
+        weak_echo_thres=weak_echo_thres,
+        min_val_used=min_dBZ_used,
+        dB_averaging=dB_averaging,
+        remove_small_objects=remove_small_objects,
+        min_km2_size=min_km2_size,
+        binary_close=False,
+        val_for_max_rad=val_for_max_conv_rad,
+        max_rad_km=max_conv_rad_km,
+        core_val=cs_core,
+        nosfcecho=nosfcecho,
+        weakecho=weakecho,
+        bkgd_val=sf,
+        feat_val=conv,
+        field=refl_field,
+        estimate_flag=estimate_flag,
+        estimate_offset=5,
+        overest_field=None,
+        underest_field=None,
+    )
+
+    return feature_dict
+
+
 def feature_detection(
     grid,
     dx=None,

From 4dde69e1267b95e59a2a3be010b1a8e78e9544f1 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Tue, 14 Nov 2023 09:54:06 -0500
Subject: [PATCH 15/21] Revert "Delete plot_convective_stratiform.py"

This reverts commit 2bff81ac0d8c33b3a51c4d6485d43994d52956ad.
---
 .../retrieve/plot_convective_stratiform.py    | 467 ++++++++++++++++++
 1 file changed, 467 insertions(+)
 create mode 100644 examples/retrieve/plot_convective_stratiform.py

diff --git a/examples/retrieve/plot_convective_stratiform.py b/examples/retrieve/plot_convective_stratiform.py
new file mode 100644
index 0000000000..da0f22cc9f
--- /dev/null
+++ b/examples/retrieve/plot_convective_stratiform.py
@@ -0,0 +1,467 @@
+"""
+=======================================
+Convective-Stratiform classification
+=======================================
+This example shows how to use the updated convective stratiform classifcation algorithm. We show 3 examples,
+a summer convective example, an example from Hurricane Ian, and an example from a winter storm.
+"""
+
+print(__doc__)
+
+# Author: Laura Tomkins (lmtomkin@ncsu.edu)
+# License: BSD 3 clause
+
+
+import cartopy.crs as ccrs
+import matplotlib.colors as mcolors
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pyart
+
+######################################
+# How the algorithm works
+# ----------
+# This first section describes how the convective-stratiform algorithm works (see references for full details). This
+# algorithm is a feature detection algorithm and classifies fields as "convective" or "stratiform". The algorithm is
+# designed to detect features in a reflectivity field but can also detect features in fields such as rain rate or
+# snow rate. In this section we describe the steps of the convective stratiform algorithm and the variables used in
+# the function.
+# The first step of the algorithm calculates a background average of the field with a circular footprint using a radius
+# provided by ``bkg_rad_km``. A larger radius will yield a smoother field. The radius needs to be at least double the
+# grid spacing, but we recommend at least three times the grid spacing. If using reflectivity, ``dB_averaging`` should be set
+# to True to convert reflectivity to linear Z before averaging, False for rescaled fields such as rain or snow rate.
+# ``calc_thres`` determines the minimum fraction of a circle that is considered in the background average calculation
+# (default is 0.75, so the points along the edges where there is less than 75% of a full circle of data,
+# the algorithm is not run).
+# Once the background average has been calculated, the original field is compared to the background average.  In
+# order for points to be considered "convective cores" they must exceed the background value by a certain value or
+# simply be greater than the ``always_core_thres``. This value is determined by either a cosine scheme, or a scalar
+# value (i.e. the reflectivity value must be X times the background value (multiplier; ``use_addition=False``),
+# or  X greater than the background value (``use_addition=True``) where X is the ``scalar_diff``).
+# ``use_cosine`` determines if a cosine scheme or scalar scheme is to be used. If ``use_cosine`` is True,
+# then the ``max_diff`` and ``zero_diff_cos_val`` come into use. These values define the cosine scheme that is used  to
+# determine the minimum difference between the background and reflectivity value in order for a core to be
+# identified. ``max_diff`` is the maximum difference between the field and the background for a core to be identified,
+# or where the cosine function crosses the y-axis. The ``zero_diff_cos_val`` is where the difference between the  field
+# and the background is zero, or where the cosine function crosses the x-axis. Note, if
+# ``always_core_thres`` < ``zero_diff_cos_val``, ``zero_diff_cos_val`` only helps define the shape of the cosine curve and
+# all values greater than ``always_core_thres`` will be considered a convective core. If
+# ``always_core_thres`` > ``zero_diff_cos_val`` then all values greater than ``zero_diff_cos_val`` will be considered a
+# convective core. We plot some examples of the schemes below:
+
+######################################
+# Example of the cosine scheme:
+pyart.graph.plot_convstrat_scheme(
+    always_core_thres=30, use_cosine=True, max_diff=5, zero_diff_cos_val=45
+)
+
+######################################
+# when zero_diff_cos_val is greater than always_core_thres, the difference becomes zero at the zero_diff_cos_val
+pyart.graph.plot_convstrat_scheme(
+    always_core_thres=55, use_cosine=True, max_diff=5, zero_diff_cos_val=45
+)
+
+######################################
+# alternatively, we can use a simpler scalar difference instead of a cosine scheme
+pyart.graph.plot_convstrat_scheme(
+    always_core_thres=40,
+    use_cosine=False,
+    max_diff=None,
+    zero_diff_cos_val=None,
+    use_addition=True,
+    scalar_diff=2,
+)
+
+######################################
+# if you are interested in picking up weak features, you can also use the scalar difference as a multiplier instead,
+# so very weak features do not have to be that different from the background to be classified as convective.
+pyart.graph.plot_convstrat_scheme(
+    always_core_thres=40,
+    use_cosine=False,
+    max_diff=None,
+    zero_diff_cos_val=None,
+    use_addition=False,
+    scalar_diff=2,
+)
+
+######################################
+# Once the cores are identified, there is an option to remove speckles (``remove_small_objects``) smaller than a  given
+# size (``min_km2_size``).
+# After the convective cores are identified, We then incorporate convective radii using
+# ``val_for_max_conv_rad`` and ``max_conv_rad_km``. The convective radii act as a dilation and are used to classify
+# additional points around the cores as convective that may not have been identified previously.  The
+# ``val_for_max_conv_rad`` is the value where the maximum convective radius is applied and the ``max_conv_rad_km`` is the
+# maximum convective radius. Values less than the ``val_for_max_conv_rad`` are assigned a convective radius using a step
+# function.
+# Finally, the points are classified as NOSFCECHO (threshold set with ``min_dBZ_used``; 0), WEAKECHO (threshold set with
+# ``weak_echo_thres``; 3), SF (stratiform; 1), CONV (convective; 2).
+
+######################################
+# Examples
+# ----------
+# **Classification of summer convective example**
+#
+# Our first example classifies echo from a summer convective event. We use a cosine scheme to classify the convective
+# points.
+
+# Now let's do a classification with our parameters
+# read in file
+filename = pyart.testing.get_test_data("swx_20120520_0641.nc")
+radar = pyart.io.read(filename)
+
+# extract the lowest sweep
+radar = radar.extract_sweeps([0])
+
+# interpolate to grid
+grid = pyart.map.grid_from_radars(
+    (radar,),
+    grid_shape=(1, 201, 201),
+    grid_limits=((0, 10000), (-50000.0, 50000.0), (-50000.0, 50000.0)),
+    fields=["reflectivity_horizontal"],
+)
+
+# get dx dy
+dx = grid.x["data"][1] - grid.x["data"][0]
+dy = grid.y["data"][1] - grid.y["data"][0]
+
+# convective stratiform classification
+convsf_dict = pyart.retrieve.conv_strat_yuter(
+    grid,
+    dx,
+    dy,
+    refl_field="reflectivity_horizontal",
+    always_core_thres=40,
+    bkg_rad_km=20,
+    use_cosine=True,
+    max_diff=5,
+    zero_diff_cos_val=55,
+    weak_echo_thres=10,
+    max_conv_rad_km=2,
+)
+
+# add to grid object
+# mask zero values (no surface echo)
+convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
+# mask 3 values (weak echo)
+convsf_masked = np.ma.masked_equal(convsf_masked, 3)
+# add dimension to array to add to grid object
+convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
+# add field
+grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
+
+# create plot using GridMapDisplay
+# plot variables
+display = pyart.graph.GridMapDisplay(grid)
+magma_r_cmap = plt.get_cmap("magma_r")
+ref_cmap = mcolors.LinearSegmentedColormap.from_list(
+    "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))
+)
+projection = ccrs.AlbersEqualArea(
+    central_latitude=radar.latitude["data"][0],
+    central_longitude=radar.longitude["data"][0],
+)
+
+# plot
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(1, 2, 1, projection=projection)
+display.plot_grid(
+    "reflectivity_horizontal",
+    vmin=5,
+    vmax=45,
+    cmap=ref_cmap,
+    transform=ccrs.PlateCarree(),
+    ax=ax1,
+)
+ax2 = plt.subplot(1, 2, 2, projection=projection)
+display.plot_grid(
+    "convsf",
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
+    ax=ax2,
+    transform=ccrs.PlateCarree(),
+    ticks=[1 / 3, 1, 5 / 3],
+    ticklabs=["", "Stratiform", "Convective"],
+)
+plt.show()
+
+
+######################################
+# In addition to the default convective-stratiform classification, the function also returns an underestimate
+# (convsf_under) and an overestimate (convsf_over) to take into consideration the uncertainty when choosing
+# classification parameters. The under and overestimate use the same parameters, but vary the input field by a
+# certain value (default is 5 dBZ, can be changed with ``estimate_offset``). The estimation can be turned off (
+# ``estimate_flag=False``), but we recommend keeping it turned on.
+
+# mask weak echo and no surface echo
+convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_masked, 3)
+convsf_dict["convsf"]["data"] = convsf_masked
+# underest.
+convsf_masked = np.ma.masked_equal(convsf_dict["convsf_under"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_masked, 3)
+convsf_dict["convsf_under"]["data"] = convsf_masked
+# overest.
+convsf_masked = np.ma.masked_equal(convsf_dict["convsf_over"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_masked, 3)
+convsf_dict["convsf_over"]["data"] = convsf_masked
+
+# Plot each estimation
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(131)
+ax1.pcolormesh(
+    convsf_dict["convsf"]["data"][0, :, :],
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
+)
+ax1.set_title("Best estimate")
+ax1.set_aspect("equal")
+ax2 = plt.subplot(132)
+ax2.pcolormesh(
+    convsf_dict["convsf_under"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
+)
+ax2.set_title("Underestimate")
+ax2.set_aspect("equal")
+ax3 = plt.subplot(133)
+ax3.pcolormesh(
+    convsf_dict["convsf_over"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
+)
+ax3.set_title("Overestimate")
+ax3.set_aspect("equal")
+plt.show()
+
+######################################
+# **Tropical example**
+#
+# Let's get a NEXRAD file from Hurricane Ian
+
+# Read in file
+nexrad_file = "s3://noaa-nexrad-level2/2022/09/28/KTBW/KTBW20220928_190142_V06"
+radar = pyart.io.read_nexrad_archive(nexrad_file)
+
+# extract the lowest sweep
+radar = radar.extract_sweeps([0])
+
+# interpolate to grid
+grid = pyart.map.grid_from_radars(
+    (radar,),
+    grid_shape=(1, 201, 201),
+    grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)),
+    fields=["reflectivity"],
+)
+
+# get dx dy
+dx = grid.x["data"][1] - grid.x["data"][0]
+dy = grid.y["data"][1] - grid.y["data"][0]
+
+# convective stratiform classification
+convsf_dict = pyart.retrieve.conv_strat_yuter(
+    grid,
+    dx,
+    dy,
+    refl_field="reflectivity",
+    always_core_thres=40,
+    bkg_rad_km=20,
+    use_cosine=True,
+    max_diff=3,
+    zero_diff_cos_val=55,
+    weak_echo_thres=5,
+    max_conv_rad_km=2,
+    estimate_flag=False,
+)
+
+# add to grid object
+# mask zero values (no surface echo)
+convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
+# mask 3 values (weak echo)
+convsf_masked = np.ma.masked_equal(convsf_masked, 3)
+# add dimension to array to add to grid object
+convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
+# add field
+grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
+
+# create plot using GridMapDisplay
+# plot variables
+display = pyart.graph.GridMapDisplay(grid)
+magma_r_cmap = plt.get_cmap("magma_r")
+ref_cmap = mcolors.LinearSegmentedColormap.from_list(
+    "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))
+)
+projection = ccrs.AlbersEqualArea(
+    central_latitude=radar.latitude["data"][0],
+    central_longitude=radar.longitude["data"][0],
+)
+# plot
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(1, 2, 1, projection=projection)
+display.plot_grid(
+    "reflectivity",
+    vmin=5,
+    vmax=45,
+    cmap=ref_cmap,
+    transform=ccrs.PlateCarree(),
+    ax=ax1,
+    axislabels_flag=False,
+)
+ax2 = plt.subplot(1, 2, 2, projection=projection)
+display.plot_grid(
+    "convsf",
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
+    axislabels_flag=False,
+    transform=ccrs.PlateCarree(),
+    ticks=[1 / 3, 1, 5 / 3],
+    ticklabs=["", "Stratiform", "Convective"],
+    ax=ax2,
+)
+plt.show()
+
+######################################
+# **Winter storm example with image muting**
+#
+# Here is a final example of the convective stratiform classification using an example from a winter storm. Before
+# doing the classification, we image mute the reflectivity to remove regions with melting or mixed precipitation. We
+# then rescale the reflectivity to snow rate (Rasumussen et al. 2003). We recommend using a rescaled reflectivity
+# to do the classification, but if you do make sure to changed dB_averaging to False because this parameter is  used
+# to convert reflectivity to a linear value before averaging (set dB_averaging to True for reflectivity fields in
+# dBZ units).
+# In this example, note how we change some of the other parameters since we are classifying snow rate instead of
+# reflecitivity.
+
+# Read in file
+nexrad_file = "s3://noaa-nexrad-level2/2021/02/07/KOKX/KOKX20210207_161413_V06"
+radar = pyart.io.read_nexrad_archive(nexrad_file)
+
+# extract the lowest sweep
+radar = radar.extract_sweeps([0])
+
+# interpolate to grid
+grid = pyart.map.grid_from_radars(
+    (radar,),
+    grid_shape=(1, 201, 201),
+    grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)),
+    fields=["reflectivity", "cross_correlation_ratio"],
+)
+
+# image mute grid object
+grid = pyart.util.image_mute_radar(
+    grid, "reflectivity", "cross_correlation_ratio", 0.97, 20
+)
+
+# convect non-muted reflectivity to snow rate
+nonmuted_ref = grid.fields["nonmuted_reflectivity"]["data"][0, :, :]
+nonmuted_ref = np.ma.masked_invalid(nonmuted_ref)
+
+nonmuted_ref_linear = 10 ** (nonmuted_ref / 10)  # mm6/m3
+snow_rate = (nonmuted_ref_linear / 57.3) ** (1 / 1.67)  #
+
+# add to grid
+snow_rate_dict = {
+    "data": snow_rate[None, :, :],
+    "standard_name": "snow_rate",
+    "long_name": "Snow rate converted from linear reflectivity",
+    "units": "mm/hr",
+    "valid_min": 0,
+    "valid_max": 40500,
+}
+grid.add_field("snow_rate", snow_rate_dict, replace_existing=True)
+
+# get dx dy
+dx = grid.x["data"][1] - grid.x["data"][0]
+dy = grid.y["data"][1] - grid.y["data"][0]
+
+# convective stratiform classification
+convsf_dict = pyart.retrieve.conv_strat_yuter(
+    grid,
+    dx,
+    dy,
+    refl_field="snow_rate",
+    dB_averaging=False,
+    always_core_thres=4,
+    bkg_rad_km=40,
+    use_cosine=True,
+    max_diff=1.5,
+    zero_diff_cos_val=5,
+    weak_echo_thres=0,
+    min_dBZ_used=0,
+    max_conv_rad_km=1,
+    estimate_flag=False,
+)
+
+# add to grid object
+# mask zero values (no surface echo)
+convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
+# mask 3 values (weak echo)
+convsf_masked = np.ma.masked_equal(convsf_masked, 3)
+# add dimension to array to add to grid object
+convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
+# add field
+grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
+
+# create plot using GridMapDisplay
+# plot variables
+display = pyart.graph.GridMapDisplay(grid)
+magma_r_cmap = plt.get_cmap("magma_r")
+ref_cmap = mcolors.LinearSegmentedColormap.from_list(
+    "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))
+)
+projection = ccrs.AlbersEqualArea(
+    central_latitude=radar.latitude["data"][0],
+    central_longitude=radar.longitude["data"][0],
+)
+# plot
+plt.figure(figsize=(10, 4))
+ax1 = plt.subplot(1, 2, 1, projection=projection)
+display.plot_grid(
+    "snow_rate",
+    vmin=0,
+    vmax=10,
+    cmap=plt.get_cmap("viridis"),
+    transform=ccrs.PlateCarree(),
+    ax=ax1,
+    axislabels_flag=False,
+)
+ax2 = plt.subplot(1, 2, 2, projection=projection)
+display.plot_grid(
+    "convsf",
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
+    axislabels_flag=False,
+    transform=ccrs.PlateCarree(),
+    ticks=[1 / 3, 1, 5 / 3],
+    ticklabs=["", "Stratiform", "Convective"],
+    ax=ax2,
+)
+plt.show()
+
+######################################
+# Summary of recommendations and best practices
+# ----------
+# * Tune your parameters to your specific purpose
+# * Use a rescaled field if possible (i.e. linear reflectivity, rain or snow rate)
+# * Keep ``estimate_flag=True`` to see uncertainty in classification
+#
+# References
+# ----------
+# Steiner, M. R., R. A. Houze Jr., and S. E. Yuter, 1995: Climatological
+# Characterization of Three-Dimensional Storm Structure from Operational
+# Radar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007.
+# https://doi.org/10.1175/1520-0450(1995)034<1978:CCOTDS>2.0.CO;2.
+#
+# Yuter, S. E., and R. A. Houze, Jr., 1997: Measurements of raindrop size
+# distributions over the Pacific warm pool and implications for Z-R relations.
+# J. Appl. Meteor., 36, 847-867.
+# https://doi.org/10.1175/1520-0450(1997)036%3C0847:MORSDO%3E2.0.CO;2
+#
+# Yuter, S. E., R. A. Houze, Jr., E. A. Smith, T. T. Wilheit, and E. Zipser,
+# 2005: Physical characterization of tropical oceanic convection observed in
+# KWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1
+#
+# Rasmussen, R., M. Dixon, S. Vasiloff, F. Hage, S. Knight, J. Vivekanandan,
+# and M. Xu, 2003: Snow Nowcasting Using a Real-Time Correlation of Radar
+# Reflectivity with Snow Gauge Accumulation. J. Appl. Meteorol. Climatol., 42, 20–36.
+# https://doi.org/10.1175/1520-0450(2003)042%3C0020:SNUART%3E2.0.CO;2

From befb2c77858843c6293ca6089f63f9c42d88df5e Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Tue, 14 Nov 2023 09:54:45 -0500
Subject: [PATCH 16/21] Update __init__.py

---
 pyart/retrieve/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py
index 685059ce4f..3b79ba5391 100644
--- a/pyart/retrieve/__init__.py
+++ b/pyart/retrieve/__init__.py
@@ -6,6 +6,7 @@
 from .advection import grid_displacement_pc, grid_shift  # noqa
 from .comp_z import composite_reflectivity  # noqa
 from .echo_class import feature_detection  # noqa
+from .echo_class import conv_strat_yuter  # noqa
 from .echo_class import get_freq_band  # noqa
 from .echo_class import hydroclass_semisupervised  # noqa
 from .echo_class import steiner_conv_strat  # noqa

From 12d6022e11d8ecc034f5057bb542b79904a83c41 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Tue, 14 Nov 2023 09:56:21 -0500
Subject: [PATCH 17/21] Update echo_class.py

---
 pyart/retrieve/echo_class.py | 54 ++++++++++++++++++------------------
 1 file changed, 27 insertions(+), 27 deletions(-)

diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py
index c312cdeeff..7a15fa448c 100644
--- a/pyart/retrieve/echo_class.py
+++ b/pyart/retrieve/echo_class.py
@@ -126,33 +126,33 @@ def steiner_conv_strat(
 
 
 def conv_strat_yuter(
-        grid,
-        dx=None,
-        dy=None,
-        level_m=None,
-        always_core_thres=42,
-        bkg_rad_km=11,
-        use_cosine=True,
-        max_diff=5,
-        zero_diff_cos_val=55,
-        scalar_diff=1.5,
-        use_addition=True,
-        calc_thres=0.75,
-        weak_echo_thres=5.0,
-        min_dBZ_used=5.0,
-        dB_averaging=True,
-        remove_small_objects=True,
-        min_km2_size=10,
-        val_for_max_conv_rad=30,
-        max_conv_rad_km=5.0,
-        cs_core=3,
-        nosfcecho=0,
-        weakecho=3,
-        sf=1,
-        conv=2,
-        refl_field=None,
-        estimate_flag=True,
-        estimate_offset=5,
+    grid,
+    dx=None,
+    dy=None,
+    level_m=None,
+    always_core_thres=42,
+    bkg_rad_km=11,
+    use_cosine=True,
+    max_diff=5,
+    zero_diff_cos_val=55,
+    scalar_diff=1.5,
+    use_addition=True,
+    calc_thres=0.75,
+    weak_echo_thres=5.0,
+    min_dBZ_used=5.0,
+    dB_averaging=True,
+    remove_small_objects=True,
+    min_km2_size=10,
+    val_for_max_conv_rad=30,
+    max_conv_rad_km=5.0,
+    cs_core=3,
+    nosfcecho=0,
+    weakecho=3,
+    sf=1,
+    conv=2,
+    refl_field=None,
+    estimate_flag=True,
+    estimate_offset=5,
 ):
     """
     Partition reflectivity into convective-stratiform using the Yuter et al. (2005)

From 2dc369c84fc3223b642c5c542282cc18f04885c6 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Tue, 14 Nov 2023 10:48:47 -0500
Subject: [PATCH 18/21] Update dictionary keys

---
 .../retrieve/plot_convective_stratiform.py    | 36 +++++++++----------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/examples/retrieve/plot_convective_stratiform.py b/examples/retrieve/plot_convective_stratiform.py
index da0f22cc9f..8623b9ef5c 100644
--- a/examples/retrieve/plot_convective_stratiform.py
+++ b/examples/retrieve/plot_convective_stratiform.py
@@ -142,13 +142,13 @@
 
 # add to grid object
 # mask zero values (no surface echo)
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_dict["feature_detection"]["data"], 0)
 # mask 3 values (weak echo)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
 # add dimension to array to add to grid object
-convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
+convsf_dict["feature_detection"]["data"] = convsf_masked[None, :, :]
 # add field
-grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
+grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True)
 
 # create plot using GridMapDisplay
 # plot variables
@@ -195,23 +195,23 @@
 # ``estimate_flag=False``), but we recommend keeping it turned on.
 
 # mask weak echo and no surface echo
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_dict["feature_detection"]["data"], 0)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-convsf_dict["convsf"]["data"] = convsf_masked
+convsf_dict["feature_detection"]["data"] = convsf_masked
 # underest.
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf_under"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_dict["feature_under"]["data"], 0)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-convsf_dict["convsf_under"]["data"] = convsf_masked
+convsf_dict["feature_under"]["data"] = convsf_masked
 # overest.
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf_over"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_dict["feature_over"]["data"], 0)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
-convsf_dict["convsf_over"]["data"] = convsf_masked
+convsf_dict["feature_over"]["data"] = convsf_masked
 
 # Plot each estimation
 plt.figure(figsize=(10, 4))
 ax1 = plt.subplot(131)
 ax1.pcolormesh(
-    convsf_dict["convsf"]["data"][0, :, :],
+    convsf_dict["feature_detection"]["data"][0, :, :],
     vmin=0,
     vmax=2,
     cmap=plt.get_cmap("viridis", 3),
@@ -220,13 +220,13 @@
 ax1.set_aspect("equal")
 ax2 = plt.subplot(132)
 ax2.pcolormesh(
-    convsf_dict["convsf_under"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
+    convsf_dict["feature_under"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
 )
 ax2.set_title("Underestimate")
 ax2.set_aspect("equal")
 ax3 = plt.subplot(133)
 ax3.pcolormesh(
-    convsf_dict["convsf_over"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
+    convsf_dict["feature_over"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
 )
 ax3.set_title("Overestimate")
 ax3.set_aspect("equal")
@@ -274,13 +274,13 @@
 
 # add to grid object
 # mask zero values (no surface echo)
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_dict["feature_detection"]["data"], 0)
 # mask 3 values (weak echo)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
 # add dimension to array to add to grid object
-convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
+convsf_dict["feature_detection"]["data"] = convsf_masked[None, :, :]
 # add field
-grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
+grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True)
 
 # create plot using GridMapDisplay
 # plot variables
@@ -393,13 +393,13 @@
 
 # add to grid object
 # mask zero values (no surface echo)
-convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0)
+convsf_masked = np.ma.masked_equal(convsf_dict["feature_detection"]["data"], 0)
 # mask 3 values (weak echo)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
 # add dimension to array to add to grid object
-convsf_dict["convsf"]["data"] = convsf_masked[None, :, :]
+convsf_dict["feature_detection"]["data"] = convsf_masked[None, :, :]
 # add field
-grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True)
+grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True)
 
 # create plot using GridMapDisplay
 # plot variables

From d2d99ac10ef176f4f05b3db20b60c95523c18eb6 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Tue, 14 Nov 2023 10:50:40 -0500
Subject: [PATCH 19/21] Update plot_convective_stratiform.py

---
 examples/retrieve/plot_convective_stratiform.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/examples/retrieve/plot_convective_stratiform.py b/examples/retrieve/plot_convective_stratiform.py
index 8623b9ef5c..4d02dbedbc 100644
--- a/examples/retrieve/plot_convective_stratiform.py
+++ b/examples/retrieve/plot_convective_stratiform.py
@@ -220,7 +220,10 @@
 ax1.set_aspect("equal")
 ax2 = plt.subplot(132)
 ax2.pcolormesh(
-    convsf_dict["feature_under"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
+    convsf_dict["feature_under"]["data"],
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3),
 )
 ax2.set_title("Underestimate")
 ax2.set_aspect("equal")

From e4414357c989723e4e148592666b5fd60fcdb268 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Tue, 14 Nov 2023 11:07:39 -0500
Subject: [PATCH 20/21] Fix issues

---
 examples/retrieve/plot_convective_stratiform.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/examples/retrieve/plot_convective_stratiform.py b/examples/retrieve/plot_convective_stratiform.py
index 4d02dbedbc..dbb21b3b04 100644
--- a/examples/retrieve/plot_convective_stratiform.py
+++ b/examples/retrieve/plot_convective_stratiform.py
@@ -146,7 +146,7 @@
 # mask 3 values (weak echo)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
 # add dimension to array to add to grid object
-convsf_dict["feature_detection"]["data"] = convsf_masked[None, :, :]
+convsf_dict["feature_detection"]["data"] = convsf_masked
 # add field
 grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True)
 
@@ -220,7 +220,7 @@
 ax1.set_aspect("equal")
 ax2 = plt.subplot(132)
 ax2.pcolormesh(
-    convsf_dict["feature_under"]["data"],
+    convsf_dict["feature_under"]["data"][0, :, :],
     vmin=0,
     vmax=2,
     cmap=plt.get_cmap("viridis", 3),
@@ -229,7 +229,10 @@
 ax2.set_aspect("equal")
 ax3 = plt.subplot(133)
 ax3.pcolormesh(
-    convsf_dict["feature_over"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3)
+    convsf_dict["feature_over"]["data"][0, :, :],
+    vmin=0,
+    vmax=2,
+    cmap=plt.get_cmap("viridis", 3)
 )
 ax3.set_title("Overestimate")
 ax3.set_aspect("equal")
@@ -281,7 +284,7 @@
 # mask 3 values (weak echo)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
 # add dimension to array to add to grid object
-convsf_dict["feature_detection"]["data"] = convsf_masked[None, :, :]
+convsf_dict["feature_detection"]["data"] = convsf_masked
 # add field
 grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True)
 
@@ -400,7 +403,7 @@
 # mask 3 values (weak echo)
 convsf_masked = np.ma.masked_equal(convsf_masked, 3)
 # add dimension to array to add to grid object
-convsf_dict["feature_detection"]["data"] = convsf_masked[None, :, :]
+convsf_dict["feature_detection"]["data"] = convsf_masked
 # add field
 grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True)
 

From eb1062c5115211e9a730e788756490ff0aef3d04 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Tue, 14 Nov 2023 11:09:10 -0500
Subject: [PATCH 21/21] Update plot_convective_stratiform.py

---
 examples/retrieve/plot_convective_stratiform.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/retrieve/plot_convective_stratiform.py b/examples/retrieve/plot_convective_stratiform.py
index dbb21b3b04..2d905ee208 100644
--- a/examples/retrieve/plot_convective_stratiform.py
+++ b/examples/retrieve/plot_convective_stratiform.py
@@ -232,7 +232,7 @@
     convsf_dict["feature_over"]["data"][0, :, :],
     vmin=0,
     vmax=2,
-    cmap=plt.get_cmap("viridis", 3)
+    cmap=plt.get_cmap("viridis", 3),
 )
 ax3.set_title("Overestimate")
 ax3.set_aspect("equal")