使用Xarray库处理nc数据

最近在做一些数据方面的研究，需要处理很多 netCDF 格式（例如 nc 格式、nc4 格式）的栅格数据。

处理这类数据，使用 Python 应该比较简单，虽然 Python 有性能的问题，但是凭借着众多的数据处理库，还是使用 Python 处理 netCDF 数据较为方便。

其中 netCDF4 库是处理 nc 和 nc4 格式数据的基本库，而 Xarray 库，则是整合了 netCDF4、NumPy 等常用库，称为了处理 netCDF 格式最方便的库了。接下来介绍下 Xarray 库的安装、使用以及踩坑的经验。

1.Xarray 库

1.1 Xarray 简介

xarray 库应该是目前用于处理 netCDF 数据最为方便的库了。xarray 可以看作是对 numpy 和 pandas 的扩展，同时也借鉴了许多 pandas 的优点，具有良好的易用性。xarray 与 netCDF4 在使用上存在许多相似之处，但是从输出的结果我们可以发现，同样是 Dataset，但是 xarray 的结果显然有更好的可读性。

1.2 netCDF 数据结构

netCDF 数据是一种多维数据，包含经度、维度、时间轴、气温、太阳辐射等各类数据。

比较好理解的方式是，将经度、维度和时间想象为自变量；而气温、太阳辐射等各类数据想象为因变量。

加入我们想要知道某一天北京市的温度，那么就可以：

f(经度, 维度, 时间) = 北京市某一天的气温

确定了相关的自变量，就等得出因变量的值。

1.3 Xarray 数据结构

xarray 有两个核心数据结构：DataArray 和 Dataset。

DataArray（一个带有标签的多维数组）
- values 获取数组的具体数值
- dims 获取维度的名字，如 ('x', 'y', 'z')
- coords 获取一个类似于字典的结果，里面包含各个坐标
- attrs 获取原始数据的属性，比如变量的名字、单位等
Dataset（可以简单的理解为由多个 DataArray 组成的集合）
- dims 获取维度的名字，结果类似于字典，如 {'x': 6, 'y': 6, 'time': 8}
- data_vars 获取物理量的名字
- coords 获取一个类似于字典的结果，里面包含各个坐标
- attrs 获取原始数据的属性，比如变量的名字、单位等

Data variables 变量，对应着真实的物理数据。比如在气象学中要做出一个气压图，也就是获取”东经 xx 度、北纬 yy 度的大气压为多少帕”，显然这就是一个二维单值函数，两个维度分别是经度和纬度，函数值为大气压

1.4 Xarray 安装

根据 Xarray 的官方文档，需要以下必要的依赖：

Python (3.8 or later)
numpy (1.18 or later)
packaging (20.0 or later)
pandas (1.1 or later)
netCDF4

根据我这段时间使用的经验来看，建议看装以下的依赖，这些依赖较为常用，剩下的就根据自己的数据处理情况自行选择安装：

pip install numpy
pip install pandas
pip install packaging
pip install netcdf4
pip install xarray
pip install dask
pip install scipy

1.5 读取数据

按需引入头部以及它们的常用缩写：

import numpy as np
import pandas as pd
import netCDF4 as nc
from netCDF4 import Dataset
import xarray as xr

加载 nc 文件（需要注意一些问题）：

# 打开 nc 文件
ds = xr.open_dataset("filename.nc")
print(ds)

# 读取多文件数据集（此函数自动将多个文件连接并合并到一个 xarray 数据集中。推荐使用 xarray 打开多个文件的方法）
# 需要安装 dask 库
xr.open_mfdataset('myfiles/*.nc', parallel=True)

💡 中文路径和文件名不能出现中文，否则 xarray 库会报错。（如果要支持中文的话，就要转而使用 netCDF4 库）。填写路径的时候， \ 为了防止转义，要使用双 \\ 才可以。

其中，读取数据使用的方法有两种，也就是上面将的两种数据结构方法。这两种 Xarray 数据结构的区别是：

xarray.open_dataset() 读取 Dataset 类型数据，即能读取多个物理量。
xarray.open_dataarray() 读取 DataArray 类型数据，即只能读取单个物理量。

💡 如果 nc 文件中含有多个物理量，用 open_dataarray() 读取会报错，因此建议统一都用 open_dataset() 来读取文件。

从文件中读取数据，假如数据中含有一个名为 longitude 的物理量：

ds = xarray.open_dataset()
ds.longitude 或 ds[longitude]

1.6 索引核心方法

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

1.6.1 根据位置索引

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

1.6.2 根据维度名字索引

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

# 通过数字索引
ds.isel(longitude=1, time=0)

# 通过标签索引
ds.sel(longitude=0.75, time='2018-01-01')

💡 如果遇到 KeyError(key) from err 错误。可能是存在小浮点差异，请尝试通过 method=”nearest” 解决。

1.7 插值

Xarray 库还带有多种插值功能，这对于我处理栅格数据并进行插值，就非常的方便。xarray 中对 scipy 的插值函数进行了进一步的封装，可以让我们方便的调用。只需要对 DataArray 或 DataSet 使用 interp() 函数就可以实现插值了，就像索引一样简单。不管是一维数据还是多维数据都可以轻松搞定。

In [1]: da = xr.DataArray(np.sin(0.3 * np.arange(12).reshape(4, 3)),
   ...:                   [('time', np.arange(4)),
   ...:                    ('space', [0.1, 0.2, 0.3])])
   ...: 

# label lookup
In [2]: da.sel(time=3)
Out[2]: 
<xarray.DataArray (space: 3)>
array([ 0.42738 ,  0.14112 , -0.157746])
Coordinates:
    time     int64 3
  * space    (space) float64 0.1 0.2 0.3

# interpolation
In [3]: da.interp(time=2.5)
Out[3]: 
<xarray.DataArray (space: 3)>
array([0.700614, 0.502165, 0.258859])
Coordinates:
  * space    (space) float64 0.1 0.2 0.3
    time     float64 2.5

da.sel(time=3)，索引时间维的值为 3（12行）。
da.interp(time=2.5)，将时间维原本不存在的 2.5 插值了出来（22行）。

数据合并也支持维度的拼接和变量的合并，具体参阅 Xarray 库的说明即可。

2.Xarray 库使用例子

最后给大家一个使用用例，方便大家入门理解。

本例采用 “全球高分辨率（3小时，10公里）地表太阳辐射数据集（1983-2018）” 的部分数据作为演示，大家可以自行下载部分数据进行测试。

Github 仓库：https://github.com/buildingdata/python-xarray

阅读: 1,668