从xarray走向netCDF处理(四)：合并与计算
前面有关xarray已经讲了3期了，介绍了数据索引，数据结构还有插值和掩膜。今天这是最后一期介绍用xarray处理nc数据了，打算聊一下如何做数据合并与计算。
数据合并
数据合并主要是两种形式
- 维度的拼接：如将日数据合成为年数据，就属于在时间维度上的合并。
- 变量的合并：如将多个物理量合到同一个Dataset中。
xarray围绕着这两种合并方式介绍了concatenate, merge, combine, update四种方法。
我在这里就挑最常用的跟大家聊聊。
维度拼接
使用 concat() 方法可以实现维度的拼接。
下面是演示数据，来源于2018年和2019年前三个月的ERA-Interim月平均数据。
1>>> ds2018
2<xarray.Dataset>
3Dimensions: (latitude: 241, longitude: 480, time: 12)
4Coordinates:
5 * longitude (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
6 * latitude (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
7 * time (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
8Data variables:
9 u10 (time, latitude, longitude) float32 ...
10 v10 (time, latitude, longitude) float32 ...
11 t2m (time, latitude, longitude) float32 ...
12Attributes:
13 Conventions: CF-1.6
14
15>>> ds2019
16<xarray.Dataset>
17Dimensions: (latitude: 241, longitude: 480, time: 3)
18Coordinates:
19 * longitude (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
20 * latitude (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
21 * time (time) datetime64[ns] 2019-01-01 2019-02-01 2019-03-01
22Data variables:
23 u10 (time, latitude, longitude) float32 ...
24 v10 (time, latitude, longitude) float32 ...
25 t2m (time, latitude, longitude) float32 ...
26Attributes:
27 Conventions: CF-1.6ds2018时间维度为12，ds2019时间维度为3，下面使用 concat() 合并后时间维度为15
1>>> xr.concat([ds2018, ds2019], dim='time')
2<xarray.Dataset>
3Dimensions: (latitude: 241, longitude: 480, time: 15)
4Coordinates:
5 * longitude (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
6 * latitude (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
7 * time (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2019-03-01
8Data variables:
9 u10 (time, latitude, longitude) float32 -0.9599868 ... 4.5229325
10 v10 (time, latitude, longitude) float32 3.1737509 ... -2.289166
11 t2m (time, latitude, longitude) float32 248.46857 ... 225.19632
12Attributes:
13 Conventions: CF-1.6变量合并
使用 merge() 方法，可以将ds2018中的u10和ds2019中的t2m合并到一起，而且在时间维上缺失会自动设置为nan。
1>>> xr.merge([ds2018.u10, ds2019.t2m])
2<xarray.Dataset>
3Dimensions: (latitude: 241, longitude: 480, time: 15)
4Coordinates:
5 * time (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2019-03-01
6 * longitude (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
7 * latitude (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
8Data variables:
9 u10 (time, latitude, longitude) float32 -0.9599868 -0.9599868 ... nan
10 t2m (time, latitude, longitude) float32 nan nan ... 225.19632数据计算
最基本的计算就是进行加减乘除，任意一个DataArray或者Dataset都可以直接进行四则运算。
除此以外，xarray还可以帮你快速地求出平均值，方差，最小值，最大值等。你可以指定具体对那个维度进行计算，如果不指定维度默认会对所有维度进行计算。
比如要对经、纬两个维度进行平均，最后的结果只有时间维的12个值。
1>>> ds2018.mean(dim=['latitude', 'longitude'])
2<xarray.Dataset>
3Dimensions: (time: 12)
4Coordinates:
5 * time (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
6Data variables:
7 u10 (time) float32 -0.120867714 0.13738841 ... -0.016953295 0.08418254
8 v10 (time) float32 -0.21417202 -0.12471106 ... -0.082666814
9 t2m (time) float32 277.57446 277.32916 276.72095 ... 278.1613 278.01758如果对时间维进行求方差，则结果会保留空间场上的信息。
1>>> ds2018.std(dim='time')
2<xarray.Dataset>
3Dimensions: (latitude: 241, longitude: 480)
4Coordinates:
5 * longitude (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
6 * latitude (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
7Data variables:
8 u10 (latitude, longitude) float32 1.9192954 1.9192954 ... 1.2133
9 v10 (latitude, longitude) float32 1.3066719 1.3066719 ... 1.577495
10 t2m (latitude, longitude) float32 9.5681305 9.5681305 ... 11.313364而且xarray在时间维上的计算还有很多贴心的用法，比如月数据转年数据，月数据转季节数据。
1>>> ds2018.groupby('time.season').min(dim='time')
2<xarray.Dataset>
3Dimensions: (latitude: 241, longitude: 480, season: 4)
4Coordinates:
5 * longitude (longitude) float32 0.0 0.75 1.5 2.25 ... 357.75 358.5 359.25
6 * latitude (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
7 * season (season) object 'DJF' 'JJA' 'MAM' 'SON'
8Data variables:
9 u10 (season, latitude, longitude) float32 -3.2505732 ... 3.4194348
10 v10 (season, latitude, longitude) float32 0.9903783 ... -3.318876
11 t2m (season, latitude, longitude) float32 248.46857 ... 214.61685一个栗子
1import numpy as np
2import xarray as xr
3import cartopy.crs as ccrs
4import cartopy.feature as cfeat
5from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
6import matplotlib.pyplot as plt
7
8def create_map(ax):
9 # 设置地图属性
10 ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8) # 加载分辨率为50的国界线
11 ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6) # 加载分辨率为50的海岸线
12 ax.add_feature(cfeat.RIVERS.with_scale('50m')) # 加载分辨率为50的河流
13 ax.add_feature(cfeat.LAKES.with_scale('50m')) # 加载分辨率为50的湖泊
14 # 设置网格点属性
15 gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
16 linewidth=1.2, color='k', alpha=0.5, linestyle='--')
17 gl.xlabels_top = False # 关闭顶端的经纬度标签
18 gl.ylabels_right = False # 关闭右侧的经纬度标签
19 gl.xformatter = LONGITUDE_FORMATTER # x轴设为经度的格式
20 gl.yformatter = LATITUDE_FORMATTER # y轴设为纬度的格式
21 return ax
22
23if __name__ == '__main__':
24 # 创建画图空间
25 proj = ccrs.PlateCarree()
26 fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(7, 14), subplot_kw={'projection': proj})
27 ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
28 temp = (ds['t2m'] - 273.15).groupby('time.season').mean('time')
29 # --画图
30 cbar_kwargs = {
31 'label': '2m temperature (℃)',
32 'ticks': np.arange(-30, 30+5, 5)
33 }
34 levels = np.arange(-30, 30+1, 1)
35 for i, season in enumerate(('DJF', 'MAM', 'JJA', 'SON')):
36 ax = create_map(axes[i])
37 temp.sel(season=season).plot.contourf(ax=ax, levels=levels, cmap='Spectral_r', extend='both',
38 cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
39 fig.show()举报
腾讯云开发者
