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

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

0 引言

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

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

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

1 資料來源和方法

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

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

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

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

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

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

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

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

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

2.2.2 繪制方法

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

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

水平風(fēng)速有明確的風(fēng)力等級劃分標(biāo)準(zhǔn)，因此在圖例設(shè)置時按照風(fēng)力等級進行分檔和顯示。但垂直氣流速度沒有標(biāo)準(zhǔn)可以參照劃分，并且垂直氣流運動上升和下沉?xí)r常交替出現(xiàn)，研究人員更關(guān)注垂直氣流是上升還是下沉，以及其發(fā)展趨勢。另外垂直氣流速度一般情況下較小，經(jīng)常為小于1.0 m/s，因此根據(jù)實際工作需要，將垂直氣流上升或下沉分開顯示，并分別具有多個分檔。另外為了便于分析，程序還設(shè)計了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)用個例

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

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

從水平風(fēng)16方位風(fēng)速風(fēng)向玫瑰圖統(tǒng)計結(jié)果可得出：7月上旬六盤山東坡山腰位置水平風(fēng)主要風(fēng)向為南風(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時，風(fēng)頻占比為2.9%；風(fēng)速在1.6～3.4 m/s時，風(fēng)頻占比為7.1%；風(fēng)速在3.4～5.5 m/s時，風(fēng)頻占比為4.0%。從8方位風(fēng)速風(fēng)向玫瑰圖中可以看出水平風(fēng)主要風(fēng)向為南風(fēng)、東南、西南，風(fēng)頻占比分別為24.0%，13.0%，16.0%左右。其中南風(fēng)風(fēng)向中風(fēng)速在0.3～1.6 m/s時，風(fēng)頻占比為6.0%；風(fēng)速在1.6～3.4 m/s時，風(fēng)頻占比為12.0%；風(fēng)速在3.4～5.5 m/s時，風(fēng)頻占比為6.0%。圖3為16方位水平風(fēng)速風(fēng)向玫瑰圖。

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

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

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

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

4 結(jié)束語

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