数据处理 | xarray的NC数据基础计算(1)

2021-04-22 10:31:30 浏览数 (1)
代码语言：javascript复制
import expectexception
# 若没有安装则需要在conda的base环境中运行下面的代码进行安装
# pip install ExpectException

import numpy as np
import xarray as xr
from matplotlib import pyplot as plt

示例数据

首先我们先导入所需的数据，本次使用的是经扩展重构的海表面温度 v5 数据集（Extended Reconstructed Sea Surface Temperature, abbr. ERSST）。这个数据集可追溯到 1854 年的海表面温度，并被广泛使用。

ERSST v5

下载完毕数据后，我们利用.open_dataset函数导入 NetCDF 数据

代码语言：javascript复制
path = "...\sst.mnmean.nc"
# 删除一些不必要的变量
ds = xr.open_dataset(path, drop_variables=["time_bnds"])
# 提取1960年~2018年的数据
ds = ds.sel(time=slice("1960", "2018")).load()
ds

ds

下面我们来做一下数据的基本制图，通过图形来检查下载数据的正确性。

代码语言：javascript复制
ds.sst.isel(time=0).plot(vmin=-2, vmax=30)

ds.sst.isel(time=0).plot(vmin=-2, vmax=30)

上述代码选取了时间维度第一个的变量 sst，同时通过vmin和vmax定义色标的绘制变量数值范围为-2 至 30. 假如我们缩小这个范围，缩小至 20 至 30，我们可以得到以下的图像。

代码语言：javascript复制
ds.sst.isel(time=0).plot(vmin=20, vmax=30)

ds.sst.isel(time=0).plot(vmin=20, vmax=30)

基本计算

xarray 的 DataArray 和 DataSet 对象可以无缝地使用计算操作符（如 , -, *, /）和 numpy 数组函数。

下面我们将上述海温数据的摄氏温度改写为开尔文温度为例说明上述问题。

两者温度转化关系：开氏度 = 摄氏度 273.15

代码语言：javascript复制
sst_kelvin = ds.sst   273.15
sst_kelvin

sst_kelvin

可以发现再进行计算操作后，数据集的维度和坐标都没有发生变化。

需注意的是，许多导入的 xarray 数据集存在单位（units）属性，这些属性可用于绘图，目前独立于 xarray 项目进行开发的包pint[1]可以实现对单位的完全感知并进行转换。

下面我们来尝试一下用更为复杂的函数进行计算。

假设转为开尔文温度的公式如下所示

f(theta) = 0.5 ln (theta^2)

则可以编写以下代码

代码语言：javascript复制
f = 0.5 * np.log(sst_kelvin ** 2)
f

f

注意到**在 python 中代表乘方，此处** 2代表平方。

apply_ufunc 函数的使用

上面可以调用np.log（ds）并使其在 xarray 中“正常工作”是非常幸运的，因为并非所有的库都能直接在 xarray 中正常工作。

numpy相关的数学函数均可以直接在 xarray 中直接运算。

我们以一个实例来开始下面的内容：用于海水热力学领域的Gibbs 海水工具包[2]。这个包提供了对 numpy 数组进行操作的 ufunc。

代码语言：javascript复制
import gsw
# 若没有安装则需要在conda的base环境中运行下面的代码进行安装
# pip install gsw

比如我们需要进行将上述数据的 IPTS-68 温度转换为 ITS-90 温度。

国际计量委员会决定从 1990 年 1 月 1 日起采用新的温度标准——“1909 年国际温度刻度”(简称 ITS-90) , 这次温度标准修订的主要目的是为了解决以前采用的温度标准 IPTS-68 存在的问题, 使温度标准更准确地反映热力学温度, 提高温度标准的复现性。这次温度标准的修订, 不是温度单位的变更，是温度标准的实现法、维持法的变更。（新的国际温度标准 ITS-90 简介，戚栋，1991）

我们先来看一下gsw.t90_from_t68函数的介绍

代码语言：javascript复制
gsw.t90_from_t68?

gsw.t90_from_t68?

同样的内容也可在网站上查询到

代码语言：javascript复制
http://www.teos-10.org/pubs/gsw/html/gsw_t90_from_t68.html

类似于上面的np.log函数，我们可以直接将 xarray 的 DataArray 对象放在函数括号里。

代码语言：javascript复制
gsw.t90_from_t68(ds.sst)

gsw.t90_from_t68(ds.sst)

当然也可以使用xr.apply_ufunc函数对于数组中的每个元素进行gsw.t90_from_t68操作。这对于一些不能直接应用于 xarray 对象的函数是非常便捷的。

代码语言：javascript复制
xr.apply_ufunc(gsw.t90_from_t68, ds.sst)

xr.apply_ufunc(gsw.t90_from_t68, ds.sst)

apply_ufunc 函数功能强大，有很多可选参数以便进行复杂操作。更多可查阅Xarray docs[3]。

参考资料

[1]

pint: https://pint.readthedocs.io/en/0.16.1/

[2]

Gibbs海水工具包: https://teos-10.github.io/GSW-Python/

[3]

Xarray docs: http://xarray.pydata.org/en/latest/generated/xarray.apply_ufunc.html

注：本文基于Xarray Tutorials进行改写，遵循Apache-2.0 License

https://github.com/xarray-contrib/xarray-tutorial

numpy https 网络安全 编程算法 github

0 人点赞