xarray （教程）第四章插值的使用

西大·阿志

已于 2024-01-15 10:10:18 修改

阅读量1.6k

点赞数 15

分类专栏： xarray使用教程文章标签： python

于 2024-01-15 10:07:55 首次发布

本文链接：https://blog.csdn.net/qq_46608464/article/details/135594634

版权

xarray使用教程专栏收录该内容

5 篇文章 2 订阅

订阅专栏

Interpolating data

Xarray提供了灵活的插值例程，与我们的索引有相似的接口。
interp需要安装scipy。

Scalar and 1-dimensional interpolation

对数据数组进行插值的工作方式很像对数据数组进行标记索引，

da = xr.DataArray(
    np.sin(0.3 * np.arange(12).reshape(4, 3)),
    [("time", np.arange(4)), ("space", [0.1, 0.2, 0.3])],
)

# label lookup
da.sel(time=3)
"""Out[2]: 
<xarray.DataArray (space: 3)>
array([ 0.427,  0.141, -0.158])
Coordinates:
    time     int64 3
  * space    (space) float64 0.1 0.2 0.3"""

# interpolation
da.interp(time=2.5)
"""Out[3]: 
<xarray.DataArray (space: 3)>
array([0.701, 0.502, 0.259])
Coordinates:
  * space    (space) float64 0.1 0.2 0.3
    time     float64 2.5"""

与索引类似，interp()也接受类似数组的，它将插值结果作为数组给出。

# label lookup
da.sel(time=[2, 3])
"""Out[4]: 
<xarray.DataArray (time: 2, space: 3)>
array([[ 0.974,  0.863,  0.675],
       [ 0.427,  0.141, -0.158]])
Coordinates:
  * time     (time) int64 2 3
  * space    (space) float64 0.1 0.2 0.3
"""
# interpolation
da.interp(time=[2.5, 3.5])
"""Out[5]: 
<xarray.DataArray (time: 2, space: 3)>
array([[0.701, 0.502, 0.259],
       [  nan,   nan,   nan]])
Coordinates:
  * space    (space) float64 0.1 0.2 0.3
  * time     (time) float64 2.5 """
3.5

要使用numpy.datetime64坐标插值数据，可以传递一个字符串。

da_dt64 = xr.DataArray(
    [1, 3], [("time", pd.date_range("1/1/2000", "1/3/2000", periods=2))]
)

da_dt64.interp(time="2000-01-02")
"""Out[7]: 
<xarray.DataArray ()>
array(2.)
Coordinates:
    time     datetime64[ns] 2000-01-02"""

通过指定所需的时间段，可以将插值数据合并到原始数据阵列中。

da_dt64.interp(time=pd.date_range("1/1/2000", "1/3/2000", periods=3))
"""Out[8]: 
<xarray.DataArray (time: 3)>
array([1., 2., 3.])
Coordinates:
  * time     (time) datetime64[ns] 2000-01-01 2000-01-02 2000-01-03"""

也允许对CFTimeIndex索引的数据进行插值。
目前，我们的插值只适用于规则的网格。因此，与sel()类似，只有沿一个维度的1D坐标可以用作要插值的原始坐标。

Multi-dimensional Interpolation

和sel()一样，interp()接受多个坐标。在这种情况下，执行多维插值。

# label lookup
da.sel(time=2, space=0.1)
"""Out[9]: 
<xarray.DataArray ()>
array(0.974)
Coordinates:
    time     int64 2
    space    float64 0.1"""

# interpolation
da.interp(time=2.5, space=0.15)
"""Out[10]: 
<xarray.DataArray ()>
array(0.601)
Coordinates:
    time     float64 2.5
    space    float64 0.15"""

也接受类似数组的坐标:

# label lookup
da.sel(time=[2, 3], space=[0.1, 0.2])
"""Out[11]: 
<xarray.DataArray (time: 2, space: 2)>
array([[0.974, 0.863],
       [0.427, 0.141]])
Coordinates:
  * time     (time) int64 2 3
  * space    (space) float64 0.1 0.2"""

# interpolation
da.interp(time=[1.5, 2.5], space=[0.15, 0.25])
"""Out[12]: 
<xarray.DataArray (time: 2, space: 2)>
array([[0.888, 0.867],
       [0.601, 0.381]])
Coordinates:
  * time     (time) float64 1.5 2.5
  * space    (space) float64 0.15 0.25"""

interp_like()方法是一个有用的快捷方式。该方法将一个xarray对象插值到另一个xarray对象的坐标上。
例如，如果我们想计算两个数据数组(da和other)在稍微不同的坐标上的差异，
首先插值da，使其保持在相同的坐标上，然后减去它。interp_like()可以用于这种情况，

