基于Python的三維風(fēng)場(chǎng)風(fēng)玫瑰圖繪制方法研究與應(yīng)用

周積強(qiáng)，魏 巍，姚肖萌，王雪妮，李龍艷，楊 勇，柳佳俊

(1.寧夏回族自治區(qū)氣象災(zāi)害防御技術(shù)中心，銀川 750002；2.銀川市氣象局，銀川 750002；3.寧夏氣象防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室，銀川 750002)

0 引言

空氣的水平運(yùn)動(dòng)稱(chēng)為風(fēng)，不僅有數(shù)值的大小，還具有方向。地面風(fēng)向一般用16方位表示，每相鄰方位間的角差為22.5°。風(fēng)玫瑰圖是統(tǒng)計(jì)和顯示風(fēng)的專(zhuān)業(yè)方法，是指極坐標(biāo)底圖上繪制出某一地區(qū)在某一時(shí)間段內(nèi)各風(fēng)向出現(xiàn)的頻率或各風(fēng)向的平均風(fēng)速的統(tǒng)計(jì)圖，相應(yīng)分為風(fēng)向玫瑰圖和風(fēng)速玫瑰圖。在真實(shí)大氣中風(fēng)還包括垂直運(yùn)動(dòng)，即上升或下沉氣流。三維風(fēng)速儀能夠?qū)崟r(shí)觀(guān)測(cè)水平、垂直風(fēng)速及風(fēng)向，廣泛應(yīng)用于測(cè)量風(fēng)場(chǎng)、湍流等觀(guān)測(cè)研究中。

風(fēng)速數(shù)據(jù)處理和玫瑰圖的繪制有很多方式。常見(jiàn)的有使用專(zhuān)用的風(fēng)統(tǒng)計(jì)軟件、Excel[1]、VB[2]、VC[3]、C#[4]、Extjs[5]、Python[6]等進(jìn)行繪制。Python具有強(qiáng)大的數(shù)據(jù)處理和分析、可視化、豐富的標(biāo)準(zhǔn)庫(kù)及第三方庫(kù)的特點(diǎn)，目前已在氣象領(lǐng)域得到廣泛應(yīng)用[7]。但使用Python對(duì)大氣三維風(fēng)繪制玫瑰圖仍然較少，包括水平風(fēng)和垂直氣流風(fēng)速風(fēng)向玫瑰圖。

根據(jù)風(fēng)速風(fēng)向玫瑰圖可視化展示特點(diǎn)，文章分析了三維風(fēng)數(shù)據(jù)和數(shù)據(jù)結(jié)構(gòu)存在的問(wèn)題，使用三倍標(biāo)準(zhǔn)差等方法對(duì)數(shù)據(jù)進(jìn)行了質(zhì)量控制。基于Python平臺(tái)使用第三方庫(kù)Windrose庫(kù)，設(shè)計(jì)了一種大氣三維風(fēng)包括水平風(fēng)和垂直氣流風(fēng)速風(fēng)向玫瑰圖的繪制方法。使用實(shí)測(cè)三維風(fēng)資料進(jìn)行處理和應(yīng)用，結(jié)果表明繪制的水平風(fēng)和垂直氣流風(fēng)玫瑰圖能夠清晰地顯示出風(fēng)的三維風(fēng)統(tǒng)計(jì)特征，為近地層大氣三維風(fēng)場(chǎng)研究提供了一種統(tǒng)計(jì)方法。

1 資料來(lái)源和方法

文章使用的三維風(fēng)速儀是超聲波風(fēng)速風(fēng)向儀，能夠記錄水平風(fēng)速、垂直氣流、風(fēng)向等信息。數(shù)據(jù)保存格式為文本格式。采樣數(shù)據(jù)存在一些問(wèn)題，主要包括存儲(chǔ)數(shù)據(jù)包含一些將風(fēng)速風(fēng)向數(shù)據(jù)記錄為NAN記錄。有些數(shù)據(jù)記錄為奇異值，還有重復(fù)記錄，如同一時(shí)間會(huì)出現(xiàn)兩條以上的數(shù)據(jù)記錄。極少情況下水平風(fēng)數(shù)據(jù)出現(xiàn)負(fù)號(hào)。最后數(shù)據(jù)在使用前需進(jìn)行粗大值篩選和剔除。使用3倍標(biāo)準(zhǔn)差法將粗大值剔除和插補(bǔ)。查找重復(fù)數(shù)據(jù)記錄，僅保留首次記錄的數(shù)據(jù)。將出現(xiàn)的個(gè)別空值使用鄰近法進(jìn)行插補(bǔ)，從而得到經(jīng)數(shù)據(jù)清洗后的數(shù)據(jù)集。

2 數(shù)據(jù)分析流程和關(guān)鍵技術(shù)

