Xarray一库通关气象编程（持续更新中）

Paper

Xarray库提供了包括数据读写、运算和绘图在内的诸多功能，非常方便灵活，下面会尽可能详细的展示Xarray库的强大美丽，一库通关气象编程。

1 必须知道的基础术语

  1. DataArray （数据数组）： 区别于numpy库中的ndarray，DataArray是具有维度标签信息或其他额外信息的数组，而ndarray仅包含数据。
  2. Dataset (数据集)：可以理解为一个或多个DataArray组成的集合。例如一个.nc文件读进来以后就是一个Dataset。
  3. Variable (变量)：每个DataArray 都包含一个变量。
  4. Dimension (维度)：对于一个气象要素场来说，例如时间，纬度，经度就是这个数据数组的维度。
  5. Coordinate (坐标)：例如在经度维度中，那些经度格点就是它的坐标。

2 DataArray和DataSet

2.1 基础

所以对于一个DataArray来说，它的数据结构应该包含了4个关键属性：

  1. values：数据的值，这里就等同于numpy库中的ndarray了；

  2. dims：每个轴（维度）的名称；

  3. coords：坐标；

  4. attrs：其它属性，例如单位等。

    所以我们可以输入这些属性，通过xr.DataArray函数来创建一个DataArray：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import xarray as xr
import numpy as np
import pandas as pd

temperature = np.random.rand(4, 3)
city = ["Beijing", "Shanghai", "Guangzhou"]
times = pd.date_range("1979-01-01", periods=4)
temperature_da1 = xr.DataArray(temperature, coords=[times, city], dims=["time", "city"])

#<xarray.DataArray (time: 4, city: 3)>
#array([[0.52175765, 0.23541586, 0.82895499],
#       [0.75949885, 0.51938664, 0.99386399],
#       [0.8609056 , 0.93325225, 0.30514099],
#       [0.51485224, 0.59548789, 0.68056358]])
#Coordinates:
#  * time     (time) datetime64[ns] 1979-01-01 1979-01-02 1979-01-03 1979-01-04
#  * city     (city) <U9 'Beijing' 'Shanghai' 'Guangzhou'

如果想偷懒，可以不指定coords和dims，仅输入temperature来快速实现：

1
2
3
4
5
6
7
temperature_da2 = xr.DataArray(temperature)
#<xarray.DataArray (dim_0: 4, dim_1: 3)>
#array([[0.68667675, 0.48846151, 0.82217817],
#       [0.61239246, 0.02478648, 0.12715034],
#       [0.27431044, 0.50271205, 0.253774  ],
#       [0.55916828, 0.72812243, 0.45839827]])
#Dimensions without coordinates: dim_0, dim_1

这时候我们再回顾这些关键属性是是么：

1
2
3
4
5
6
7
8
9
10
11
12
13
print(temperature_da1.values)
#[[0.52175765 0.23541586 0.82895499]
# [0.75949885 0.51938664 0.99386399]
# [0.8609056  0.93325225 0.30514099]
# [0.51485224 0.59548789 0.68056358]]

print(temperature_da1.dims)
#('time', 'city')

print(temperature_da1.coords)
#Coordinates:
#  * time     (time) datetime64[ns] 1979-01-01 1979-01-02 1979-01-03 1979-01-04
#  * city     (city) <U9 'Beijing' 'Shanghai' 'Guangzhou'

当我们有了DataArray后，便可以轻松的创建Dataset从而来生成.nc文件。

1
2
3
4
5
6
7
8
temperature_ds = xr.Dataset({'temperature' : temperature_da1})
#<xarray.Dataset>
#Dimensions:      (time: 4, city: 3)
#Coordinates:
#  * time         (time) datetime64[ns] 1979-01-01 1979-01-02 ... 1979-01-04
#  * city         (city) <U9 'Beijing' 'Shanghai' 'Guangzhou'
#Data variables:
#    temperature  (time, city) float64 0.5218 0.2354 0.829 ... 0.5955 0.6806

可以看到先前的DataArray作为一个variable出现在了Dataset里。

在实际使用中，我们读取到的每个数据文件都是一个DataSet，数据文件中的气象要素则是一个DataArray，同样，我们要生成.nc文件时，也要将数据按DataArray<DataSet<.nc的顺序来操作。这部分内容将放在文件读写中一起介绍。

3 文件读写

3.1 netCDF文件读写

先来写一个.nc文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np
import xarray as xr

# 创建随机的 3D NumPy 数组
shape = (12, 10, 15)
t2m = np.random.rand(*shape)

# 将 NumPy 数组转换为 xarray DataArray
data_array = xr.DataArray(t2m, dims=('month', 'latitude', 'longitude'))

# 将 DataArray 转换为 Dataset
dataset = data_array.to_dataset(name='t2m')

# 将 Dataset 写入为 .nc 文件
dataset.to_netcdf('t2m_example.nc')

这样就可以实现从ndarray到.nc文件的输出了，对于需要大量运算的研究来说，建议通过这种方法输出一些中间变量。

