Find peaks

[dopplershift]

Description Of Changes

This adds find_peaks (and persistence) as API's for finding peaks, as well as makes scattertext() a public API for plotting this information. This is an alternative to #3110 that makes use of persistent homology as the mathematical underpinning for peak-finding in 2D. This runs in 30s on my mac for a 1377x2145 HRRR grid, which isn't too bad. It's much better for smaller grids.

TODO:

• Improve example
• Clean up scattertext() a bit to handle formatting and broadcasting values

In the meanwhile I wanted to get this up if anyone had API thoughts. Here's a quick example:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
import metpy.plots as mpplots
from metpy.units import units
import numpy as np
import xarray as xr

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

h_y, h_x = mpcalc.find_peaks(mslp.values)
l_y, l_x = mpcalc.find_peaks(mslp.values, maxima=False)

fig = plt.figure(figsize=(8., 8.))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-95))

clevmslp = np.arange(0., 2000., 30)
ax.contour(mslp.metpy.x, mslp.metpy.y, mslp,
           clevmslp, colors='k', linewidths=1.25, linestyles='solid', transform=dataproj)

mpplots.scattertext(ax, mslp.metpy.x[h_x], mslp.metpy.y[h_y], ['H']*2, size=20, color='red',
                    fontweight='bold', transform=dataproj)
mpplots.scattertext(ax, mslp.metpy.x[l_x], mslp.metpy.y[l_y], ['L']*4, size=20, color='blue',
                    fontweight='bold', transform=dataproj)
mpplots.scattertext(ax, mslp.metpy.x[l_x], mslp.metpy.y[l_y], [format(f, '.0f') for f in mslp.values[l_y, l_x]],
                    size=12, color='blue', fontweight='bold', loc=(0, -15), transform=dataproj)

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)

Checklist

• Closes Clean up High/Low Plotting Function A Bit #1234
• Tests added
• Fully documented

[kgoebber]

Generally, I think I like the user experience of the new function. The only thing I'm wondering is whether it would be useful to include the number of peaks in the output of the function. I only say this because when you go to plot with the now public scattertext, you need to know the number of peaks for when you need to plot a given number of H or L "symbols". Clearly this could be gotten after the fact from a simple len() call, but don't know if it would be more convenient to have that number be a third output from the function.

[dopplershift]

@kgoebber My plan is to tweak scattertext() so that it broadcasts arguments correctly and you don't need to supply that for plotting. Given that, I don't think there's a need to directly return it and any other usage can use len().

[kgoebber]

@dopplershift That was another thought I had too, which I think is an even better solution! And if I could actually read, I would have seen that it was you put on the to do list. D'oh!

[dopplershift]

Ok, now the example (which is included) reads as:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
import metpy.plots as mpplots
from metpy.units import units
import numpy as np
import xarray as xr

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

h_y, h_x = mpcalc.find_peaks(mslp.values)
l_y, l_x = mpcalc.find_peaks(mslp.values, maxima=False)

fig = plt.figure(figsize=(8., 8.))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(central_longitude=-95))

clevmslp = np.arange(0., 2000., 30)
ax.contour(mslp.metpy.x, mslp.metpy.y, mslp,
           clevmslp, colors='k', linewidths=1.25, linestyles='solid', transform=dataproj)

mpplots.scattertext(ax, mslp.metpy.x[h_x], mslp.metpy.y[h_y], 'H', size=20, color='red',
                    fontweight='bold', transform=dataproj)
mpplots.scattertext(ax, mslp.metpy.x[l_x], mslp.metpy.y[l_y], 'L', size=20, color='blue',
                    fontweight='bold', transform=dataproj)
mpplots.scattertext(ax, mslp.metpy.x[l_x], mslp.metpy.y[l_y], mslp.values[l_y, l_x],
                    formatter='.0f',
                    size=12, color='blue', fontweight='bold', loc=(0, -15), transform=dataproj)

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)