如何用Xarray提取全球数据的美国经纬度子集并绘图？

我用Python处理Xarray的global_data数据（维度：lat=73、lon=144；纬度范围90.0至-90.0，经度范围0.0至357.5），想要提取美国区域子集（纬度25_{50，经度-125}-66）并绘图。尝试过.sel、.isel和slice方法，但要么无数据可绘，要么绘图位置错误（比如显示非洲区域）。目前靠试错写出一段.isel代码能生成美国区域图，但不确定正确性和原理，求解答。

相关数据结构：

<xarray.DataArray (lat: 73, lon: 144)>
array([[ 1.1423118 ,  1.1418152 ,  1.1437986 , ...,  1.1403227 ,
         1.1362497 ,  1.1439507 ],
       [ 0.49379024,  0.42622158,  0.3565474 , ...,  0.70473236,
         0.63061965,  0.5644286 ],
       [ 0.1380711 ,  0.19678137,  0.2836361 , ...,  0.24298143,
         0.16488086,  0.12564908],
       ...,
       [ 0.18887411,  0.25456694,  0.30384657, ..., -0.08306076,
         0.01299069,  0.10468105],
       [ 0.37176454,  0.4389612 ,  0.50888765, ...,  0.16366327,
         0.23381352,  0.30179456],
       [ 0.6100794 ,  0.61286193,  0.6167843 , ...,  0.6154521 ,
         0.6117071 ,  0.610261  ]], dtype=float32)
Coordinates:
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5

试错代码：

usa_data = global_data.isel(lon=slice(-55,-23),lat=slice(13,29))
ax = plt.axes(projection=ccrs.PlateCarree())
contourf = anomspeed_us.plot.contourf(
    ax=ax, levels=np.arange(-5,5.1,0.5),extend=
    'both',cmap='RdBu_r',
    add_colorbar=True,
    cbar_kwargs={'label': 'Anomalous Wind Speed'})
ax.coastlines()
states = cfeature.NaturalEarthFeature(category='cultural', 
                                      name='admin_1_states_provinces_lines',
                                      scale='50m', facecolor='none')
ax.add_feature(states, linewidth=0.5, edgecolor='black')
ax.add_feature(cfeature.BORDERS, linestyle='-')
ax.gridlines()
plt.title('Anomalous Wind Speed for Week')
plt.show()

问题根源 #

你的数据经度是**0°_{360°的东经格式**，但美国的经度是西经（-125°}-66°），直接用sel(lon=slice(-125,-66))会匹配不到数据——因为数据里没有负的经度值，自然会返回空或者错误区域。

试错代码的有效性原理 #

拆解你用的isel参数：

• 纬度slice(13,29)：数据纬度从90°开始递减（步长2.5°），索引n对应的纬度为90 - n*2.5。索引13对应57.5°，索引29对应17.5°，这个范围刚好覆盖美国本土的25°~50°，所以能显示目标区域。
• 经度slice(-55,-23)：经度共144个点，负索引对应从末尾往前数：索引-55等价于144-55=89，对应经度89*2.5=222.5°；索引-23等价于144-23=121，对应经度121*2.5=302.5°。而美国西经范围转成东经是360-125=235°到360-66=294°，你的索引范围刚好覆盖了这个区间，所以能提取到正确数据。

但这种依赖索引的方法完全绑定当前数据的维度和分辨率，一旦数据参数变化就会失效，非常不稳健。

正确的稳健提取方法 #

方法1：手动转换经度范围后用.sel #

把西经转成东经（西经值 + 360），同时注意纬度是递减的，slice要从高纬度到低纬度：

# 转换美国西经范围为东经：-125°→235°，-66°→294°
usa_data = global_data.sel(
    lat=slice(50, 25),  # 纬度从90到-90递减，所以slice(高纬度, 低纬度)
    lon=slice(235, 294)
)

# 修正绘图变量名并绘图
ax = plt.axes(projection=ccrs.PlateCarree())
contourf = usa_data.plot.contourf(
    ax=ax, levels=np.arange(-5,5.1,0.5),
    extend='both', cmap='RdBu_r',
    add_colorbar=True,
    cbar_kwargs={'label': 'Anomalous Wind Speed'}
)
ax.coastlines()
states = cfeature.NaturalEarthFeature(
    category='cultural', 
    name='admin_1_states_provinces_lines',
    scale='50m', facecolor='none'
)
ax.add_feature(states, linewidth=0.5, edgecolor='black')
ax.add_feature(cfeature.BORDERS, linestyle='-')
ax.gridlines()
plt.title('Anomalous Wind Speed for Week')
plt.show()

方法2：转换经度坐标为-180°~180°格式 #

如果经常需要处理跨0°经度的区域，可以先把数据经度转换为-180°~180°的西经/东经格式，之后就能直接用直观的西经范围筛选：

# 将经度从0-360转换为-180-180
global_data_180 = global_data.assign_coords(lon=(((global_data.lon + 180) % 360) - 180))
# 直接用西经范围筛选美国区域
usa_data = global_data_180.sel(lat=slice(50,25), lon=slice(-125,-66))

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