Python气象数据处理进阶之Xarray(4)：计算（完结）
2020-04-12 本文已影响0人
摸鱼咯
最近开始忙了。更新的速度会放缓一些，见谅。
这部分主要讲一下Xarray内部的计算函数，可以直接对dataarray进行操作。
还是以举例子的方式来展示。
先创建一个dataarray
import xarray as xr
da = xr.DataArray(np.random.RandomState(0).randn(2, 3),[('x', ['a', 'b']), ('y', [10, 20, 30])])
print(da)
#<xarray.DataArray (x: 2, y: 3)>
#array([[ 1.76405235, 0.40015721, 0.97873798],
# [ 2.2408932 , 1.86755799, -0.97727788]])
#Coordinates:
# * x (x) <U1 'a' 'b'
# * y (y) int64 10 20 30
ex1 基本运算
首先是基本的加减乘除(以减法为例)
print(da - 3)
#<xarray.DataArray (x: 2, y: 3)>
#array([[-1.23594765, -2.59984279, -2.02126202],
# [-0.7591068 , -1.13244201, -3.97727788]])
#Coordinates:
# * x (x) <U1 'a' 'b'
# * y (y) int64 10 20 30
同时还支持python自身的一些函数
print(abs(da))
#<xarray.DataArray (x: 2, y: 3)>
#array([[1.76405235, 0.40015721, 0.97873798],
# [2.2408932 , 1.86755799, 0.97727788]])
#Coordinates:
# * x (x) <U1 'a' 'b'
# * y (y) int64 10 20 30
也支持numpy和scipy的一些计算函数（但并不是所有，很多涉及到维度或者什么的就不能了）
print(np.sin(arr))
#<xarray.DataArray (x: 2, y: 3)>
#array([[ 0.9813841 , 0.38956314, 0.82979374],
# [ 0.78376151, 0.95628847, -0.82897801]])
#Coordinates:
# * x (x) <U1 'a' 'b'
# * y (y) int64 10 20 30
可以使用where函数进行条件更改
print(xr.where(da > 0, 'positive', 'negative'))
#<xarray.DataArray (x: 2, y: 3)>
#array([['positive', 'positive', 'positive'],
# ['positive', 'positive', 'negative']], dtype='<U8')
#Coordinates:
# * x (x) <U1 'a' 'b'
# * y (y) int64 10 20 30
使用@进行矩阵乘法，使用round进行四舍五入，使用.T数组转置等等。这部分略去代码了，可以自己试验
ex2 缺测处理
这部分是针对缺测的多种处理方法：
首先建立一个带有缺测的dataarray
x = xr.DataArray([0, 1, np.nan, np.nan, 2], dims=['x'])
print(x)
#<xarray.DataArray (x: 5)>
#array([ 0., 1., nan, nan, 2.])
#Dimensions without coordinates: x
判断是否缺测,还有notnull()函数
print(x.isnull())
#<xarray.DataArray (x: 5)>
#array([False, False, True, True, False])
#Dimensions without coordinates: x
判断有几个未缺测
print(x.count())
#<xarray.DataArray ()>
#array(3)
去除缺测
print(x.dropna(dim='x'))
#<xarray.DataArray (x: 3)>
#array([0., 1., 2.])
#Dimensions without coordinates: x
将缺测填为-1
print(x.fillna(-1))
#<xarray.DataArray (x: 5)>
#array([ 0., 1., -1., -1., 2.])
#Dimensions without coordinates: x
还有一种方法是上一篇文章将的插值，将缺测插补。
ex3 改变shape的运算
当涉及到运算后维度shape发生变化后，xarray提供了自己的一些方法。
先建立数组：
da = xr.DataArray(np.random.RandomState(0).randn(2, 3),[('x', ['a', 'b']), ('y', [10, 20, 30])])
print(da)
#<xarray.DataArray (x: 2, y: 3)>
#array([[ 1.76405235, 0.40015721, 0.97873798],
# [ 2.2408932 , 1.86755799, -0.97727788]])
#Coordinates:
# * x (x) <U1 'a' 'b'
# * y (y) int64 10 20 30
如果直接当作np数组进行sum,,mean等运算会报错，因为运算后数据结构发生了变化。
正确的操作是：
print(da.sum(dim='x'))
#<xarray.DataArray (y: 3)>
#array([4.00494555e+00, 2.26771520e+00, 1.46010423e-03])
#Coordinates:
# * y (y) int64 10 20 30
要指明求和的方式，而不是类似np.sum()，直接填充第几个轴。
特别需要注意的是求平均时缺测的处理
print(xr.DataArray([1, 2, np.nan, 3]).mean())
#<xarray.DataArray ()>
#array(2.)
print(xr.DataArray([1, 2, np.nan, 3]).mean(skipna=False))
#<xarray.DataArray ()>
#array(nan)
默认是跳过nan的
ex4 滑动运算
这部分主要是计算滑动平均或者滑动标准差的，这部分比较抽象，我尽可能表述的通俗一些。
新建数组：
da = xr.DataArray(np.arange(0, 7.5, 0.5).reshape(3, 5),dims=('x', 'y'))
print(da)
#<xarray.DataArray (x: 3, y: 5)>
#array([[0. , 0.5, 1. , 1.5, 2. ],
# [2.5, 3. , 3.5, 4. , 4.5],
# [5. , 5.5, 6. , 6.5, 7. ]])
#Dimensions without coordinates: x, y
通过指令进行滑动平均，mean换位std即为滑动标准差：
print(da.rolling(y=3).mean())
#<xarray.DataArray (x: 3, y: 5)>
#array([[nan, nan, 0.5, 1. , 1.5],
# [nan, nan, 3. , 3.5, 4. ],
# [nan, nan, 5.5, 6. , 6.5]])
#Dimensions without coordinates: x, y
print(da.rolling(y=3, center=True).mean())
#<xarray.DataArray (x: 3, y: 5)>
#array([[nan, 0.5, 1. , 1.5, nan],
# [nan, 3. , 3.5, 4. , nan],
# [nan, 5.5, 6. , 6.5, nan]])
#Dimensions without coordinates: x, y
这指的是对y轴（即列）取滑动步长为3进行平均，注意这个xr.rolling和np.roll是不一样的。xarray的rolling函数是指定滑动步长，即窗口，而不是指数组滚动。
以第一行数据为例子，原数组第一行为【0 , 0.5, 1 , 1.5, 2】
滑动后为[(0+0.5+1)/3=0.5，(0.5+1+1.5)/3=1，(1+1.5+2)/3=1.5]
而np.roll的效果是变成[1，1.5，2，0，0.5]
但是这种滑动平均的缺点是两侧变成了缺测。
xarray给出了一个参数可以规定最短窗口，也就是说在数据两端，窗口不足时，可以减小窗口，延长滑动后的结果。
print(da.rolling(y=3, center=True, min_periods=2).mean())
#<xarray.DataArray (x: 3, y: 5)>
#array([[0.25, 0.5 , 1. , 1.5 , 1.75],
# [2.75, 3. , 3.5 , 4. , 4.25],
# [5.25, 5.5 , 6. , 6.5 , 6.75]])
#Dimensions without coordinates: x, y
实现了滑动后长度不改变，但是两侧数据是存在误差的。
ex5 权重运算
依然还是新建一组数据：
coords = dict(month=('month', [1, 2, 3]))
print(coords)
#{'month': ('month', [1, 2, 3])}
prec = xr.DataArray([1.1, 1.0, 0.9], dims=('month', ), coords=coords)
print(prec)
#<xarray.DataArray (month: 3)>
#array([1.1, 1. , 0.9])
#Coordinates:
# * month (month) int64 1 2 3
weights = xr.DataArray([31, 28, 31], dims=('month', ), coords=coords)
print(weights)
#<xarray.DataArray (month: 3)>
#array([31, 28, 31])
#Coordinates:
# * month (month) int64 1 2 3
prec比如是降水数据，weight是权重数组
对prec进行权重处理：
#这两种方法等价
weighted_prec = prec.weighted(weights)
#weighted_sum = (prec * weights).sum()
加权求和：
print(weighted_prec.sum())
#<xarray.DataArray ()>
#array(90.)
结果相当于311.1+281+31*0.9
ex6 聚块
新建数组：
import pandas as pd
x = np.linspace(0, 10, 300)
t = pd.date_range('15/12/1999', periods=364)
da = xr.DataArray(np.sin(x) * np.cos(np.linspace(0, 1, 364)[:, np.newaxis]),
dims=['time', 'x'], coords={'time': t, 'x': x})
print(da)
#<xarray.DataArray (time: 364, x: 300)>
#array([[ 0. , 0.03343858, 0.06683976, ..., -0.48672119,
# -0.51565952, -0.54402111],
# [ 0. , 0.03343845, 0.06683951, ..., -0.48671934,
# -0.51565756, -0.54401905],
# [ 0. , 0.03343807, 0.06683875, ..., -0.4867138 ,
# -0.51565169, -0.54401285],
# ...,
# [ 0. , 0.0182217 , 0.03642301, ..., -0.26522911,
# -0.28099849, -0.29645358],
# [ 0. , 0.01814439, 0.03626848, ..., -0.26410385,
# -0.27980632, -0.29519584],
# [ 0. , 0.01806694, 0.03611368, ..., -0.26297658,
# -0.27861203, -0.29393586]])
#Coordinates:
# * time (time) datetime64[ns] 1999-12-15 1999-12-16 ... 2000-12-12
# * x (x) float64 0.0 0.03344 0.06689 0.1003 ... 9.9 9.933 9.967 10.0
比如说气象上需要求逐周，逐旬，逐候的数据，处理起来很麻烦，我们则可以使用这个方法先将其聚块，再进行平均。
print(da.coarsen(time=7).mean())
#<xarray.DataArray (time: 52, x: 300)>
#array([[ 0. , 0.03343693, 0.06683647, ..., -0.48669718,
# -0.51563408, -0.54399428],
# [ 0. , 0.03342539, 0.06681339, ..., -0.48652913,
# -0.51545604, -0.54380644],
# [ 0. , 0.03340141, 0.06676547, ..., -0.48618016,
# -0.51508632, -0.54341639],
# ...,
# [ 0. , 0.0193641 , 0.03870654, ..., -0.28185754,
# -0.29861556, -0.31503961],
# [ 0. , 0.01883484, 0.03764862, ..., -0.2741539 ,
# -0.29045391, -0.30642905],
# [ 0. , 0.01829859, 0.03657671, ..., -0.26634833,
# -0.28218424, -0.29770455]])
#Coordinates:
# * time (time) datetime64[ns] 1999-12-18 1999-12-25 ... 2000-12-09
# * x (x) float64 0.0 0.03344 0.06689 0.1003 ... 9.9 9.933 9.967 10.0
可以看到，coarsen(time=7)将时间轴每7天固定在一起，通过平均函数，得到了逐周的平均结果。
但是这里364是可以整除7的，如果不能刚好整除的话，还提供了一些参数来处理边界结果。
boundary='trim'，或者boundary='pad'.
trim是指末尾多余的不参与计算，或者叫修建多余条目。
pad指填充nan到不足的地方，补全为可以整除的更大数据再进行聚块。
还可以通过coord_func函数返回不同的坐标信息。
比如说默认下，7天平均，结果的时间是第4天，就是中间那天，但是通过coord_func={'time': 'min'}可以指定时间为最小的那天。
ex6 差分和积分计算
新建数据：
da = xr.DataArray(np.arange(8).reshape(4, 2), dims=['x', 'y'],coords={'x': [0.1, 0.11, 0.2, 0.3]})
print(da)
#<xarray.DataArray (x: 4, y: 2)>
#array([[0, 1],
# [2, 3],
# [4, 5],
# [6, 7]])
#Coordinates:
# * x (x) float64 0.1 0.11 0.2 0.3
#Dimensions without coordinates: y
中央差计算差分：
print(da.differentiate('x'))
#<xarray.DataArray (x: 4, y: 2)>
#array([[200. , 200. ],
# [182.22222222, 182.22222222],
# [ 21.16959064, 21.16959064],
# [ 20. , 20. ]])
#Coordinates:
# * x (x) float64 0.1 0.11 0.2 0.3
#Dimensions without coordinates: y
这里补充说一句，那些可选参数大部分是通用的，你可以通过参数的设置达到不同的目的。
梯形原则计算积分：
print(da.integrate('x'))
#<xarray.DataArray (y: 2)>
#array([0.78, 0.98])
#Dimensions without coordinates: y
ex7 广播计算
广播的意思大致可以理解为扩展，主要是用于两个不同大小(指shape的不同)数组的计算，原理同线性代数的矩阵计算。并不是所有不同大小的数组都可以计算，广播也要遵循一定的规则。
我们先建立两个大小不等的数组：
a = xr.DataArray([1, 2], [('x', ['a', 'b'])])
b = xr.DataArray([-1, -2, -3], [('y', [10, 20, 30])])
print(a)
#<xarray.DataArray (x: 2)>
#array([1, 2])
#Coordinates:
# * x (x) <U1 'a' 'b'
print(b)
#<xarray.DataArray (y: 3)>
#array([-1, -2, -3])
#Coordinates:
# * y (y) int64 10 20 30
我们可以计算矩阵乘法：
print(a * b)
#<xarray.DataArray (x: 2, y: 3)>
#array([[-1, -2, -3],
# [-2, -4, -6]])
#Coordinates:
# * x (x) <U1 'a' 'b'
# * y (y) int64 10 20 30
实际上原理是这样的。a(n,m) * b(m) = a(n,m) * b(1,m) b默认扩充了，也就是广播了。
a2, b2 = xr.broadcast(a, b)
print(a2,b2)
#<xarray.DataArray (x: 2, y: 3)>
#array([[1, 1, 1],
# [2, 2, 2]])
#Coordinates:
# * x (x) <U1 'a' 'b'
# * y (y) int64 10 20 30 <xarray.DataArray (x: 2, y: 3)>
#array([[-1, -2, -3],
# [-1, -2, -3]])
#Coordinates:
# * y (y) int64 10 20 30
# * x (x) <U1 'a' 'b'
通过broadcast函数可以更清晰的看出广播的原理。
