NCEP NAM-12km GRIB2数据使用MetPy无法为xarray.DataArray分配x/y坐标

核心原因 #

你遇到的assign_y_x报错是因为手动填写的Lambert投影参数和NCEP NAM 218网格实际使用的投影参数存在微小偏差，导致经纬度转投影坐标后无法按公差要求降为1D数组。

解决方案1：使用GRIB自带的投影参数（推荐） #

不要手动硬编码投影参数，直接从读取的GRIB文件属性中提取官方参数，避免偏差：

# 读取数据步骤和你原有代码一致
data = xr.open_dataset(filename, engine="cfgrib",filter_by_keys={'typeOfLevel': 'isobaricInhPa','name':'Temperature'})

# 从GRIB属性提取正确的投影参数赋值
data1 = data.metpy.assign_crs(
    grid_mapping_name='lambert_conformal_conic',
    latitude_of_projection_origin=data.attrs['GRIB_LaD'],
    longitude_of_central_meridian=data.attrs['GRIB_LoV'],
    standard_parallel=(data.attrs['GRIB_Latin1'], data.attrs['GRIB_Latin2']),
    earth_radius=6371229.0
)

# 调大公差阈值（单位：米），避免浮点误差导致的分配失败
data2 = data1.metpy.assign_y_x(tolerance=1e3)

解决方案2：手动赋值固定网格坐标 #

NAM 12km CONUS的218号网格的x/y坐标是固定的，直接生成赋值即可，无需投影转换：

# 生成218网格的标准x/y坐标（单位：米）
x_coords = np.linspace(-4223618.75, 4136381.25, 614)
y_coords = np.linspace(-831630.0, 4296370.0, 428)

# 赋值给数据集并添加必要的units属性
data2 = data1.assign_coords(
    x = ('x', x_coords, {'units': 'meters'}),
    y = ('y', y_coords, {'units': 'meters'})
)

完成上述操作后，你就可以正常调用metpy.interpolate.cross_section生成横截面了。

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