文章使用Python第三方庫(kù)Windrose來(lái)繪制三維風(fēng)(水平風(fēng)和垂直風(fēng))風(fēng)速風(fēng)向玫瑰圖。僅安裝Windrose第三方庫(kù)并不能順利使用其函數(shù)進(jìn)行風(fēng)玫瑰圖繪制，同時(shí)還需要安裝相關(guān)的依賴(lài)庫(kù)，如Pandas，Numpy，Matplotlib，Scipy等。

2.1 數(shù)據(jù)分析流程

三維風(fēng)數(shù)據(jù)存儲(chǔ)格式為文本格式，在多個(gè)文檔中存儲(chǔ)。處理之前需選擇待分析的多個(gè)文檔合并、數(shù)據(jù)清洗。然后將需要繪制的風(fēng)向、水平風(fēng)、垂直氣流信息提取出，使用Windrose第三方庫(kù)分別完成水平風(fēng)和垂直氣流速度的風(fēng)速風(fēng)向玫瑰圖。數(shù)據(jù)分析流程如圖1所示。

2.2 關(guān)鍵技術(shù)

2.2.1 數(shù)據(jù)預(yù)處理

使用os和tkinter庫(kù)函數(shù)選擇需要處理的數(shù)據(jù)文件。讀取數(shù)據(jù)使用函數(shù)pd.read_csv()。將循環(huán)讀取的多個(gè)DataFrame表分別追加暫存入列表，使用concat()函數(shù)將列表內(nèi)的DataFrame表進(jìn)行垂直合并。按照數(shù)據(jù)質(zhì)量控制方法使用Pandas庫(kù)中的replace()函數(shù)將數(shù)據(jù)表中的NAN信息替換為程序能夠識(shí)別的np.nan。使用astypes()函數(shù)強(qiáng)制轉(zhuǎn)換數(shù)據(jù)類(lèi)型。待數(shù)據(jù)空值和類(lèi)型全部確定后，使用3倍標(biāo)準(zhǔn)差方式分別篩選出水平風(fēng)和垂直氣流速度的粗大值并將其定義為空值，將出現(xiàn)水平風(fēng)為負(fù)值的數(shù)據(jù)定義為空值。最后將篩選出的空值使用就近填充法對(duì)其進(jìn)行填充。

圖1 數(shù)據(jù)處理流程

2.2.2 繪制方法

直接使用Windrose中WindroseAxes()函數(shù)繪制風(fēng)速風(fēng)向玫瑰圖時(shí)，其顯示內(nèi)容不符合使用習(xí)慣和需求，要將風(fēng)向信息由度數(shù)調(diào)整為風(fēng)向表示方式。繪制出的0°方向指向右側(cè)，需要將其旋轉(zhuǎn)設(shè)置為正上方為N方向或0°方向。數(shù)據(jù)旋轉(zhuǎn)方式為逆時(shí)針?lè)较?，需將風(fēng)向旋轉(zhuǎn)方向調(diào)整為順時(shí)針?lè)较?。最重要的是直接使用顯示數(shù)據(jù)結(jié)果存在錯(cuò)誤，如將0°方向繪制在90°方向，將90°方向繪制在0°方向。使用ax.set_theta_zero_location('N') 函數(shù)調(diào)整方位，使用ax.set_theta_direction('clockwise')將數(shù)據(jù)旋轉(zhuǎn)方式調(diào)整為順時(shí)針?lè)绞斤@示。對(duì)風(fēng)向數(shù)據(jù)進(jìn)行轉(zhuǎn)換處理，經(jīng)調(diào)整后數(shù)據(jù)顯示為正確結(jié)果，如圖2所示。

圖2 調(diào)整后風(fēng)速風(fēng)向玫瑰圖

水平風(fēng)速有明確的風(fēng)力等級(jí)劃分標(biāo)準(zhǔn)，因此在圖例設(shè)置時(shí)按照風(fēng)力等級(jí)進(jìn)行分檔和顯示。但垂直氣流速度沒(méi)有標(biāo)準(zhǔn)可以參照劃分，并且垂直氣流運(yùn)動(dòng)上升和下沉?xí)r常交替出現(xiàn)，研究人員更關(guān)注垂直氣流是上升還是下沉，以及其發(fā)展趨勢(shì)。另外垂直氣流速度一般情況下較小，經(jīng)常為小于1.0 m/s，因此根據(jù)實(shí)際工作需要，將垂直氣流上升或下沉分開(kāi)顯示，并分別具有多個(gè)分檔。另外為了便于分析，程序還設(shè)計(jì)了8方位、16方位和帶有中文或英文兩種不同顯示風(fēng)向方式的玫瑰圖。繪制方法程序詳細(xì)如下：

風(fēng)向數(shù)據(jù)轉(zhuǎn)換程序：