da = xr.DataArray(
    np.arange(12).reshape(4, 3),
    [("time", np.arange(4)), ("space", [0.1, 0.2, 0.3])],
)
other = xr.DataArray(
    np.arange(9).reshape(3, 3),
    [("time", [0.9, 1.9, 2.9]), ("space", [0.15, 0.25, 0.35])],
)
print(da)
print(other)
# interpolate da along other's coordinates
interpolated = da.interp_like(other)
print(interpolated)
"""
<xarray.DataArray (time: 4, space: 3)>
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11]])
Coordinates:
  * time     (time) int32 0 1 2 3
  * space    (space) float64 0.1 0.2 0.3
<xarray.DataArray (time: 3, space: 3)>
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
Coordinates:
  * time     (time) float64 0.9 1.9 2.9
  * space    (space) float64 0.15 0.25 0.35
<xarray.DataArray (time: 3, space: 3)>
array([[ 3.2,  4.2,  nan],
       [ 6.2,  7.2,  nan],
       [ 9.2, 10.2,  nan]])
Coordinates:
  * time     (time) float64 0.9 1.9 2.9
  * space    (space) float64 0.15 0.25 0.35
"""

Interpolation methods

我们使用scipy.interpolate.interp1d进行一维插值。对于多维插值，首先尝试将插值分解为一系列一维插值，在这种情况下，使用scipy.interpolate.interp1d。如果无法进行分解(
例如使用高级插值)，则使用scipy.interpolate.interpn()。
插值方法可以由可选的方法参数指定。

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="original")
# Out[17]: [<matplotlib.lines.Line2D at 0x7f6e93697ee0>]

da.interp(x=np.linspace(0, 1, 100)).plot.line(label="linear (default)")
# Out[18]: [<matplotlib.lines.Line2D at 0x7f6e936969b0>]

da.interp(x=np.linspace(0, 1, 100), method="cubic").plot.line(label="cubic")
# Out[19]: [<matplotlib.lines.Line2D at 0x7f6e93696b90>]

plt.legend()
# Out[20]: <matplotlib.legend.Legend at 0x7f6e83d6d810>

在这里插入图片描述

附加的关键字参数可以传递给scipy的函数。

# fill 0 for the outside of the original coordinates.
da.interp(x=np.linspace(-0.5, 1.5, 10), kwargs={"fill_value": 0.0})
"""Out[21]: 
<xarray.DataArray (x: 10)>
array([ 0.   ,  0.   ,  0.   ,  0.814,  0.604, -0.604, -0.814,  0.   ,  0.   ,  0.   ])
Coordinates:
  * x        (x) float64 -0.5 -0.2778 -0.05556 0.1667 ... 0.8333 1.056 1.278 1.5"""

# 1-dimensional extrapolation
da.interp(x=np.linspace(-0.5, 1.5, 10), kwargs={"fill_value": "extrapolate"})
"""Out[22]: 
<xarray.DataArray (x: 10)>
array([-2.893, -1.607, -0.321,  0.814,  0.604, -0.604, -0.814,  0.321,  1.607,  2.893])
Coordinates:
  * x        (x) float64 -0.5 -0.2778 -0.05556 0.1667 ... 0.8333 1.056 1.278 1.5"""

# multi-dimensional extrapolation
da = xr.DataArray(
    np.sin(0.3 * np.arange(12).reshape(4, 3)),
    [("time", np.arange(4)), ("space", [0.1, 0.2, 0.3])],
)


da.interp(
    time=4, space=np.linspace(-0.1, 0.5, 10), kwargs={"fill_value": "extrapolate"}
)

"""Out[24]: 
<xarray.DataArray (space: 10)>
array([ 0.805,  0.497,  0.189, -0.119, -0.427, -0.718, -0.991, -1.264, -1.538, -1.811])
Coordinates:
    time     int64 4
  * space    (space) float64 -0.1 -0.03333 0.03333 0.1 ... 0.3 0.3667 0.4333 0.5"""

Advanced Interpolation

interp()接受类似于sel()的DataArray，这使我们能够进行更高级的插值。基于传递给interp()的新坐标的尺寸，结果的尺寸被确定。例如，如果您想沿着特定的维度插入一个二维数组，如下所示，您可以传递两个具有相同维度的一维数据数组作为新坐标。

在这里插入图片描述

da = xr.DataArray(
    np.arange(20).reshape(5, 4),
    [("x", np.arange(5)), ("y", [1,2,3,4])],
)

print(da)
# advanced indexing
x = xr.DataArray([0, 2, 4], dims="z")

y = xr.DataArray([1,2,3], dims="z")

