xarray专题再次开讲,错过第一部分的可以先去补个课从xarray走向netCDF处理:数据结构及数据读取。 今天要介绍的就是xarray的索引功能,通过索引你可以对数据进行切片,从整体中提取你所关注的区域、高度或者时间。

索引核心方法

在xarray的官方文档中给出了如下几种索引方式

索引演示

对如下数据进行索引演示:名为dsDataSet,名为tempDataArray

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
>>>ds
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, time: 12)
Coordinates:
  * longitude  (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
  * latitude   (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
  * time       (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Data variables:
    u10        (time, latitude, longitude) float32 ...
    v10        (time, latitude, longitude) float32 ...
    t2m        (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2019-03-28 02:03:39 GMT by grib_to_netcdf-2.12.0: grib_to_n...
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
>>>temp = ds['t2m']
>>>temp
<xarray.DataArray 't2m' (time: 12, latitude: 241, longitude: 480)>
[1388160 values with dtype=float32]
Coordinates:
  * longitude  (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
  * latitude   (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
  * time       (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Attributes:
    units:      K
    long_name:  2 metre temperature

根据位置索引

位置索引是最直接也是最简单的索引方式,但是位置索引只对DataArray有效,对DataSet无效。 下面用两种不同方法获取相同的值。

  • 通过数字索引
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>>>temp[:,1,1]
<xarray.DataArray 't2m' (time: 12)>
array([249.14844, 256.4179 , 247.45125, 254.26143, 267.49988, 273.03592,
       273.91003, 273.30096, 267.9147 , 264.65695, 255.07956, 251.31808],
      dtype=float32)
Coordinates:
    longitude  float32 0.75
    latitude   float32 89.25
  * time       (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Attributes:
    units:      K
    long_name:  2 metre temperature
  • 通过标签索引
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>>>temp.loc[:, 89.25, 0.75]
<xarray.DataArray 't2m' (time: 12)>
array([249.14844, 256.4179 , 247.45125, 254.26143, 267.49988, 273.03592,
       273.91003, 273.30096, 267.9147 , 264.65695, 255.07956, 251.31808],
      dtype=float32)
Coordinates:
    longitude  float32 0.75
    latitude   float32 89.25
  * time       (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Attributes:
    units:      K
    long_name:  2 metre temperature

根据维度名字索引

通过维度的名字就可以不必按照指定的维度顺序对数据进行切片了。对DataArrayDataSet都有效,且方法一致。

  • 通过数字索引
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
>>>ds.isel(longitude=1, time=0)
<xarray.Dataset>
Dimensions:    (latitude: 241)
Coordinates:
    longitude  float32 0.75
  * latitude   (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
    time       datetime64[ns] 2018-01-01
Data variables:
    u10        (latitude) float32 ...
    v10        (latitude) float32 ...
    t2m        (latitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2019-03-28 02:03:39 GMT by grib_to_netcdf-2.12.0: grib_to_n...
  • 通过标签索引
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
>>>ds.sel(longitude=0.75, time='2018-01-01')
<xarray.Dataset>
Dimensions:    (latitude: 241)
Coordinates:
    longitude  float32 0.75
  * latitude   (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
    time       datetime64[ns] 2018-01-01
Data variables:
    u10        (latitude) float32 ...
    v10        (latitude) float32 ...
    t2m        (latitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2019-03-28 02:03:39 GMT by grib_to_netcdf-2.12.0: grib_to_n...

索引及可视化实战

首先引入所需要的库。

1
2
3
4
5
6
7
import arrow
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt

定义一个create_map函数,该函数主要用来创建地图相关的信息,后面画图的时候可以随时调用。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
def create_map():
    # 创建画图空间
    proj = ccrs.PlateCarree()  # 创建坐标系
    fig = plt.figure()  # 创建页面
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  # 创建子图
    # 设置地图属性
    ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8)  # 加载分辨率为50的国界线
    ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6)  # 加载分辨率为50的海岸线
    ax.add_feature(cfeat.RIVERS.with_scale('50m'))  # 加载分辨率为50的河流
    ax.add_feature(cfeat.LAKES.with_scale('50m'))  # 加载分辨率为50的湖泊
    # 设置网格点属性
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                      linewidth=1.2, color='k', alpha=0.5, linestyle='--')
    gl.xlabels_top = False  # 关闭顶端的经纬度标签
    gl.ylabels_right = False  # 关闭右侧的经纬度标签
    gl.xformatter = LONGITUDE_FORMATTER  # x轴设为经度的格式
    gl.yformatter = LATITUDE_FORMATTER  # y轴设为纬度的格式
    return fig, ax

对数据中感兴趣的区域进行提取并简单的可视化。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
# 生成地图
fig, ax = create_map()
# 数据读取及时间平均处理
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
time = ds.time
temp = (ds['t2m'] - 273.15) # 把温度转换为℃
# 区域选择
lon_range = lon[(lon>70) & (lon<140)]
lat_range = lat[(lat>0) & (lat<60)]
temp_region = temp.sel(longitude=lon_range, latitude=lat_range, time='2018-02-01')
# 画图
cbar_kwargs = {
   'label': '2m temperature (℃)',
   'ticks': np.arange(-30, 30+5, 5)
}
levels = np.arange(-30, 30+1, 1)
temp_region.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r', extend='both',
    cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
fig.show()