如何绘制3维xarray.Dataset的2D剖面并获取对应绘图数组

优化实现方案 #

你可以直接利用xarray的广播特性跳过堆叠遍历步骤，几行代码就能得到需要的三个一维数组：

# 1. 计算格点相对于左下角原点的距离（示例为球面近似距离，单位km）
lon_orig = da.lon.min().item()
lat_orig = da.lat.min().item()
# 纬度修正：经度每度距离随纬度升高降低
lat_cos = np.cos(np.deg2rad(lat_orig))
dist = np.sqrt(
    ((da.lon - lon_orig) * lat_cos) ** 2 + 
    (da.lat - lat_orig) ** 2
) * 111  # 经纬度每度近似对应111km

# 2. 直接广播后展平得到三个形状一致的一维数组
z = da.time.broadcast_like(da.air).values.ravel()
dist_arr = dist.broadcast_like(da.air).values.ravel()
air_val = da.air.values.ravel()

*如果你的实际需求是提取沿任意斜线的垂直剖面而非整个矩形区域所有点，可以在初始选点时就把经纬度设为相同维度，直接得到二维剖面数据，后续处理更简单：

# 提取从(220,30)到(280,50)的斜剖面，经纬度共用profile维度
tgt_lon = xr.DataArray(np.linspace(220, 280, num=15), dims="profile")
tgt_lat = xr.DataArray(np.linspace(30, 50, num=15), dims="profile")
da_profile = ds.sel(lon=tgt_lon, lat=tgt_lat, method="nearest")

# 直接计算沿剖面的距离
dist = np.sqrt(
    ((da_profile.lon - lon_orig) * lat_cos) ** 2 + 
    (da_profile.lat - lat_orig) ** 2
) * 111
# 此时da_profile.air维度为(time, profile)，广播后展平即可得到目标数组
z = da_profile.time.broadcast_like(da_profile.air).values.ravel()
dist_arr = dist.broadcast_like(da_profile.air).values.ravel()
air_val = da_profile.air.values.ravel()

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