Adding annotations for High and Low pressure symbols in Mean Sea Level Pressure field

CGHA 0
Viewed 47

I'm currently plotting the Mean Sea Level Pressure (MSLP) data from ERA5 Reanalysis (sample netcdf data is included) using the 'recipe' from Unidata.

The reproducible code (I tried to be as minimal as I can) runs like this

```python
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
from scipy.ndimage import gaussian_filter
import xarray as xr

#This function is from Unidata recipe as indicated above used for plotting max/min points 
def plot_maxmin_points(lon, lat, data, extrema, nsize, symbol, color='k',
                       plotValue=True, transform=None):
    
    from scipy.ndimage import maximum_filter, minimum_filter

    if (extrema == 'max'):
        data_ext = maximum_filter(data, nsize, mode='nearest')
    elif (extrema == 'min'):
        data_ext = minimum_filter(data, nsize, mode='nearest')
    else:
        raise ValueError('Value for hilo must be either max or min')

    mxy, mxx = np.where(data_ext == data)

    for i in range(len(mxy)):
        ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]], symbol, color=color, size=24,
                clip_on=True, horizontalalignment='center', verticalalignment='center',
                transform=transform)
        ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]],
                '\n' + str(np.int(data[mxy[i], mxx[i]])),
                color=color, size=12, clip_on=True, fontweight='bold',
                horizontalalignment='center', verticalalignment='top', transform=transform)

#Open the included ERA5 data in this post
ds = xr.open_dataset('EU_Single Levels.nc')

lats = ds.latitude
lons = ds.longitude

# Select and grab data
mslp = ds['msl']

# Apply gaussian filter to the MSLP
mslp_surf = gaussian_filter(mslp.data[30], sigma=3.0) * units.pascal
mslp_hpa = mslp_surf.to(units.hectopascal)

mapcrs = ccrs.PlateCarree()

datacrs = ccrs.PlateCarree()

# Create plot axes with proper projection
fig = plt.figure(1, figsize=(14, 12))
ax = plt.subplot(111, projection=datacrs)
ax.set_extent([105, 140, 0, 25], mapcrs)

# Plot MSLP Contour
MSLP_range = np.arange(995, 1010, 2)
prs = ax.contour(lons, lats, mslp_hpa, MSLP_range, colors='k',
                 transform=datacrs)
ax.clabel(prs, fmt='%d')

#Should annotate highs and low pressure areas...not working?
plot_maxmin_points(lons, lats, mslp_hpa, 'max', 50, symbol='H', color='b', transform=datacrs)
plot_maxmin_points(lons, lats, mslp_hpa, 'min', 25, symbol='L', color='r', transform=datacrs)

ax.add_feature(cfeature.COASTLINE.with_scale('10m'), linewidth=1.15)

plt.title('ERA5 Reanalysis: MSLP (hPa)', loc='left')

plt.subplots_adjust(bottom=0, top=1)
plt.show()
```

Which produces a plot like this:

enter image description here

Unfortunately, I ran into some problems when the MSLP plot never 'displayed' the annotations of the MSLP max/min points from the function as created by Unidata. In fact, it also shows a:

DeprecationWarning: elementwise comparison failed; this will raise an error in the future. mxy, mxx = np.where(data_ext == data)

Is there a way to possibly sort out this small issue? Thanks in advance.

metpy matplotlib python-3.x
1 Answers

It does not plot anything because the line mxy, mxx = np.where(data_ext == data) return empty arrays.

You are making an elementwise comparison between arrays whose elements are respectively float32 (in data_ext, returned by maximum_filter) and Quantity objects created by the pint package, dealing with units (in data). Eventhough numbers can be the same, the elementwise comparison will fail (indicated by the warning), and therefore return False everywhere.

To fix this, you can change mxy, mxx = np.where(data_ext == data) to mxy, mxx = np.where(data_ext == np.array(data)).

Two sides note:

  • Having imports inside a function is a bad habit, as it will run the imoprt each time you run the function. Therefore, move from scipy.ndimage import maximum_filter, minimum_filter with the other imports.

  • The function you use to catch the extrema actually return the max/min of box centered on a pixel, and whose size is fixed by your parameter nsize. Therefore, nsize should reflect the 'localness' of what you are looking for, here your pressure structures. As your image size is (101, 141), using a box size of 50 will be too big to capture these structures.

Liris0
Related