参考:Interpolating data (pydata.org)
Python气象数据处理进阶之Xarray(3):插值 - 简书 (jianshu.com)
1. 标量和1维插值
对一个DataArray的插值有一些类似于对DataArray的索引
首先,我们建立一个DataArray:
da = xr.DataArray(np.sin(0.3 * np.arange(12).reshape(4,3)),
[('time', np.arange(4)), ('space', [0.1,0.2,0.3])])
print(da)
#xarray.DataArraytime: 4space: 3
#array([[ 0. , 0.29552021, 0.56464247],
# [ 0.78332691, 0.93203909, 0.99749499],
# [ 0.97384763, 0.86320937, 0.67546318],
# [ 0.42737988, 0.14112001, -0.15774569]])
#Coordinates:
#time (time) int64 0 1 2 3
#space (space) float64 0.1 0.2 0.3
#Attributes: (0)
- 我们已知的是第0,1,2,3天的数据,想得到第2.5天的数据。
print(da.sel(time=3))
print(da.interp(time=2.5))
#xarray.DataArray (space: 3)
#array([ 0.42737988, 0.14112001, -0.15774569])
#Coordinates:
#time () int64 3
#space (space) float64 0.1 0.2 0.3
#Attributes: (0)
#xarray.DataArray (space: 3)
#array([0.70061376, 0.50216469, 0.25885874])
#Coordinates:
#space (space) float64 0.1 0.2 0.3
#time () float64 2.5
#Attributes: (0)
- 与索引类似,
interp()
函数也可以接收一个数组形式的index,插值后的输出结果也为一个数组。
# lable look up
print(da.sel(time=[2,3]))
#xarray.DataArray (time: 2, space: 3)
#array([[ 0.97384763, 0.86320937, 0.67546318],
# [ 0.42737988, 0.14112001, -0.15774569]])
#Coordinates:
#time (time) int64 2 3
#space (space) float64 0.1 0.2 0.3
#Attributes: (0)
# interpolation
print(da.interp(time=[2.5, 3.5]))
#xarray.DataArray (time: 2, space: 3)
#array([[0.70061376, 0.50216469, 0.25885874],
# [ nan, nan, nan]])
#Coordinates:
#space (space) float64 0.1 0.2 0.3
#time (time) float64 2.5 3.5
#Attributes: (0)
这里看到,如果向外插,得到的就会是一些缺测值。
- 如果是对numpy.datetime64的时间维进行插值,可以输入字符串作为索引
da_dt64.interp(time = pd.date_range("1/1/2000", "1/3/2000")
。 - 也可以指定对多个维度插值
da.interp(time=[1.5, 2.5], space=[0.15,0.25])
。以此可以实现不同网格间的插值,甚至是站点数据向格点的插值。
2. 插值方法介绍
da = xr.DataArray(np.sin(np.linspace(0,2*np.pi,10)), dims="x", coords={"x":np.linspace(0,1,10)})
da.plot.line("o",label="oringinal")
da.interp(x=np.linspace(0,1,100)).plot.line(label="linear (default)")
da.interp(x=np.linspace(0,1,100), method="cubic").plot.line(label="cubic")
plt.legend()
通过其他的参数传递可以实现不同的插值方法,并对interp的结果进行调整。
比如,我们可以调整外推的结果:
## 用0填充外推结果中的缺测值
da.interp(x=np.linspace(-0.5, 1.5, 10), kwargs={'fill_value': 0.0})
#xarray.DataArray (x: 10)
#array([ 0. , 0. , 0. , 0.81379768, 0.60402277,
# -0.60402277, -0.81379768, 0. , 0. , 0. ])
#Coordinates:
#x (x) float64 -0.5 -0.2778 -0.05556 ... 1.278 1.5
#Attributes: (0)
## 外推
da.interp(x=np.linspace(-0.5, 1.5, 10), kwargs={"fill_value":"extrapolate"})
xarray.DataArray (x: 10)
#array([-2.89254424, -1.60696902, -0.3213938 , 0.81379768, 0.60402277,
# -0.60402277, -0.81379768, 0.3213938 , 1.60696902, 2.89254424])
#Coordinates:
#x (x) float64 -0.5 -0.2778 -0.05556 ... 1.278 1.5
#Attributes: (0)
3. 高端操作
上图中,图左是通过索引,得到一条红色的斜线,图右表示通过插值得到一条红色的斜线。
以下代码实现左图操作:
da = xr.DataArray(
np.sin(0.3 * np.arange(20).reshape(5, 4)),
[("x", np.arange(5)), ("y", [0.1, 0.2, 0.3, 0.4])],)
# advance indexing
x = xr.DataArray([0,2,4], dims = "z")
y = xr.DataArray([0.1,0.2,0.3], dims = "z")
da.sel(x=x,y=y)
#xarray.DataArray (z: 3)
#array([ 0. , 0.42737988, -0.77276449])
#Coordinates:
#x (z) int64 0 2 4
#y (z) float64 0.1 0.2 0.3
#Attributes: (0)
以下代码实现右图操作:
# advanced interpolation
x = xr.DataArray([0.5, 1.5, 2.5], dims = "z")
y = xr.DataArray([0.15, 0.25, 0.35], dims = "z")
da.interp(x=x,y=y)
#xarray.DataArray (z: 3)
#array([ 0.55626357, 0.63496063, -0.46643289])
#Coordinates:
#x (z) float64 0.5 1.5 2.5
#y (z) float64 0.15 0.25 0.35
#Attributes: (0)
4. 缺测值
最近用的服务器上没有安装ncl,为了插值一套数据被迫开始研究python中的插值方法,被缺测值折磨良久终于看到这里。下次一定谨记,遇事先看官方文档TAT。
我们先看默认方法插值出来的结果:
da = xr.DataArray([0,2,np.nan,3,3.25], dims = "x", coords = {"x": range(5)})
da.interp(x=[0.5, 1.5, 2.5])
#xarray.DataArray (x: 3)
#array([ 1., nan, nan])
#Coordinates:
#x (x) float64 0.5 1.5 2.5
#Attributes: (0)
可以看到,linear
方法会返回包含缺测值的数组(nearest
也是这样)。
da.interp(x=[0.5, 1.5, 2.5], method="cubic")
#xarray.DataArray (x: 3)
#array([nan, nan, nan])
#Coordinates:
#x (x) float64 0.5 1.5 2.5
#Attributes: (0)
而cubic
方法(包括quadratic
)则会返回全部为缺测值的数组。
为了避免这种情况,我们可以选择在插值前去除缺测值: dropna()
。
dropped = da.dropna('x')
print(dropped)
print(dropped.interp(x=[0.5, 1.5, 2.5], method = "cubic"))
#<xarray.DataArray (x: 4)>
#array([0. , 2. , 3. , 3.25])
#Coordinates:
# * x (x) int64 0 1 3 4
#<xarray.DataArray (x: 3)>
#array([1.19010417, 2.5078125 , 2.9296875 ])
#Coordinates:
# * x (x) float64 0.5 1.5 2.5
但如果去除缺测值,会发现原来的数组形状发生改变;若是一个多维数组,那么采取这种方法会损失大量的信息。
在这里,Xarray提供了填补缺测值的方法interpolate_na()
,利用插值的方法将缺测的位置填补,该方法与pandas.Series.interpolate()
相似。
filled = da.interpolate_na(dim="x")
print(filled)
#<xarray.DataArray (x: 5)>
#array([0. , 2. , 2.5 , 3. , 3.25])
#Coordinates:
# * x (x) int64 0 1 2 3 4
在填补缺测值后即可进行插值:
print(filled.interp(x=[0.5, 1.5, 2.5], method = "cubic"))
#<xarray.DataArray (x: 3)>
#array([1.30859375, 2.31640625, 2.73828125])
#Coordinates:
# * x (x) float64 0.5 1.5 2.5
总结
Xarray关于插值的用法还是非常好用的,可以实现包括:
1、填补缺测
2、站点插格点
3、任意剖面截取
4、不同分辨率格点互插
5、规则网格插非规则网格
6、数据延长