add function to find and plot local extrema

docs/_templates/overrides/metpy.calc.rst

@@ -205,6 +205,7 @@ Other
       azimuth_range_to_lat_lon
       find_bounding_indices
       find_intersections
+      find_local_extrema
       get_layer
       get_layer_heights
       get_perturbation

src/metpy/calc/tools.py

@@ -781,6 +781,41 @@ def take(indexer):
     return take
 
 
+@exporter.export
+def find_local_extrema(var, nsize, extrema):
+    r"""Find the local extreme (max/min) values of a 2D array.
+
+    Parameters
+    ----------
+    var : `xarray.DataArray`
+        The variable to locate the local extrema using the nearest method
+        from the maximum_filter or minimum_filter from the scipy.ndimage module.
+    nsize : int
+        The minimum number of grid points between each local extrema.
+    extrema: str
+        The value 'max' for local maxima or 'min' for local minima.
+
+    Returns
+    -------
+    var_extrema: `xarray.DataArray`
+        The values of the local extrema with other values as NaNs
+
+    See Also
+    --------
+    :func:`~metpy.plots.plot_local_extrema`
+
+    """
+    from scipy.ndimage import maximum_filter, minimum_filter
+    if extrema not in ['max', 'min']:
+        raise ValueError('Invalid input for "extrema". Valid options are "max" or "min".')
+
+    if extrema == 'max':
+        extreme_val = maximum_filter(var.values, nsize, mode='nearest')
+    elif extrema == 'min':
+        extreme_val = minimum_filter(var.values, nsize, mode='nearest')
+    return var.where(extreme_val == var.values)
+
+
 @exporter.export
 @preprocess_and_wrap()
 def lat_lon_grid_deltas(longitude, latitude, x_dim=-1, y_dim=-2, geod=None):

src/metpy/plots/__init__.py

@@ -7,7 +7,7 @@
 from . import _mpl  # noqa: F401
 from . import cartopy_utils, plot_areas
 from ._util import (add_metpy_logo, add_timestamp, add_unidata_logo,  # noqa: F401
-                    convert_gempak_color)
+                    convert_gempak_color, plot_local_extrema)
 from .ctables import *  # noqa: F403
 from .declarative import *  # noqa: F403
 from .patheffects import *  # noqa: F403
@@ -23,7 +23,7 @@
 __all__.extend(station_plot.__all__)  # pylint: disable=undefined-variable
 __all__.extend(wx_symbols.__all__)  # pylint: disable=undefined-variable
 __all__.extend(['add_metpy_logo', 'add_timestamp', 'add_unidata_logo',
-                'convert_gempak_color'])
+                'convert_gempak_color', 'plot_local_extrema'])
 
 set_module(globals())
 

src/metpy/plots/_util.py

@@ -286,3 +286,62 @@ def normalize(x):
     except TypeError:
         res = cols[normalize(c)]
     return res
+
+
+def plot_local_extrema(ax, extreme_vals, symbol, plot_val=True, **kwargs):
+    """Plot the local extreme (max/min) values of an array.
+
+    The behavior of the plotting will have the symbol horizontal/vertical alignment
+    be center/bottom and any value plotted will be center/top. The text size of plotted
+    values is 0.65 of the symbol size.
+
+    Parameters
+    ----------
+    ax : `matplotlib.axes`
+        The axes which to plot the local extrema
+    extreme_vals : `xarray.DataArray`
+        The DataArray that contains the variable local extrema
+    symbol : str
+        The text or other string to plot at the local extrema location
+    plot_val : bool
+        Whether to plot the local extreme value (default is True)
+
+    Returns
+    -------
+    Plots local extrema on the plot axes
+
+    Other Parameters
+    ----------------
+    kwargs : `matplotlib.pyplot.text` properties.
+        Other valid `matplotlib.pyplot.text` kwargs can be specified
+        except verticalalalignment if plotting both a symbol and the value.
+
+        Default kwargs:
+            size : 20
+            color : 'black'
+            fontweight : 'bold'
+            horizontalalignment : 'center'
+            verticalalignment : 'center'
+            transform : None
+
+    """
+    defaultkwargs = {'size': 20, 'color': 'black', 'fontweight': 'bold',
+                     'horizontalalignment': 'center', 'verticalalignment': 'center'}
+    kwargs = {**defaultkwargs, **kwargs}
+    if plot_val:
+        kwargs.pop('verticalalignment')
+    size = kwargs.pop('size')
+    textsize = size * .65
+
+    stack_vals = extreme_vals.stack(x=[extreme_vals.metpy.x.name, extreme_vals.metpy.y.name])
+    for extrema in stack_vals[stack_vals.notnull()]:
+        x = extrema[extreme_vals.metpy.x.name].values
+        y = extrema[extreme_vals.metpy.y.name].values
+        if plot_val:
+            ax.text(x, y, symbol, clip_on=True, clip_box=ax.bbox, size=size,
+                    verticalalignment='bottom', **kwargs)
+            ax.text(x, y, f'{extrema.values:.0f}', clip_on=True, clip_box=ax.bbox,
+                    size=textsize, verticalalignment='top', **kwargs)
+        else:
+            ax.text(x, y, symbol, clip_on=True, clip_box=ax.bbox, size=size,
+                    **kwargs)

tests/calc/test_calc_tools.py

