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