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

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、数据延长
等等操作,这边的例子仅是参考官方文档的一些基础举例。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,711评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,932评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,770评论 0 330
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,799评论 1 271
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,697评论 5 359
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,069评论 1 276
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,535评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,200评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,353评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,290评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,331评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,020评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,610评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,694评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,927评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,330评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,904评论 2 341

推荐阅读更多精彩内容