接下来读取这个文件：

1
f = xr.open_dataset("t2m_example.nc")

3.2 HDF5文件读取

HDF5格式文件的区域与netCDF文件的读取方法类似：

1
f = xr.open_dataset("file.h5")

我没有相关的HDF5数据文件，这里就不具体展示了。本质上没有任何区别，读进来以后都是DataSet结构。

如果报错，请查看是否缺少h5netcdf库，可通过：

1
conda install h5netcdf

进行安装。

3.3 tif文件读取

tif格式数据常见于各类地理遥感类数据中，如人口，灯光，土地，农作物，地形等等。

读取这类数据需要安装rioxarray库：

1
2
3
pip install rioxarray
#或
conda install -c conda-forge rioxarray

使用方法类似xarray：

1
2
import rioxarray
rds = rioxarray.open_rasterio("example.tif")

读取结果为DataArray格式，后续操作即同netCDF数据一致。

3.4 grib文件读取

需要安装cfgrib库：

1
conda install -c conda-forge cfgrib

之后通过xarray调用cfgrib引擎即可读取：

1
f = xr.open_dataset("example.grib", engine="cfgrib")

后续操作即同netCDF数据一致。

4 索引和筛选数据

索引实际上是Xarray中最重要的部分之一，精通索引可以通过简单的几个Xarray指令便实现用普通代码很复杂的功能（例如闰年2月29号的处理，如何筛选指定位置数据等等）。

4.1 位置索引（下标索引）

使用index来进行索引（我喜欢叫这种索引下标索引，就是用它下标的数字位置来进行索引）的方法与Numpy中是一样的，例如：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
print(temperature_da1[:2])
# <xarray.DataArray (time: 2, city: 3)>
# array([[0.33777998, 0.23731543, 0.74981315],
#        [0.4033871 , 0.05466693, 0.80505634]])
# Coordinates:
#   * time     (time) datetime64[ns] 1979-01-01 1979-01-02
#   * city     (city) <U9 'Beijing' 'Shanghai' 'Guangzhou'

print(temperature_da1[1,2])
#<xarray.DataArray ()>
#array(0.80505634)
#Coordinates:
#    time     datetime64[ns] 1979-01-02
#    city     <U9 'Guangzhou'

print(temperature_da1[:,[2,1]])
#<xarray.DataArray (time: 4, city: 2)>
#array([[0.74981315, 0.23731543],
#       [0.80505634, 0.05466693],
#       [0.66774707, 0.28434184],
#       [0.74922789, 0.26928012]])
#Coordinates:
#  * time     (time) datetime64[ns] 1979-01-01 1979-01-02 1979-01-03 1979-01-04
#  * city     (city) <U9 'Guangzhou' 'Shanghai'

数据数组的属性在所有索引操作中都会保留。

4.2 标签索引

Xarray 还支持基于标签的索引，就像 pandas 一样。

1
2
3
4
5
6
print(temperature_da1.loc['1979-01-02','Shanghai'])
#<xarray.DataArray ()>
#array(0.05466693)
#Coordinates:
#    time     datetime64[ns] 1979-01-02
#    city     <U9 'Shanghai'

这种操作会广泛应用于气象编程中，例如索引出某个区域的数据，某个时间断的数据等等，这种索引是十分灵活的。

1
2
3
# 以一个NCEP数据为例
f_u = xr.open_dataset('uwnd_1948_2022.nc')
u200 = f_u['uwnd'].loc[f_u.time.dt.month.isin([6,7,8])].loc['1979-05-01':'2021-09-30',200]

这里选出了月份在6,7,8月的200hPa风速数据。还可以更复杂：

1
u200 = f_u['uwnd'].loc[(f_u.time.dt.month.isin([12,1,2]))].loc['1979-12-01':'2021-02-28']

这样就选出了1979年到2020年冬季12,1,2月的数据，并且考虑到了冬季数据的跨年问题。

我先前的博客也有写到用索引方式剔除闰年2月29号的方式：

1
u200 = f_u['uwnd'].loc[~((f_u.time.dt.month==2)&(f_u.time.dt.day==29))]

~符号表示非，所以这里是选择了非2月29号的全部数据。

因此这种索引方式有很高的上限，取决于对基础编程的掌握水平。

4.3 模糊索引

如果我们想索引某个坐标点的数据，但是这个坐标在数据中又不存在（也可能是显示问题导致数据的小数位显示不全），比如我们索引300hPa数据，但实际上在数据中为300.1hPa，那么一般索引会报错。因此可以使用模糊索引：

1
2
3
4
5
6
7
8
9
10
11
da = xr.DataArray([1, 2, 3], [("x", [0, 1, 2])])
#<xarray.DataArray (x: 3)>
#array([1, 2, 3])
#Coordinates:
#  * x        (x) int64 0 1 2

da.sel(x=[1.1, 1.9], method="nearest")
#<xarray.DataArray (x: 2)>
#array([2, 3])
#Coordinates:
#  * x        (x) int64 1 2