defwind_direction_change(wd)：

wind_direction = np.array(wd)

wind_direction[np.where((wind_direction >= 0) & (wind_direction <= 135))[0]] =90 - wind_direction[np.where((wind_direction>=0)& (wind_direction <= 135))[0]]

wind_direction[np.where((wind_direction>= 315) & (wind_direction < 360))[0]] = 90 - wind_direction[np.where((wind_direction >= 315) & (wind_direction < 360))[0]]

wind_direction[np.where((wind_direction > 135) & (wind_direction< 315))[0]] = 450 - wind_direction[np.where((wind_direction>135)& (wind_direction < 315))[0]]

wind_direction[np.where(wind_direction < 0)[0]] = 360 + wind_direction[np.where(wind_direction < 0)[0]]

returnwind_direction

風(fēng)速風(fēng)向玫瑰圖繪制程序：

defwind_rose_plt(wd，ws，levl=16，bin_num=6，seris=0)：

wd_new = wind_direction_change(wd)

ax =WindroseAxes.from_ax()

ax.bar(direction=wd_new，var=ws，nsector=levl，bins=bin_num，normed=True，cmap=cm.tab10，opening=0.8，edgecolor='white')

ax.set_theta_zero_location('N')

ax.set_theta_direction('clockwise')

ax.set_legend(loc='lower center'，title='風(fēng)速(單位：m/s)'，bbox_to_anchor=[0，0]，ncol=1)

ax.set_xticks([(i / (levl / 2)) * np.pi for i in range(levl)])

iflevl == 16 and seris == 0:

ax.set_xticklabels(['N'，'NNE'，'NE'，'ENE'，'E'，'ESE'，'SE'，'SSE'，'S'，'SSW'，'SW'，'WSW'，'W'，'WNW'，'NW'，'NNW']，fontdict={'family'：'SimHei'，'size'：15，'color'：'black'})

elif levl == 16 and seris == 1:

ax.set_xticklabels(['北'，'東北偏北'，'東北'，'東北偏東'，'東'，'東南偏東'，'東南'，'東南偏南'，'南'，'西南偏南'，'西南'，'西南偏西'，'西'，'西北偏西'，'西北'，'西北偏北']，fontdict={'family'：'SimHei'，'size'：15，'color'：'black'})

elif levl == 8 and seris == 0:

ax.set_xticklabels(['N'，'NE'，'E'，'SE'，'S'，'SW'，'W'，'NW']，fontdict={'family'：'SimHei'，'size'：15，'color'：'red'})

elif levl == 8 and seris == 1:

ax.set_xticklabels(['北'，'東北'，'東'，'東南'，'南'，'西南'，'西'，'西北']，fontdict={'family'：'SimHei'，'size'：15，'color'：'red'})

3 應(yīng)用個(gè)例

文章使用數(shù)據(jù)的時(shí)間段為2021-07-01—2021-07-10。按照設(shè)計(jì)的程序，對(duì)個(gè)例數(shù)據(jù)進(jìn)行測(cè)試應(yīng)用。

3.1 水平風(fēng)速風(fēng)向玫瑰圖

從水平風(fēng)16方位風(fēng)速風(fēng)向玫瑰圖統(tǒng)計(jì)結(jié)果可得出：7月上旬六盤(pán)山東坡山腰位置水平風(fēng)主要風(fēng)向?yàn)槟巷L(fēng)、西南偏南、西南風(fēng)、東南偏南，風(fēng)頻占比分別為14.0%，9.4%，8.8%，9.0%左右。其中南風(fēng)風(fēng)向中風(fēng)速在0.3～1.6 m/s時(shí)，風(fēng)頻占比為2.9%；風(fēng)速在1.6～3.4 m/s時(shí)，風(fēng)頻占比為7.1%；風(fēng)速在3.4～5.5 m/s時(shí)，風(fēng)頻占比為4.0%。從8方位風(fēng)速風(fēng)向玫瑰圖中可以看出水平風(fēng)主要風(fēng)向?yàn)槟巷L(fēng)、東南、西南，風(fēng)頻占比分別為24.0%，13.0%，16.0%左右。其中南風(fēng)風(fēng)向中風(fēng)速在0.3～1.6 m/s時(shí)，風(fēng)頻占比為6.0%；風(fēng)速在1.6～3.4 m/s時(shí)，風(fēng)頻占比為12.0%；風(fēng)速在3.4～5.5 m/s時(shí)，風(fēng)頻占比為6.0%。圖3為16方位水平風(fēng)速風(fēng)向玫瑰圖。

3.2 垂直氣流風(fēng)玫瑰圖

