add function to find and plot local extrema

src/metpy/calc/tools.py

@@ -781,6 +781,38 @@ def take(indexer):
     return take
 
 
+@exporter.export
+def find_local_extrema(var, nsize, extrema):
+    r"""Find the local extreme (max/min) values of an 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
+    --------
+    plot_local_extrema
+
+    """
+    from scipy.ndimage import maximum_filter, minimum_filter
+    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,67 @@ 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
+
+    See Also
+    --------
+    find_local_extrema
+
+    """
+    defaultkwargs = {'size': 20, 'color': 'black', 'fontweight': 'bold',
+                     'horizontalalignment': 'center', 'verticalalignment': 'center',
+                     'transform': None}
+    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[stack_vals.metpy.x.name].values
+        y = extrema[stack_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,14 +13,16 @@
 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,
-                              DIR_STRS, nominal_lat_lon_grid_deltas, parse_grid_arguments, UND)
+                              DIR_STRS, nominal_lat_lon_grid_deltas,
+                              parse_grid_arguments, UND)
 from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_array_equal,
                            get_test_data)
 from metpy.units import units
@@ -475,6 +477,47 @@ 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."""
+    mslp_data = local_extrema_data
+    local_max = find_local_extrema(mslp_data, 3, 'max')
+    local_min = find_local_extrema(mslp_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)
+
+
 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

@@ -5,13 +5,17 @@
 
 from datetime import datetime
 
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 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 +158,42 @@ 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.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(1, figsize=(17., 11.))
+    ax = plt.subplot(111, projection=ccrs.LambertConformal(central_latitude=35,
+                                                           central_longitude=-100))
+
+    # Set extent and plot map lines
+    ax.set_extent([-124., -70, 20., 60.], ccrs.PlateCarree())
+    ax.add_feature(cfeature.COASTLINE.with_scale('50m'),
+                edgecolor='grey', linewidth=0.75)
+    ax.add_feature(cfeature.STATES.with_scale('50m'),
+                edgecolor='grey', linewidth=0.5)
+
+    # Plot thickness with multiple colors
+    clevs = (np.arange(0, 5400, 60),
+            np.array([5400]),
+            np.arange(5460, 7000, 60))
+
+    # Plot MSLP
+    clevmslp = np.arange(800., 1120., 4)
+    cs2 = ax.contour(mslp.lon, mslp.lat, mslp.metpy.convert_units('hPa'),
+                    clevmslp, colors='k', linewidths=1.25,
+                    linestyles='solid', transform=ccrs.PlateCarree())
+    plt.clabel(cs2)
+
+    plot_local_extrema(ax, relmax2D, 'H', plot_val=False, color='tab:red',
+                       transform=ccrs.PlateCarree())
+    plot_local_extrema(ax, relmin2D, 'L', color='tab:blue', transform=ccrs.PlateCarree())
+
+    return fig