Find peaks

docs/_templates/overrides/metpy.calc.rst

@@ -215,12 +215,14 @@ Other
       azimuth_range_to_lat_lon
       find_bounding_indices
       find_intersections
+      find_peaks
       get_layer
       get_layer_heights
       get_perturbation
       isentropic_interpolation
       isentropic_interpolation_as_dataset
       nearest_intersection_idx
       parse_angle
+      peak_persistence
       reduce_point_density
       resample_nn_1d

docs/api/references.rst

@@ -114,6 +114,11 @@ References
 .. [Holton2004] Holton, J. R., 2004: *An Introduction to Dynamic Meteorology*. 4th ed.
            Academic Press, 535 pp.
 
+.. [Huber2021] Huber, S., 2020: Persistent Homology in Data Science.
+           *Proc. Int. Data Sci. Conf.*, doi:`10.1007/978-3-658-32182-6_13
+           <https://doi.org/10.1007/978-3-658-32182-6_13>`_.
+           `[PDF] <https://www.sthu.org/research/publications/files/Hub20.pdf>`_
+
 .. [IAPWS1995] The International Association for the Properties of Water and Steam, 1995:
            `Revised Release on the IAPWS Formulation 1995 for the Thermodynamic Properties
            of Ordinary Water Substance for General and Scientific Use (updated

docs/userguide/gempak.rst

@@ -486,11 +486,11 @@ blue is uncertain of parity, and white is unevaluated.
         <td></td>
       </tr>
       <tr>
-        <td>HIGH(S, RADIUS)</td>
-        <td>Relative maxima over a grid</td>
-        <td></td>
-        <td></td>
-        <td></td>
+        <td class="tg-implemented">HIGH(S, RADIUS)</td>
+        <td class="tg-implemented">Relative maxima over a grid</td>
+        <td class="tg-implemented"><a href="../api/generated/metpy.calc.find_peaks.html#metpy.calc.find_peaks">metpy.calc.find_peaks</a></td>
+        <td class="tg-yes">Yes</td>
+        <td class="tg-no">No</td>
         <td></td>
       </tr>
       <tr>
@@ -534,11 +534,11 @@ blue is uncertain of parity, and white is unevaluated.
         <td></td>
       </tr>
       <tr>
-        <td>LOWS(S, RADIUS)</td>
-        <td>Relative minima over a grid</td>
-        <td></td>
-        <td></td>
-        <td></td>
+        <td class="tg-implemented">LOWS(S, RADIUS)</td>
+        <td class="tg-implemented">Relative minima over a grid</td>
+        <td class="tg-implemented"><a href="../api/generated/metpy.calc.find_peaks.html#metpy.calc.find_peaks">metpy.calc.find_peaks</a></td>
+        <td class="tg-yes">Yes</td>
+        <td class="tg-no">No</td>
         <td></td>
       </tr>
       <tr>

examples/calculations/High_Low_Analysis.py

@@ -0,0 +1,62 @@
+# Copyright (c) 2025 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+=================
+High/Low Analysis
+=================
+
+This uses MetPy's `find_peaks` function to automatically identify locations of high and low
+centers, and then plots them on a map.
+"""
+
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import matplotlib.pyplot as plt
+import xarray as xr
+
+from metpy.calc import find_peaks
+from metpy.cbook import get_test_data
+from metpy.plots import add_metpy_logo, scattertext
+from metpy.units import units
+
+###########################################
+# Start by loading some data from our sample GFS model dataset. Pull out the geopotential
+# heights field for the 850 hPa layer, as well as grab the projection metadata.
+data = xr.open_dataset(get_test_data('GFS_test.nc', as_file_obj=False)).metpy.parse_cf()
+mslp = data.Geopotential_height_isobaric.metpy.sel(vertical=850 * units.hPa).squeeze()
+dataproj = mslp.metpy.cartopy_crs
+
+
+###########################################
+# Here we use `find_peaks` to find the locations of the highs and then the lows
+h_y, h_x = find_peaks(mslp.values)
+l_y, l_x = find_peaks(mslp.values, maxima=False)
+
+###########################################
+# Plot the analyzed locations on top of the contours of height on a map
+fig = plt.figure(figsize=(11., 8.))
+ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-95))
+ax.contour(mslp.metpy.x, mslp.metpy.y, mslp, range(0, 2000, 30),
+           colors='k', linewidths=1.25, linestyles='solid', transform=dataproj)
+
+# Using scattertext() plot the high centers using a red 'H' and put the height value
+# below the 'H' using a smaller font.
+scattertext(ax, mslp.metpy.x[h_x], mslp.metpy.y[h_y], 'H', size=20, color='red',
+            fontweight='bold', transform=dataproj)
+scattertext(ax, mslp.metpy.x[h_x], mslp.metpy.y[h_y], mslp.values[h_y, h_x], formatter='.0f',
+            size=12, color='red', loc=(0, -15), fontweight='bold', transform=dataproj)
+
+# Now do the same for the lows using a blue 'L'
+scattertext(ax, mslp.metpy.x[l_x], mslp.metpy.y[l_y], 'L', size=20, color='blue',
+            fontweight='bold', transform=dataproj)
+scattertext(ax, mslp.metpy.x[l_x], mslp.metpy.y[l_y], mslp.values[l_y, l_x], formatter='.0f',
+            size=12, color='blue', loc=(0, -15), fontweight='bold', transform=dataproj)
+
+ax.add_feature(cfeature.OCEAN)
+ax.add_feature(cfeature.LAND)
+ax.add_feature(cfeature.COASTLINE)
+
+ax.set_title('Automated 850hPa High and Low Locations')
+add_metpy_logo(fig, 275, 295, size='large')
+plt.show()

examples/plots/Plotting_Surface_Analysis.py

@@ -17,8 +17,8 @@
 
 from metpy.cbook import get_test_data
 from metpy.io import parse_wpc_surface_bulletin
-from metpy.plots import (add_metpy_logo, ColdFront, OccludedFront, StationaryFront,
-                         StationPlot, WarmFront)
+from metpy.plots import (add_metpy_logo, ColdFront, OccludedFront, scattertext,
+                         StationaryFront, WarmFront)
 
 ###########################################
 # Define a function that can be used to readily plot a bulletin that has been parsed into a
@@ -45,9 +45,10 @@ def plot_bulletin(ax, data):
     for field in ('HIGH', 'LOW'):
         rows = data[data.feature == field]
         x, y = zip(*((pt.x, pt.y) for pt in rows.geometry), strict=False)
-        sp = StationPlot(ax, x, y, transform=ccrs.PlateCarree(), clip_on=True)
-        sp.plot_text('C', [field[0]] * len(x), **complete_style[field])
-        sp.plot_parameter('S', rows.strength, **complete_style[field])
+        scattertext(ax, x, y, field[0],
+                    **complete_style[field], transform=ccrs.PlateCarree(), clip_on=True)
+        scattertext(ax, x, y, rows.strength, formatter='.0f', loc=(0, -10),
+                    **complete_style[field], transform=ccrs.PlateCarree(), clip_on=True)
 
     # Handle all the boundary types
     for field in ('WARM', 'COLD', 'STNRY', 'OCFNT', 'TROF'):

src/metpy/calc/tools.py

@@ -5,6 +5,7 @@
 import contextlib
 import functools
 from inspect import Parameter, signature
+import itertools
 from operator import itemgetter
 import textwrap
 
@@ -1946,3 +1947,114 @@ def _remove_nans(*variables):
     for v in variables:
         ret.append(v[~mask])
     return ret
+
+
+def _neighbor_inds(y, x):
+    """Generate index (row, col) pairs for each neighbor of (x, y)."""
+    incs = (-1, 0, 1)
+    for dx, dy in itertools.product(incs, incs):
+        yield y + dy, x + dx
+
+
+def _find_uf(uf, item):
+    """Find the root item for ``item`` in the union find structure ``uf``."""
+    # uf is a dictionary mapping item->parent. Loop until we find parent=parent.
+    while (next_item := uf[item]) != item:
+        uf[item] = uf[next_item]
+        item = next_item
+    return item
+
+
+@exporter.export
+def peak_persistence(arr, *, maxima=True):
+    """Return all local extrema points ordered by their persistence.
+
+    This uses the concept of persistent homology to find peaks ordered by their persistence.
+    [Huber2021]_ This essentially works akin to a falling water level and seeing how long a
+    peak persists until the falling level allows connection to a larger, neighboring peak.
+    NaN points are ignored.
+
+    Parameters
+    ----------
+    arr : array-like
+        2-dimensional array of numeric values to search
+    maxima : bool, optional
+        Whether to find local maxima, defaults to True. If False, local minima will be found
+        instead.
+
+    Returns
+    -------
+    list[((int, int), float)]
+        Point indices and persistence values in descending order of persistence
+
+    See Also
+    --------
+    find_peaks
+
+    """
+    # Get the collection of points and values in descending strength of peak.
+    points = sorted((item for item in np.ndenumerate(arr) if not np.isnan(item[1])),
+                    key=lambda i: i[1], reverse=maxima)
+
+    # The global max will never be merged and thus should be added to the final list
+    # of persisted points
+    per = {points[0][0]: np.inf}
+
+    # Loop over all points and add them to the set (a dict storing as a union-find data
+    # structure) one-by-one
+    peaks = {}
+    for pt, val in points:
+        # Look to see if any neighbors of this point are attached to any existing peaks
+        if already_done := {_find_uf(peaks, n) for n in _neighbor_inds(*pt) if n in peaks}:
+
+            # Get these peaks in order of value
+            biggest, *others = sorted(already_done, key=lambda i: arr[i], reverse=maxima)
+
+            # Connect our point to the biggest peak
+            peaks[pt] = biggest
+
+            # Any other existing peaks will join to this biggest and will end their
+            # persistence, denoted as the difference between their original level and the
+            # current level
+            for neighbor in others:
+                peaks[neighbor] = biggest
+                if arr[neighbor] != val:
+                    per[neighbor] = abs(arr[neighbor] - val)
+        else:
+            peaks[pt] = pt
+
+    # Return the points in descending order of persistence
+    return sorted(per.items(), key=lambda i: i[1], reverse=True)
+
+
+@exporter.export
+def find_peaks(arr, *, maxima=True, iqr_ratio=2):
+    """Search array for significant peaks (or valleys).
+
+    Parameters
+    ----------
+    arr: array-like
+        2-dimensional array of numeric values to search
+    maxima: bool, optional
+        Whether to find local maxima, defaults to True. If False, local minima will be found
+        instead.
+    iqr_ratio: float, optional
+        Ratio of interquantile range (Q1 - Q3) of peak persistence values to use as a
+        threshold (when added to Q3) for significance of peaks. Defaults to 2.
+
+    Returns
+    -------
+    row indices, column indices
+        Locations of significant peaks
+
+    See Also
+    --------
+    peak_persistence
+
+    """
+    peaks = peak_persistence(arr, maxima=maxima)
+    q1, q3 = np.percentile([p[-1] for p in peaks], (25, 75))
+    thresh = q3 + iqr_ratio * (q3 - q1)
+    return map(list,
+               zip(*(it[0] for it in itertools.takewhile(lambda i: i[1] > thresh, peaks)),
+                   strict=True))

src/metpy/plots/__init__.py

@@ -3,8 +3,6 @@
 # SPDX-License-Identifier: BSD-3-Clause
 r"""Contains functionality for making meteorological plots."""
 
-# Trigger matplotlib wrappers
-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)
@@ -13,6 +11,7 @@
 from .patheffects import *  # noqa: F403
 from .skewt import *  # noqa: F403
 from .station_plot import *  # noqa: F403
+from .text import *  # noqa: F403
 from .wx_symbols import *  # noqa: F403
 from ..package_tools import set_module
 
@@ -21,6 +20,7 @@
 __all__.extend(patheffects.__all__)  # pylint: disable=undefined-variable
 __all__.extend(skewt.__all__)  # pylint: disable=undefined-variable
 __all__.extend(station_plot.__all__)  # pylint: disable=undefined-variable
+__all__.extend(text.__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'])