如何在Cartopy中绘制高斯网格数据？解决xarray重塑后NaN问题

你的问题出在错误地对已有网格维度的数据做了不必要的重塑操作，ARCO ERA5的单级别预报数据集本身已经自带latitude和longitude作为维度，不需要手动用MultiIndex和unstack来转换，这反而会破坏原有的网格结构，生成大量NaN。

修正后的代码 #

import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

# 读取数据
sl_forecasts = xr.open_zarr(
    'gs://gcp-public-data-arco-era5/co/single-level-forecast.zarr/', 
    chunks={'time': 48},
    consolidated=True,
)

# 筛选目标时间、步长和变量，直接保留原维度
tp_data = sl_forecasts['tp'].sel(
    time='2021-02-27T06:00:00.000000000',
    step='3600000000000'
)

# 绘图
proj = ccrs.PlateCarree()
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=proj))

# 直接用xarray的plot方法，自动处理维度和投影
tp_data.plot.contourf(
    ax=ax,
    transform=proj,
    cmap='Blues',
    add_colorbar=True
)
ax.coastlines()
plt.show()

为什么原来的代码会产生大量NaN？ #

你创建的MultiIndex是将所有经度和纬度值一一配对，但原数据的latitude和longitude是两个独立的维度（形成二维网格）。当你执行unstack('values')时，会生成一个包含所有经纬度笛卡尔积的新网格，而原数据只有对应网格点的有效值，其余非对应位置自然被填充为NaN，导致绘图失败。

额外提示 #

xarray和Cartopy有很好的集成，直接使用xarray的plot方法可以省去手动提取values、lon、lat的步骤，同时自动处理坐标匹配，避免维度混乱。

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