從垂直氣流16方位風(fēng)速風(fēng)向玫瑰圖統(tǒng)計(jì)結(jié)果可得出：7月上旬在六盤(pán)山東坡山腰位置，當(dāng)風(fēng)向?yàn)闁|南風(fēng)時(shí)，垂直氣流表現(xiàn)為上升運(yùn)動(dòng)，當(dāng)風(fēng)向?yàn)槲髂?、西風(fēng)、西北時(shí)，垂直氣流表現(xiàn)為下沉和上升運(yùn)動(dòng)均存在。其中當(dāng)西風(fēng)時(shí)，下沉氣流占比大于下沉氣流。其中當(dāng)風(fēng)向?yàn)槟巷L(fēng)時(shí)，下沉氣流風(fēng)速在-2.7～0.0 m/s時(shí)，風(fēng)頻占比為5.8%；上升氣流風(fēng)速區(qū)間在0.0～1.1 m/s時(shí)，風(fēng)頻占比為8.2%；上升氣流風(fēng)速在1.1～2.1 m/s時(shí)，風(fēng)頻占比為0.3%。當(dāng)風(fēng)向?yàn)闁|南風(fēng)時(shí)，垂直氣流大部分為上升運(yùn)動(dòng)，下沉氣流風(fēng)速在-2.7～0.0 m/s時(shí)，風(fēng)頻占比為1.0%；上升氣流風(fēng)速在0.0～1.1 m/s，風(fēng)頻占比為4.7%；上升氣流風(fēng)速在1.1～2.1 m/s，風(fēng)頻占比為6.2%。從8方位風(fēng)速風(fēng)向玫瑰圖中可更加明顯地看出大氣垂直運(yùn)動(dòng)上升氣流與風(fēng)向的關(guān)系。當(dāng)風(fēng)向?yàn)槟巷L(fēng)時(shí)，垂直氣流下沉運(yùn)動(dòng)風(fēng)頻占比為10.1%，上升運(yùn)動(dòng)風(fēng)頻占比為14.0%。當(dāng)風(fēng)向?yàn)闁|南風(fēng)時(shí)，垂直氣流下沉運(yùn)動(dòng)風(fēng)頻占比為2.0%，上升運(yùn)動(dòng)風(fēng)頻占比為12.0%。當(dāng)風(fēng)向?yàn)闁|風(fēng)時(shí)，垂直氣流下沉運(yùn)動(dòng)風(fēng)頻占比小于1.0%，上升運(yùn)動(dòng)風(fēng)頻占比為9.0%。因六盤(pán)山山脈總體呈南北走向，當(dāng)近地層風(fēng)從西向東吹時(shí)，近地層大氣翻越山體后在山體東側(cè)下沉，觀(guān)測(cè)得到的垂直氣流表現(xiàn)為下沉運(yùn)動(dòng)。當(dāng)近地層風(fēng)從東向西吹時(shí)，近地層大氣由于山體強(qiáng)迫抬升，觀(guān)測(cè)得到的垂直氣流表現(xiàn)為上升運(yùn)動(dòng)。16方位垂直氣流風(fēng)速風(fēng)向玫瑰圖如圖4所示。

圖3 16方位水平風(fēng)速風(fēng)向玫瑰圖

圖4 16方位垂直氣流風(fēng)速風(fēng)向玫瑰圖

4 結(jié)束語(yǔ)

大氣除了水平運(yùn)動(dòng)以外還有大氣的垂直運(yùn)動(dòng)，三維風(fēng)速儀能夠觀(guān)測(cè)大氣三維運(yùn)動(dòng)情況。繪制三維風(fēng)玫瑰圖對(duì)分析研究大氣三維運(yùn)動(dòng)具有重要作用。文章分析了三維風(fēng)速儀資料特點(diǎn)，基于Python平臺(tái)使用第三方Windrose函數(shù)設(shè)計(jì)了三維風(fēng)玫瑰圖繪制程序，并使用該程序繪制寧夏六盤(pán)山地區(qū)東坡山腰實(shí)際觀(guān)測(cè)三維風(fēng)玫瑰圖，包括水平風(fēng)風(fēng)速風(fēng)向玫瑰圖和垂直氣流風(fēng)速風(fēng)向玫瑰圖。文章解決了直接使用Windrose庫(kù)函數(shù)繪制風(fēng)速風(fēng)向玫瑰圖存在風(fēng)向0°方向沒(méi)有指向北方、風(fēng)向旋轉(zhuǎn)方式為逆時(shí)針，顯示數(shù)據(jù)存在錯(cuò)誤等問(wèn)題。繪制得到水平風(fēng)和垂直氣流風(fēng)速風(fēng)向玫瑰圖能夠清晰地反應(yīng)出當(dāng)?shù)亟貙哟髿膺\(yùn)動(dòng)情況及與地形之間的關(guān)系，為相關(guān)研究和玫瑰圖自動(dòng)繪制提供了一種方法。