@@ -13,10 +13,11 @@
 import xarray as xr
 
 from metpy.calc import (angle_to_direction, find_bounding_indices, find_intersections,
-                        first_derivative, geospatial_gradient, get_layer, get_layer_heights,
-                        gradient, laplacian, lat_lon_grid_deltas, nearest_intersection_idx,
-                        parse_angle, pressure_to_height_std, reduce_point_density,
-                        resample_nn_1d, second_derivative, vector_derivative)
+                        find_local_extrema, first_derivative, geospatial_gradient, get_layer,
+                        get_layer_heights, gradient, laplacian, lat_lon_grid_deltas,
+                        nearest_intersection_idx, parse_angle, pressure_to_height_std,
+                        reduce_point_density, resample_nn_1d, second_derivative,
+                        vector_derivative)
 from metpy.calc.tools import (_delete_masked_points, _get_bound_pressure_height,
                               _greater_or_close, _less_or_close, _next_non_masked_element,
                               _remove_nans, azimuth_range_to_lat_lon, BASE_DEGREE_MULTIPLIER,
@@ -475,6 +476,48 @@ def test_get_layer_heights_agl_bottom_no_interp():
     assert_array_almost_equal(data_true, data, 6)
 
 
+@pytest.fixture
+def local_extrema_data():
+    """Test data for local extrema finding."""
+    data = xr.DataArray(
+        np.array([[101628.24, 101483.67, 101366.06, 101287.55, 101233.45],
+                  [101637.19, 101515.555, 101387.164, 101280.32, 101210.15],
+                  [101581.78, 101465.234, 101342., 101233.22, 101180.25],
+                  [101404.31, 101318.4, 101233.18, 101166.445, 101159.93],
+                  [101280.586, 101238.445, 101195.234, 101183.34, 101212.8]]),
+        name='mslp',
+        dims=('lat', 'lon'),
+        coords={'lat': xr.DataArray(np.array([45., 43., 41., 39., 37.]),
+                                    dims=('lat',), attrs={'units': 'degrees_north'}),
+                'lon': xr.DataArray(np.array([265., 267., 269., 271., 273.]),
+                                    dims=('lon',), attrs={'units': 'degrees_east'})},
+        attrs={'units': 'Pa'}
+    )
+    return data
+
+
+def test_find_local_extrema(local_extrema_data):
+    """Test find_local_extrema function for maximum."""
+    local_max = find_local_extrema(local_extrema_data, 3, 'max')
+    local_min = find_local_extrema(local_extrema_data, 3, 'min')
+
+    max_truth = np.array([[np.nan, np.nan, np.nan, np.nan, np.nan],
+                          [101637.19, np.nan, np.nan, np.nan, np.nan],
+                          [np.nan, np.nan, np.nan, np.nan, np.nan],
+                          [np.nan, np.nan, np.nan, np.nan, np.nan],
+                          [np.nan, np.nan, np.nan, np.nan, 101212.8]])
+    min_truth = np.array([[np.nan, np.nan, np.nan, np.nan, np.nan],
+                          [np.nan, np.nan, np.nan, np.nan, np.nan],
+                          [np.nan, np.nan, np.nan, np.nan, np.nan],
+                          [np.nan, np.nan, np.nan, np.nan, 101159.93],
+                          [np.nan, np.nan, np.nan, np.nan, np.nan]])
+    assert_array_almost_equal(local_max.data, max_truth)
+    assert_array_almost_equal(local_min.data, min_truth)
+
+    with pytest.raises(ValueError):
+        find_local_extrema(local_extrema_data, 3, 'large')
+
+
 def test_lat_lon_grid_deltas_1d():
     """Test for lat_lon_grid_deltas for variable grid."""
     lat = np.arange(40, 50, 2.5)

tests/plots/test_util.py

@@ -11,7 +11,9 @@
 import pytest
 import xarray as xr
 
-from metpy.plots import add_metpy_logo, add_timestamp, add_unidata_logo, convert_gempak_color
+from metpy.calc import find_local_extrema
+from metpy.plots import (add_metpy_logo, add_timestamp, add_unidata_logo, convert_gempak_color,
+                         plot_local_extrema)
 from metpy.testing import get_test_data
 
 MPL_VERSION = matplotlib.__version__[:3]
@@ -154,3 +156,26 @@ def test_gempak_color_scalar():
     mplc = convert_gempak_color(6)
     truth = 'cyan'
     assert mplc == truth
+
+
+@pytest.mark.mpl_image_compare(remove_text=True)
+def test_plot_extrema():
+    """Test plotting of local max/min values."""
+    data = xr.open_dataset(get_test_data('GFS_test.nc', as_file_obj=False))
+
+    mslp = data.metpy.parse_cf('Pressure_reduced_to_MSL_msl').squeeze()
+    relmax2d = find_local_extrema(mslp, 10, 'max').metpy.convert_units('hPa')
+    relmin2d = find_local_extrema(mslp, 15, 'min').metpy.convert_units('hPa')
+
+    fig = plt.figure(figsize=(8., 8.))
+    ax = fig.add_subplot(1, 1, 1)
+
+    # Plot MSLP
+    clevmslp = np.arange(800., 1120., 4)
+    ax.contour(mslp.lon, mslp.lat, mslp.metpy.convert_units('hPa'),
+               clevmslp, colors='k', linewidths=1.25, linestyles='solid')
+
+    plot_local_extrema(ax, relmax2d, 'H', plot_val=False, color='tab:red')
+    plot_local_extrema(ax, relmin2d, 'L', color='tab:blue')
+
+    return fig