Python氣象數據處理進階之Xarray(3):插值

Xarray庫關于插值的操作內容不多,但是很實用叮叹,需要安裝scipy庫鸯乃。
首先我們還是創(chuàng)建一個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

左圖表示通過索引,選取出了一條斜線涨享。
右圖表示通過插值筋搏,選取出了一條不在格點上的斜線。
這種方法灰伟,在我們繪制不規(guī)則垂直剖面時非常需要拆又。

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發(fā)生了變化十分不方便。
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兑障、規(guī)則網格插非規(guī)則網格
6侄非、數據延長
等等操作,這邊的例子僅是參考官方文檔的一些基礎舉例流译。

?著作權歸作者所有,轉載或內容合作請聯系作者
  • 序言:七十年代末逞怨,一起剝皮案震驚了整個濱河市,隨后出現的幾起案子福澡,更是在濱河造成了極大的恐慌叠赦,老刑警劉巖,帶你破解...
    沈念sama閱讀 216,470評論 6 501
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件革砸,死亡現場離奇詭異除秀,居然都是意外死亡,警方通過查閱死者的電腦和手機业岁,發(fā)現死者居然都...
    沈念sama閱讀 92,393評論 3 392
  • 文/潘曉璐 我一進店門鳞仙,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人笔时,你說我怎么就攤上這事≌贪叮” “怎么了允耿?”我有些...
    開封第一講書人閱讀 162,577評論 0 353
  • 文/不壞的土叔 我叫張陵,是天一觀的道長扒怖。 經常有香客問我较锡,道長,這世上最難降的妖魔是什么盗痒? 我笑而不...
    開封第一講書人閱讀 58,176評論 1 292
  • 正文 為了忘掉前任蚂蕴,我火速辦了婚禮,結果婚禮上俯邓,老公的妹妹穿的比我還像新娘骡楼。我一直安慰自己,他們只是感情好稽鞭,可當我...
    茶點故事閱讀 67,189評論 6 388
  • 文/花漫 我一把揭開白布鸟整。 她就那樣靜靜地躺著,像睡著了一般朦蕴。 火紅的嫁衣襯著肌膚如雪篮条。 梳的紋絲不亂的頭發(fā)上弟头,一...
    開封第一講書人閱讀 51,155評論 1 299
  • 那天,我揣著相機與錄音涉茧,去河邊找鬼赴恨。 笑死,一個胖子當著我的面吹牛伴栓,可吹牛的內容都是我干的嘱支。 我是一名探鬼主播,決...
    沈念sama閱讀 40,041評論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼挣饥,長吁一口氣:“原來是場噩夢啊……” “哼除师!你這毒婦竟也來了?” 一聲冷哼從身側響起扔枫,我...
    開封第一講書人閱讀 38,903評論 0 274
  • 序言:老撾萬榮一對情侶失蹤汛聚,失蹤者是張志新(化名)和其女友劉穎,沒想到半個月后短荐,有當地人在樹林里發(fā)現了一具尸體倚舀,經...
    沈念sama閱讀 45,319評論 1 310
  • 正文 獨居荒郊野嶺守林人離奇死亡,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內容為張勛視角 年9月15日...
    茶點故事閱讀 37,539評論 2 332
  • 正文 我和宋清朗相戀三年忍宋,在試婚紗的時候發(fā)現自己被綠了痕貌。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點故事閱讀 39,703評論 1 348
  • 序言:一個原本活蹦亂跳的男人離奇死亡糠排,死狀恐怖舵稠,靈堂內的尸體忽然破棺而出,到底是詐尸還是另有隱情入宦,我是刑警寧澤哺徊,帶...
    沈念sama閱讀 35,417評論 5 343
  • 正文 年R本政府宣布,位于F島的核電站乾闰,受9級特大地震影響落追,放射性物質發(fā)生泄漏。R本人自食惡果不足惜涯肩,卻給世界環(huán)境...
    茶點故事閱讀 41,013評論 3 325
  • 文/蒙蒙 一轿钠、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧病苗,春花似錦疗垛、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,664評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至,卻和暖如春花履,著一層夾襖步出監(jiān)牢的瞬間芽世,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 32,818評論 1 269
  • 我被黑心中介騙來泰國打工诡壁, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留济瓢,地道東北人。 一個月前我還...
    沈念sama閱讀 47,711評論 2 368
  • 正文 我出身青樓妹卿,卻偏偏與公主長得像旺矾,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子夺克,可洞房花燭夜當晚...
    茶點故事閱讀 44,601評論 2 353

推薦閱讀更多精彩內容