数据处理 | xarray的计算距平、重采样、时间窗

点击下方公众号，回复资料，收获惊喜

引入相关包和导入数据：

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

# 数据导入
path = "...\\sst.mnmean.nc"
# 丢弃一个不必要导入的变量
ds = xr.open_dataset(path, drop_variables=["time_bnds"])
ds = ds.sel(time=slice("1960", "2018")).load()

背景知识：距平

现代气候学认为在相当长的时间段（一般认为是 30 年）中，变量多年平均是一个稳定的值。因此在一个时间段中，如果能够充分认识变量随平均状态的变化趋势，那么对于预测未来情况是非常有利的。那么这个所谓随着平均态的偏移值便可称为距平（异常，anomaly）.

距平

下面便提出一个问题：为什么要费尽心思研究变量的距平而非变量的原始数据？若针对于温度这个变量而言，即为什么要使用温度距平（偏离平均值的值）而不非研究绝对温度的变化？

出于以下几个原因，很难对全球平均表面温度以绝对温度的形式进行计算。

某些地域的气象观测站点分布稀少（如撒哈拉沙漠地区、偏远的密林），这就意味着为取得格点数据（栅格数据）必须对离散的站点数据值在较大且站点分布稀疏区域内进行插值。这会带来很大的数据不真实性。

对于那些山区中的数据（山区中的的气象观测大多是有人居住地区），必须考虑海拔高度对区域平均温度的影响。例如，对于一个地区的夏季而言，无论是在山顶还是山下，都可能比往年的平均温度低，然而若考虑绝对温度，这两个地方有很大的不同（一般认为山顶气温比山下温度低）。

在同一时间范围内在一个更小的尺度下（即格点分辨率）考虑变量变化的基准参考值，然后基于这个基准参考值（多年平均值）计算相对于这个基准参考值的异常变化（距平）。在这种情况下，整合了数据，使得不同地域的变量能够得以进行比较，以便反映一个区域内不同地方的变量分布形式。

来源：https://www.ncdc.noaa.gov/monitoring-references/faq/anomalies.php

Groupby（Ⅲ）

Transformations 转换

下面需从数据集中删除气候平均，从而得到变量随气候平均态变化的残差。一般将这个残差称为距平。

对转换（Transformations）操作而言，消除数据的气候平均是一个很好的例子。转换操作对分组的对象进行操作，但不改变原数据的维度尺寸。

xarray 通过使用Groupby 算法使这些类型的转换变得容易。下面给出了计算去除月份温度差异的海温月数据。

def remove_time_mean(x):
    return x - x.mean(dim="time")

ds_anom = ds.groupby("time.month").map(remove_time_mean)
ds_anom

也可以简写为下面这种形式

gb = ds.groupby("time.month")
ds_anom = gb - gb.mean(dim="time")
ds_anom

ds_anom

gb是分好月份后的海温数据（12 组），gb.mean(dim="time")是各月的平均海温（12 组），那么gb - gb.mean(dim="time")即为对 12 组中的对应组的海温数据（这个组内的每一天的海温数据）减去平均的海温数据。这个结果即为距平。

当经过上述去除季节性周期的影响后，便很容易发现气候变率的信号。

北大西洋单点的时间序列

ds_anom.sst.sel(lon=300, lat=50).plot()

北大西洋单点的时间序列

2018 年 1 月 1 日与 1960 年 1 月 1 日之间 SST 之间的差异

(ds_anom.sel(time="2018-01-01") - ds_anom.sel(time="1960-01-01")).sst.plot()

2018年1月1日与1960年1月1日之间SST之间的差异

Resample（重采样）

xarray 中的Resample（重采样）的处理方法与 Pandas 包几乎相同。就本质而言，Resample 也是一个分割数据的操作。它与分割操作的基本语法类似。应当注意，对于 Resample 操作而言，其作用对象必须是时间维度。

为说明 Resample 的用法，下面给出一个例子计算逐五年的平均值曲线。

resample_obj = ds_anom.resample(time="5Y")
resample_obj

resample_obj

可以看到对于 Resample 操作而言，与 Groupby 操作非常类似，首先也创建了一个DatasetResample对象。.resample(time="5Y")是对如何对时间进行重采样进行设置，维度为time，设置的时间间隔为 5 年。应当指出这里的时间间隔写法与之前pd.date_range函数中的freq的时间间隔的关键词是一致的。

ds_anom_resample = resample_obj.mean(dim="time")
ds_anom_resample