da=da.sel(x=x, y=y)
print(da)
"""
<xarray.DataArray (x: 5, y: 4)>
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15],
       [16, 17, 18, 19]])
Coordinates:
  * x        (x) int32 0 1 2 3 4
  * y        (y) int32 1 2 3 4
<xarray.DataArray (z: 3)>
array([ 0,  9, 18])
Coordinates:
    x        (z) int32 0 2 4
    y        (z) int32 1 2 3
Dimensions without coordinates: z
"""
da = xr.DataArray(
    np.arange(20).reshape(5, 4),
    [("x", np.arange(5)), ("y", [1,2,3,4])],
)

print(da)
# advanced interpolation, without extrapolation
x = xr.DataArray([0.5, 1.5, 2.5, 3.5], dims="z")

y = xr.DataArray([1.5, 2.5, 3.5, 4.5], dims="z")

da=da.interp(x=x, y=y)
print(da)
"""
<xarray.DataArray (x: 5, y: 4)>
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15],
       [16, 17, 18, 19]])
Coordinates:
  * x        (x) int32 0 1 2 3 4
  * y        (y) int32 1 2 3 4
<xarray.DataArray (z: 4)>
array([ 2.5,  7.5, 12.5,  nan])
Coordinates:
    x        (z) float64 0.5 1.5 2.5 3.5
    y        (z) float64 1.5 2.5 3.5 4.5
Dimensions without coordinates: z
"""

其中原始坐标(x，y) = ((0.5，1.5)，(1.5，2.5)，(2.5，3.5)，(3.5，4.5))上的值通过二维插值获得，并沿着新的维度z映射。
由于没有关键字参数被传递到插值例程，因此不执行外插，从而产生nan值。
如果要向新维度z添加坐标，可以向DataArray s提供坐标。外推可以通过向SciPy的interpnd函数传递附加参数来实现，

da = xr.DataArray(
    np.arange(20).reshape(5, 4),
    [("x", np.arange(5)), ("y", [1,2,3,4])],
)

print(da)
x = xr.DataArray([0.5, 1.5, 2.5, 3.5], dims="z", coords={"z": ["a", "b", "c", "d"]})

y = xr.DataArray(
    [1.5, 2.5, 3.5, 4.5], dims="z", coords={"z": ["a", "b", "c", "d"]}
)


da=da.interp(x=x, y=y, kwargs={"fill_value": None})
print(da)
"""
<xarray.DataArray (x: 5, y: 4)>
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15],
       [16, 17, 18, 19]])
Coordinates:
  * x        (x) int32 0 1 2 3 4
  * y        (y) int32 1 2 3 4
<xarray.DataArray (z: 4)>
array([ 2.5,  7.5, 12.5, 17.5])
Coordinates:
    x        (z) float64 0.5 1.5 2.5 3.5
    y        (z) float64 1.5 2.5 3.5 4.5
  * z        (z) <U1 'a' 'b' 'c' 'd'
"""

Interpolating arrays with NaN

我们的interp()处理NaN数组的方式与scipy.interpolate.interp1d和scipy.interpolate.interpn相同。线性和最近方法返回包括NaN的数组，而其他方法如三次或二次方法返回所有NaN数组。

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])
"""Out[36]: 
<xarray.DataArray (x: 3)>
array([ 1., nan, nan])
Coordinates:
  * x        (x) float64 0.5 1.5 2.5"""

da.interp(x=[0.5, 1.5, 2.5], method="cubic")
"""Out[37]: 
<xarray.DataArray (x: 3)>
array([nan, nan, nan])
Coordinates:
  * x        (x) float64 0.5 1.5 2.5"""

为了避免这种情况，您可以通过dropna()删除NaN，然后进行插值

dropped = da.dropna("x")

dropped
"""Out[39]: 
<xarray.DataArray (x: 4)>
array([0.  , 2.  , 3.  , 3.25])
Coordinates:
  * x        (x) int64 0 1 3 4"""

dropped.interp(x=[0.5, 1.5, 2.5], method="cubic")
"""Out[40]: 
<xarray.DataArray (x: 3)>
array([1.19 , 2.508, 2.93 ])
Coordinates:
  * x        (x) float64 0.5 1.5 2.5"""

如果nan在多维数组中随机分布，那么dropna()删除包含多个nan的所有列可能会丢失大量信息。在这种情况下，可以通过interpolate_na()填充NaN，类似于pandas.Series.interpolate()。
这将通过沿指定尺寸插值来填充NaN。填充NaNs后，您可以进行插值:

filled = da.interpolate_na(dim="x")

filled
"""Out[42]: 
<xarray.DataArray (x: 5)>
array([0.  , 2.  , 2.5 , 3.  , 3.25])
Coordinates:
  * x        (x) int64 0 1 2 3 4"""
filled.interp(x=[0.5, 1.5, 2.5], method="cubic")
"""Out[43]: 
<xarray.DataArray (x: 3)>
array([1.309, 2.316, 2.738])
Coordinates:
  * x        (x) float64 0.5 1.5 2.5"""