参考：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. 高端操作

advanced_selection_interpolation.png

上图中，图左是通过索引，得到一条红色的斜线，图右表示通过插值得到一条红色的斜线。
以下代码实现左图操作：

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、数据延长