ds_anom_resample

之后就需要对这些分割好的 Resample 对象进行取平均，以便获得每一个分组好的 Resample 对象中的平均值。

假如第一个 Resample 对象的时间范围为 2010 年-2014 年，那么需要对这五年进行平均后，以便得到第一个进行重采样后的值。往后的时间范围类似。

为了说明进行重采样后的效果，下面来看一下(50°N, 60°E)的海温变化情况

ds_anom.sst.sel(lon=300, lat=50).plot()
ds_anom_resample.sst.sel(lon=300, lat=50).plot(marker="o")

(50°N, 60°E) 的海温变化

第一行代码将原始海温变化的时间序列画了出来，第二行画了经逐 5 年平均后的海温变化的时间序列。

.plot(marker="o")中的参数marker="o"指定了点的记号类型。如下图[1]是关于marker参数的可设置类型

matplotlib.markers

注意：resample 仅能用于正确的日期、时间索引。

Rolling（时间窗移动）

Pandas Rolling (Source: forgifs.com)

Rolling 方法也与pandas 包[2]中的类似，但是稍有不同的是，它可适用于任意维度。如果将其作用于时间维度，也可称之为滑动平均。

ds_anom_rolling = ds_anom.rolling(time=12, center = True).mean()
ds_anom_rolling

ds_anom_rolling

参数time=12指定了对维度time以 12 个月为周期（月数据）变动时间窗，center参数表明以当前窗的两侧筛选数据，否则是以当前窗的前 12 个月作为筛选目标（包括本身）。.mean()表明对每一个 Rolling 对象取平均。

为了更好的说明 Rolling 的作用，下面举一个简单的例子说明其功能。

da = xr.DataArray(
     np.linspace(0, 11, num=12),
     dims="time",
     coords=[
     pd.date_range("1999-12-15",periods=12,freq=pd.DateOffset(months=1))
     ]
    )
da

da

此处创建 DataArray 类型 da 的方法与之前创建 DataArray 稍有不同。np.linspace(0, 11, num=12)代表创建数组的初始值为 0，终末值为 11，并且在这个范围内均匀间隔生成 12 个样本。关于这个函数的说明，可参考numpy.linspace[3].

dims的创建与之前的类似，但coords就有着明显的区别，此处的coords是一个元组列表（用方括号包裹，List），而之前的教程中创建的是一个字典（用花括号包裹，dict）. 两者创建的区别在于如果用列表创建 DataArray 的话，坐标名称和维度名称是重名的（Coordinates 项会加粗或者在名称前加*）。若要创建非索引坐标，则必须通过字典创建。

对于多个维度的创建，列表的创建方法也与之前的字典创建方法类似

foo = xr.DataArray(
   np.random.rand(4, 3),
   dims=("time", "space"),
   coords=[
       pd.date_range("2000-01-01", periods=4),
       ["IA", "IL", "IN"],
       ]
   )
foo

foo

多个维度dims需用小括号或者方括号包裹。不同的 coords 之间的参数用逗号间隔，因为用列表创建坐标维度的特性，无需写坐标维度名称。坐标维度的名称将沿用维度名称的名字。

首先先创建一个时间窗对象

da.rolling(time=5, center=True)

时间窗对象

然后对这个时间窗对象施以平均的方法，即.mean()

da.rolling(time=5, center=True).mean()

5年滑动平均

我们可以通过下图来理解时间窗是如何操作数据的

时间窗操作过程理解

若时间窗为偶数值，那么对应中心位置将会在平均位置偏右侧

da.rolling(time=4, center=True).mean()

可用下图理解上述原理

若不指定参数center=True，则采用从当前元素往上筛选的方法，否则采用以当前元素为中心，从两个方向上进行筛选。

da.rolling(time=5).mean()

da.rolling(time=5).mean()

具体操作如下图所示

当然和 grouby 对象类似，也可用 list 来访问每一个时间窗移动对象

rolling_obj = da.rolling(time=5, center=True)
list(rolling_obj)[4][1]

list(rolling_obj)[4][1]

关于 pandas 中 rolling 方法的深入理解可参见详解pandas 中的 rolling[4]

参考资料

[1]

下图: https://matplotlib.org/stable/api/markers_api.html?highlight=marker

[2]

pandas包: https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.rolling.html

[3]

numpy.linspace: https://numpy.org/doc/stable/reference/generated/numpy.linspace.html

[4]

pandas中的rolling: https://www.cnblogs.com/traditional/p/13776180.html

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

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