前面有关xarray已经讲了3期了,介绍了数据索引数据结构还有插值和掩膜。今天这是最后一期介绍用xarray处理nc数据了,打算聊一下如何做数据合并与计算

数据合并

数据合并主要是两种形式

  1. 维度的拼接:如将日数据合成为年数据,就属于在时间维度上的合并。

  2. 变量的合并:如将多个物理量合到同一个*Dataset*中。

xarray围绕着这两种合并方式介绍了concatenate, merge, combine, update四种方法。我在这里就挑最常用的跟大家聊聊。

维度拼接

使用 concat() 方法可以实现维度的拼接。

下面是演示数据,来源于2018年和2019年前三个月的ERA-Interim月平均数据。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
>>> ds2018
<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

>>> ds2019
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, time: 3)
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] 2019-01-01 2019-02-01 2019-03-01
Data variables:
    u10        (time, latitude, longitude) float32 ...
    v10        (time, latitude, longitude) float32 ...
    t2m        (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6

ds2018时间维度为12,ds2019时间维度为3,下面使用 concat() 合并后时间维度为15

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
>>> xr.concat([ds2018, ds2019], dim='time')
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, time: 15)
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 ... 2019-03-01
Data variables:
    u10        (time, latitude, longitude) float32 -0.9599868 ... 4.5229325
    v10        (time, latitude, longitude) float32 3.1737509 ... -2.289166
    t2m        (time, latitude, longitude) float32 248.46857 ... 225.19632
Attributes:
    Conventions:  CF-1.6

变量合并

使用 merge() 方法,可以将ds2018中的u10ds2019中的t2m合并到一起,而且在时间维上缺失会自动设置为nan

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
>>> xr.merge([ds2018.u10, ds2019.t2m])
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, time: 15)
Coordinates:
  * time       (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2019-03-01
  * 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
Data variables:
    u10        (time, latitude, longitude) float32 -0.9599868 -0.9599868 ... nan
    t2m        (time, latitude, longitude) float32 nan nan ... 225.19632

数据计算

最基本的计算就是进行加减乘除,任意一个DataArray或者Dataset都可以直接进行四则运算。

除此以外,xarray还可以帮你快速地求出平均值,方差,最小值,最大值等。你可以指定具体对那个维度进行计算,如果不指定维度默认会对所有维度进行计算。

比如要对经、纬两个维度进行平均,最后的结果只有时间维的12个值。

1
2
3
4
5
6
7
8
9
>>> ds2018.mean(dim=['latitude', 'longitude'])
<xarray.Dataset>
Dimensions:  (time: 12)
Coordinates:
  * time     (time) datetime64[ns] 2018-01-01 2018-02-01 ... 2018-12-01
Data variables:
    u10      (time) float32 -0.120867714 0.13738841 ... -0.016953295 0.08418254
    v10      (time) float32 -0.21417202 -0.12471106 ... -0.082666814
    t2m      (time) float32 277.57446 277.32916 276.72095 ... 278.1613 278.01758

如果对时间维进行求方差,则结果会保留空间场上的信息。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
>>> ds2018.std(dim='time')
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480)
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
Data variables:
    u10        (latitude, longitude) float32 1.9192954 1.9192954 ... 1.2133
    v10        (latitude, longitude) float32 1.3066719 1.3066719 ... 1.577495
    t2m        (latitude, longitude) float32 9.5681305 9.5681305 ... 11.313364

而且xarray在时间维上的计算还有很多贴心的用法,比如月数据转年数据,月数据转季节数据。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
>>> ds2018.groupby('time.season').min(dim='time')
<xarray.Dataset>
Dimensions:    (latitude: 241, longitude: 480, season: 4)
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
  * season     (season) object 'DJF' 'JJA' 'MAM' 'SON'
Data variables:
    u10        (season, latitude, longitude) float32 -3.2505732 ... 3.4194348
    v10        (season, latitude, longitude) float32 0.9903783 ... -3.318876
    t2m        (season, latitude, longitude) float32 248.46857 ... 214.61685

一个栗子

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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

def create_map(ax):
    # 设置地图属性
    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 ax
    
if __name__ == '__main__':
    # 创建画图空间
    proj = ccrs.PlateCarree()
    fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(7, 14), subplot_kw={'projection': proj})
    ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
    temp = (ds['t2m'] - 273.15).groupby('time.season').mean('time')
    # --画图
    cbar_kwargs = {
       'label': '2m temperature (℃)',
       'ticks': np.arange(-30, 30+5, 5)
    }
    levels = np.arange(-30, 30+1, 1)
    for i, season in enumerate(('DJF', 'MAM', 'JJA', 'SON')):
        ax = create_map(axes[i])
        temp.sel(season=season).plot.contourf(ax=ax, levels=levels, cmap='Spectral_r', extend='both',
            cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
    fig.show()