MetPy绘制HRRR剖面时如何将经纬度转换为米/千米距离单位

AI 小北

2026-8-27

问题说明

基于归档HRRR Grib2数据制作自定义垂直剖面，已参考MetPy剖面示例解决全部文件格式问题，目前可正常生成x轴为纬度、y轴为等压面气压的初步剖面结果：
目标是将x轴替换为断面线起点到终点的沿线距离（单位米/千米），用于判定湖风、出流等近地面气象特征的水平尺度，目标效果如下：
现有实现代码：

# 经纬度定义断面起止点
startpoint = (42.857, -85.381)
endpoint = (43.907, -83.910)

# 数据读取
grib = pygrib.open('file.grib2')
msg = grib.message(1) 

# 读取grib为xarray数据集，HRRR默认使用兰伯特圆锥投影
ds = xr.open_dataset(file, engine="cfgrib",
                     backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa'}})
ds = ds.metpy.assign_crs(CRS(msg.projparams).to_cf()).metpy.assign_y_x()

# 调用MetPy接口生成断面
cross = cross_section(ds, startpoint, endpoint).set_coords(('latitude', 'longitude'))
 
# 计算诊断变量
temperature = cross['t']
pressure = cross['isobaricInhPa']
cross['Potential_temperature'] = mpcalc.potential_temperature(cross['isobaricInhPa'],cross['t'])
cross['u_wind'] = cross['u'].metpy.convert_units('knots')
cross['v_wind'] = cross['v'].metpy.convert_units('knots')
cross['t_wind'], cross['n_wind'] = mpcalc.cross_section_components(cross['u_wind'],cross['v_wind'])
cross['qv'] = cross['q'] *1000* units('g/kg')

# 初步绘图
fig = plt.figure(1, figsize=(20,9))
ax = plt.axes()
temp = ax.contourf(cross['latitude'], cross['isobaricInhPa'], cross['qv'], 100, cmap='rainbow')
clb = fig.colorbar(temp)
clb.set_label('g $\mathregular{kg^{-1}}$')
theta_contour = ax.contour(cross['latitude'], cross['isobaricInhPa'], cross['Potential_temperature'],
                           400, colors='k', linewidths=2)
theta_contour.clabel(theta_contour.levels[1::2], fontsize=8, colors='k', inline=1,
                     inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True)
ax.set_ylim(775,1000)
ax.invert_yaxis()
plt.title('HRRR contour fill of Mixing ratio(g/kg), contour of Potential Temperature (K),\n Tangential/Normal Winds (knots)')
plt.title('Run: '+date+'', loc='left', fontsize='small')
plt.title('Valid: '+date+' '+f_hr, loc='right', fontsize='small')
plt.xlabel('Latitude')
plt.ylabel('Pressure (hPa)')
wind_slc_vert = list(range(0, 19, 2)) + list(range(19, 29))
wind_slc_horz = slice(5, 100, 5)
ax.barbs(cross['latitude'][wind_slc_horz], cross['isobaricInhPa'][wind_slc_vert],
         cross['t_wind'][wind_slc_vert, wind_slc_horz],
         cross['n_wind'][wind_slc_vert, wind_slc_horz], color='k')
 
# y轴设置为对数坐标
ax.set_yscale('symlog')
ax.set_yticklabels(np.arange(1000, 775,-100))
ax.set_yticks(np.arange(1000, 775,-100))

plt.show()

实现方法

cross_section返回的断面数据已经按起点到终点的顺序排列所有格点，不需要额外做坐标插值，直接计算每个格点到起点的沿程距离替换原有纬度x轴即可，操作分三步：

生成断面数据集后，新增沿程距离变量。优先用投影平面坐标计算，HRRR为兰伯特投影，平面距离精度满足分析需求：

# 提取投影坐标，计算到起点的平面距离，单位转换为千米
proj_x = cross['x'].values
proj_y = cross['y'].values
cross['along_track_km'] = np.sqrt((proj_x - proj_x[0])**2 + (proj_y - proj_y[0])**2) / 1000

如果存在版本兼容问题找不到投影坐标，直接用经纬度计算大圆距离即可：

lats = cross['latitude'].values
lons = cross['longitude'].values
cross['along_track_km'] = mpcalc.great_circle_distance(lats[0], lons[0], lats, lons).metpy.convert_units('km').values

将所有绘图语句中作为x轴输入的cross['latitude']全部替换为cross['along_track_km']，涉及contourf、contour、barbs三个绘图函数的第一个坐标参数。
修改x轴标签，将原plt.xlabel('Latitude')替换为plt.xlabel('沿断面距离 (km)')，按需调整x轴刻度间隔即可。

如果需要同时展示对应位置的经纬度，可在x轴上方添加次坐标轴，映射对应距离点的经纬度刻度，不影响主距离轴的使用。

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