海平面气压场高低压标注失效及警告问题求助

问题原因 #

• 弃用警告根源：mslp_hpa是带metpy单位的数组，而scipy.ndimage的滤波函数返回纯numpy数组，两者直接用==做元素比较会触发类型不匹配警告，且无法正确匹配极值点。
• 高低压标注不显示：除了上述类型冲突，调用plot_maxmin_points时传入的nsize参数（50、25）过大，超出数据网格实际范围，导致极值检测失效。

修复方案 #

1. 核心修改点 #

• 传入极值检测函数时，提取带单位数组的纯数值部分（用.magnitude属性），确保和滤波后的numpy数组类型一致。
• 调整nsize参数为合理值（根据ERA5数据分辨率，推荐5-10之间），让极值检测窗口适配数据网格。

2. 修改后的完整代码 #

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

# 极值点标注函数
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('extrema参数只能是max或min')

    # 用np.isclose处理浮点精度问题，避免直接==比较的误差
    mxy, mxx = np.where(np.isclose(data_ext, data, atol=1e-3))

    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)
        if plotValue:
            ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]],
                    '\n' + str(np.int(np.round(data[mxy[i], mxx[i]]))),
                    color=color, size=12, clip_on=True, fontweight='bold',
                    horizontalalignment='center', verticalalignment='top', transform=transform)

# 读取ERA5数据
ds = xr.open_dataset('EU_Single Levels.nc')

lats = ds.latitude
lons = ds.longitude

# 提取海平面气压数据
mslp = ds['msl']

# 高斯滤波并转换单位
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()

# 创建绘图对象
fig = plt.figure(1, figsize=(14, 12))
ax = plt.subplot(111, projection=datacrs)
ax.set_extent([105, 140, 0, 25], mapcrs)

# 绘制海平面气压等值线
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')

# 调用极值标注函数：传入纯数值数组，调整nsize为8
plot_maxmin_points(lons, lats, mslp_hpa.magnitude, 'max', 8, symbol='H', color='b', transform=datacrs)
plot_maxmin_points(lons, lats, mslp_hpa.magnitude, 'min', 8, 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()

3. 关键修改说明 #

• 用np.isclose替代直接==比较：解决浮点数值精度问题，避免因滤波后微小数值差异导致匹配失败。
• 传入mslp_hpa.magnitude：提取纯数值数组，消除单位导致的类型冲突，解决弃用警告。
• 调整nsize=8：适配ERA5数据的网格分辨率，确保极值检测窗口能正确识别高低压中心。

内容的提问来源于stack exchange，提问作者CGHA