美文网首页xarrayPython与大气科学python
Python气象数据处理进阶之Xarray(3):插值

Python气象数据处理进阶之Xarray(3):插值

作者: 摸鱼咯 | 来源:发表于2020-04-07 21:59 被阅读0次

    Xarray库关于插值的操作内容不多,但是很实用,需要安装scipy库。
    首先我们还是创建一个DataArray:

    da = xr.DataArray(np.sin(0.3 * np.arange(12).reshape(4, 3)),[('time', np.arange(4)),('lat', [30, 60, 90])])
    print(da)
    #<xarray.DataArray (time: 4, lat: 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
    #  * lat      (lat) int64 30 60 90
    

    exp1

    我们先对时间维插值,比如说这里有第0天,1天,2天和3天,我们想获取第1.5天的数据:

    print(da.interp(time=1.5))
    #<xarray.DataArray (lat: 3)>
    #array([0.87858727, 0.89762423, 0.83647908])
    #Coordinates:
    #  * lat      (lat) int64 30 60 90
    #    time     float64 1.5
    

    可以看出实际上是对第1天和第2天的数据取了平均.
    如果要往外插,比如说得到第3.5天的数据,得到的结果全是np.nan
    也可以通过datetime格式的字符串指定插值时间:da.interp(time='2000-01-02')
    也可以同时指定对多个维度插值da.interp(time=1.5,lat=45)
    通过这个功能,就能实现不同网格的插值,甚至是站点数据到格点的插值。

    exp2

    这部分主要 是介绍几种插值的方法

    import matplotlib.pyplot as plt
    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')
    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()函数通过参数设置,可以实用不同的插值方法,比如linear线性,cubic立方等等
    为了解决EXP1中提到的外插全是缺测的问题,可以通过如下参数的设置进行调整:

    # 用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.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
    # 外推
    da.interp(x=np.linspace(-0.5, 1.5, 10), kwargs={'fill_value': 'extrapolate'})
    #<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
    

    原本小于0和大于1的部分是没用数值的,但是通过这种方法,填充或外推出了对应的数值。

    exp3

    接下来实现更高端一点的操作:


    插值exp3

    左图表示通过索引,选取出了一条斜线。
    右图表示通过插值,选取出了一条不在格点上的斜线。
    这种方法,在我们绘制不规则垂直剖面时非常需要。

    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])])
    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)
    print(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
    #Dimensions without coordinates: z
    

    以上代码实现了左图的操作。

    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])])
    x = xr.DataArray([0.5, 1.5, 2.5], dims='z')
    y = xr.DataArray([0.15, 0.25, 0.35], dims='z')
    print(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
    #Dimensions without coordinates: z
    

    以上代码实现了右图的操作。
    Metpy库中有类似操作的实现,本质就是对这个函数的封装。

    exp4

    接下来部分涉及到缺测值的问题。

    da = xr.DataArray([0, 2, np.nan, 3, 3.25], dims='x',coords={'x': range(5)})
    print(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
    

    可以看到默认设置是无法处理缺测值的,可以首先去除缺测部分:

    dropped = da.dropna('x')
    print(dropped)
    #<xarray.DataArray (x: 4)>
    #array([0.  , 2.  , 3.  , 3.25])
    #Coordinates:
    #  * x        (x) int64 0 1 3 4
    print(dropped.interp(x=[0.5, 1.5, 2.5], method='cubic'))
    #<xarray.DataArray (x: 3)>
    #array([1.19010417, 2.5078125 , 2.9296875 ])
    #Coordinates:
    #  * x        (x) float64 0.5 1.5 2.5
    

    但是这样之后数组的shape发生了变化十分不方便。
    XARRAY提供了填补缺测数据的方法:

    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
    

    在填补之后,在进行插值便可。

    总结

    Xarray关于插值的用法还是非常好用的,可以实现包括:
    1、填补缺测
    2、站点插格点
    3、任意剖面截取
    4、不同分辨率格点互插
    5、规则网格插非规则网格
    6、数据延长
    等等操作,这边的例子仅是参考官方文档的一些基础举例。

    相关文章

      网友评论

        本文标题:Python气象数据处理进阶之Xarray(3):插值

        本文链接:https://www.haomeiwen.com/subject/qhyfphtx.html