还可以使用tolerance参数来限制模糊索引的有效距离：

1
2
3
4
5
print(da.reindex(x=[1.1, 1.5], method="nearest", tolerance=0.2))
#<xarray.DataArray (x: 2)>
#array([ 2., nan])
#Coordinates:
#  * x        (x) float64 1.1 1.5

4.4 where方法

数据数组也可以使用类似np中的where方法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
da = xr.DataArray(np.arange(16).reshape(4, 4), dims=["x", "y"])
#<xarray.DataArray (x: 4, y: 4)>
#array([[ 0,  1,  2,  3],
#       [ 4,  5,  6,  7],
#       [ 8,  9, 10, 11],
#       [12, 13, 14, 15]])
#Dimensions without coordinates: x, y

da.where(da>4,0)
#<xarray.DataArray (x: 4, y: 4)>
#array([[ 0,  0,  0,  0],
#       [ 0,  5,  6,  7],
#       [ 8,  9, 10, 11],
#       [12, 13, 14, 15]])
#Dimensions without coordinates: x, y

使用这个方法可以进行一些索引替换赋值等等的灵活操作。

4.5 气象实例

这里通过一个实例来展示这几种索引方法在实际科研中的常见用法(数据下载)：

1
2
3
4
5
6
7
8
9
10
11
import xarray as xr
f = xr.open_dataset('/mnt/f/NCEP1/monthly/surface/slp.mon.mean.nc')
print(f)
#<xarray.Dataset>
#Dimensions:  (lat: 73, lon: 144, time: 897)
#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
#  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 ... 2022-09-01
#Data variables:
#    slp      (time, lat, lon) float32 ...
  1. 通过下标索引：假设我们想选出前10个时间记录的数据：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
f.slp[:10]
#<xarray.DataArray 'slp' (time: 10, lat: 73, lon: 144)>
#[105120 values with 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
#  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 ... 1948-10-01
#Attributes:
#    long_name:     Sea Level Pressure
#    valid_range:   [ 870. 1150.]
#    units:         millibars
#    precision:     1
#    var_desc:      Sea Level Pressure
#    level_desc:    Sea Level
#    statistic:     Mean
#    parent_stat:   Other
#    dataset:       NCEP Reanalysis Derived Products
#    actual_range:  [ 955.56085 1082.5582 ]
  1. 通过坐标索引：假设我们想选出1979年，北半球的数据：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
f.slp.loc[f.time.dt.year==1979,90:0,:]
#<xarray.DataArray 'slp' (time: 12, lat: 37, lon: 144)>
#[63936 values with dtype=float32]
#Coordinates:
#  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... 10.0 7.5 5.0 2.5 0.0
#  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
#  * time     (time) datetime64[ns] 1979-01-01 1979-02-01 ... 1979-12-01
#Attributes:
#    long_name:     Sea Level Pressure
#    valid_range:   [ 870. 1150.]
#    units:         millibars
#    precision:     1
#    var_desc:      Sea Level Pressure
#    level_desc:    Sea Level
#    statistic:     Mean
#    parent_stat:   Other
#    dataset:       NCEP Reanalysis Derived Products
#    actual_range:  [ 955.56085 1082.5582 ]
  1. 模糊索引：我们想选出离37.2°E，44.6°N附近的数据，但是这个坐标点不在数据中存在：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
f.slp.sel(lat=44.6,lon=37.2,method='nearest')
#<xarray.DataArray 'slp' (time: 897)>
#[897 values with dtype=float32]
#Coordinates:
#    lat      float32 45.0
#    lon      float32 37.5
#  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 ... 2022-09-01
#Attributes:
#    long_name:     Sea Level Pressure
#    valid_range:   [ 870. 1150.]
#    units:         millibars
#    precision:     1
#    var_desc:      Sea Level Pressure
#    level_desc:    Sea Level
#    statistic:     Mean
#    parent_stat:   Other
#    dataset:       NCEP Reanalysis Derived Products
#    actual_range:  [ 955.56085 1082.5582 ]

这些索引方法可以组合起来使用：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
f.slp.sel(lat=44.6,lon=37.2,method='nearest')[f.time.dt.year==1979]
#<xarray.DataArray 'slp' (time: 12)>
#[12 values with dtype=float32]
#Coordinates:
#    lat      float32 45.0
#    lon      float32 37.5
#  * time     (time) datetime64[ns] 1979-01-01 1979-02-01 ... 1979-12-01
#Attributes:
#    long_name:     Sea Level Pressure
#    valid_range:   [ 870. 1150.]
#    units:         millibars
#    precision:     1
#    var_desc:      Sea Level Pressure
#    level_desc:    Sea Level
#    statistic:     Mean
#    parent_stat:   Other
#    dataset:       NCEP Reanalysis Derived Products
#    actual_range:  [ 955.56085 1082.5582 ]

所以说Xarray的索引方法十分灵活，兼容性很强，发挥程度完全取决于你的python基础（所以还是那句话，不要想着速成这个速成那个，欠的债都要还的）。