添加新云天氣象小編微信或QQ：130188121，及時獲取或發(fā)布?xì)庀笊龑W(xué)、就業(yè)、會議、征稿及學(xué)術(shù)動態(tài)信息！最新熱點(diǎn)文章：行業(yè)動態(tài) | 2020中國氣象現(xiàn)代化建設(shè)科技博覽會行業(yè)動態(tài) | “遠(yuǎn)征杯”中國信用防雷宣傳活動開鑼了！氣象招聘 | 中南空管局氣象中心2020年校園招聘公告

今天我們來講講風(fēng)數(shù)據(jù)的簡單處理和分析。通過這個案例，你將知道一些與風(fēng)數(shù)據(jù)處理和分析有關(guān)的Python函數(shù)庫和繪圖方法。具體內(nèi)容如下：

1)如何從網(wǎng)絡(luò)上獲取數(shù)據(jù)并讀取格式化數(shù)據(jù)(reading data)

2)如何處理時間變量(datetime variable)

3)如何判別異常值和去除異常值(Outlier)

4)如何使用Pandas.DataFrame來篩選、整合、選擇數(shù)據(jù)

5)如何繪制風(fēng)玫瑰圖

6)如何繪制數(shù)據(jù)的統(tǒng)計(jì)分布曲線

7)如何從數(shù)據(jù)中獲取Weibull分布參數(shù)

8)如何繪制16位風(fēng)向統(tǒng)計(jì)柱狀圖

9)如何繪制單點(diǎn)風(fēng)矢量動畫和GIF輸出

10)如何繪制箱式統(tǒng)計(jì)圖(boxplot)

11)風(fēng)的矢量分解函數(shù)的使用

Python編程語言已經(jīng)成為世界上最流行和受歡迎的編程語言。它的靈活性和易用性已經(jīng)強(qiáng)大的函數(shù)庫吸引了數(shù)據(jù)分析愛好者。今天我們一起來看一例使用Python語言來處理和分析風(fēng)數(shù)據(jù)的一個案例。希望這個案例能夠幫助到大家。

風(fēng)數(shù)據(jù)是我們應(yīng)用氣象專業(yè)的同行們經(jīng)常要接觸的。相比溫度、降水量、太陽輻射等氣象要素，風(fēng)要素的重要特征是它是一個矢量要素，有方向有大小。風(fēng)方向從正北方向?yàn)?°開始，順時針轉(zhuǎn)向360°，因此正東方向就對應(yīng)風(fēng)向90°，正南為180°，正西為270°，正北為0°。在氣象上通常將風(fēng)向劃分為16個方向，用符號N，NNE，NE，...來表示。風(fēng)速則根據(jù)風(fēng)速大小劃分等級，比如國際通用的風(fēng)力等級：蒲福風(fēng)級(https://baike.baidu.com/item/蒲福風(fēng)級)。因此有了風(fēng)向、風(fēng)速的分類分級，我們就可以對風(fēng)要素進(jìn)行分析，獲取不同地方的風(fēng)場特征。其中最重要的風(fēng)場特征表現(xiàn)形式是風(fēng)玫瑰圖。下面我們就風(fēng)數(shù)據(jù)的讀取、預(yù)處理和分析來做一個案例。

# 加載一些常用的Python函數(shù)庫

# numpy是支持多維數(shù)據(jù)存儲、運(yùn)算、分析函數(shù)庫

# pandas是結(jié)構(gòu)化數(shù)據(jù)的存儲、運(yùn)算、分析函數(shù)庫

# matplotlib是比用繪圖函數(shù)庫??

# seaborn是繪圖函數(shù)庫??

# warnings是運(yùn)行警告函數(shù)庫

import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as snsimport warnings

# 運(yùn)行環(huán)境設(shè)置：在Notebook中抑制警告輸出? ? ? ? ? ? ? ??

warnings.filterwarnings('ignore')

# 在notebook中顯示繪制的圖像

%matplotlib

# 加載風(fēng)數(shù)據(jù)。風(fēng)數(shù)據(jù)存放在網(wǎng)絡(luò)上。注：本案例的風(fēng)數(shù)據(jù)是在原始數(shù)據(jù)基礎(chǔ)上增加了噪聲

url?=?'https://raw.githubusercontent.com/yangsbin/Meteo_data_analysis/master/data/wind_data.txt'

# 通過Pandas的read_csv函數(shù)讀取格式化數(shù)據(jù)，參數(shù)中sep是讀取的原數(shù)據(jù)列的分隔符參數(shù)定義， header是頭行(None就是無頭行)，names定義了每個數(shù)據(jù)列的名稱；WD定義為風(fēng)向；WS定義為風(fēng)速

wind_data?=?pd.read_csv(url,?sep?=?r'\s+',?header?=?None,?names?=?['Year','Month','Day','Hour',?'WD','WS'])

# 原始數(shù)據(jù)中風(fēng)速是10倍值，因此將讀取后的數(shù)據(jù)中的風(fēng)速列數(shù)據(jù)除10轉(zhuǎn)化為正常值

wind_data.WS = wind_data.WS/10.0

# 瀏覽wind_data數(shù)據(jù)變量的前5行的數(shù)值

wind_data.head()

# 加載日期和時間函數(shù)庫，處理與日期和時間有關(guān)的變量或數(shù)值

from datetime import datetime

# 定義一個空列表變量，用于存儲時間變量

datetimes = []

# 定義year_v變量存放原數(shù)據(jù)中的Year數(shù)據(jù)列

year_v = wind_data['Year'].values

# 定義month_v變量存放原數(shù)據(jù)中的Month數(shù)據(jù)列

month_v = wind_data['Month'].values

# 定義day_v變量存放原數(shù)據(jù)中的Day數(shù)據(jù)列

day_v = wind_data['Day'].values

# 定義hour_v變量存放原數(shù)據(jù)中的Hour數(shù)據(jù)列? ??

hour_v = wind_data['Hour'].values

# 下面的循環(huán)是整合各時間要素變量，每個時間值都是“年月日時”

# List數(shù)據(jù)類型的append方法就是把每個時間值串起來形成一個列表

for?i?in?range(0,?len(year_v)): datetimes.append(datetime(year_v[i], month_v[i], day_v[i], hour_v[i]))

# 把整合后的時間變量作為導(dǎo)入數(shù)據(jù)的行索引，在氣象數(shù)據(jù)的處理上很有用，因?yàn)闅庀髷?shù)據(jù)多是時間序列數(shù)據(jù)

wind_data.index = pd.Series(datetimes)

# 給這個剛定義的索引命名

wind_data.index.name = 'Datetime'

# 瀏覽前5行的數(shù)值

wind_data.head()

# 根據(jù)年和月進(jìn)行分組統(tǒng)計(jì)，WD代表風(fēng)向，WS代表風(fēng)速

wind_data[['Year','Month','WD', 'WS']].groupby(['Year', 'Month']).describe().head()

# 下面繪制圖形，看看有沒有奇異值(outliers)

# 首先創(chuàng)建一個新的繪圖區(qū)

fg = plt.figure()

# 創(chuàng)建組圖中的第一個圖

plt.subplot(121)plt.plot(wind_data.WD)plt.ylabel('Values')plt.title('WD')

# 創(chuàng)建組圖中的第二個

plt.subplot(122)plt.plot(wind_data.WS)plt.ylabel('Values')plt.title('WS')

# 組圖布局設(shè)置

plt.subplots_adjust(top=0.92, bottom=0.41, left=0.10, right=0.95, hspace=0.25, wspace=0.55)

# 從上圖可以看到有很多異常值。定義異常值：第一個是風(fēng)向角度超過360°的，第二個是風(fēng)速異常大的超過1000的

outlier_index = np.logical_or(wind_data['WD'] > 360, wind_data['WS'] > 1000)wind_data = wind_data.loc[~outlier_index, :]print("{0} rows of data were removed".format(len(outlier_index)) )17544 rows of data were removed

# 對剔除了異常值后的數(shù)據(jù)進(jìn)行再次統(tǒng)計(jì)，可以看出，風(fēng)向和風(fēng)速統(tǒng)計(jì)值正常了

wind_data[['Year','Month','WD', 'WS']].groupby(['Year', 'Month']).describe().head()

# 再畫圖看看，這回?cái)?shù)值都正常了

fg = plt.figure()plt.subplot(121)plt.plot(wind_data['WD'],?linewidth=0.2)plt.ylabel('Values')plt.title('WD')plt.subplot(122)plt.plot(wind_data['WS'])plt.ylabel('Values')plt.title('WS')plt.subplots_adjust(top=0.92,?bottom=0.41,?left=0.10,?right=0.95,?hspace=0.25,?wspace=0.55)

# 下面以2012年數(shù)據(jù)為例，我們做一些分析和統(tǒng)計(jì)

# 首先，選擇2012年的5列數(shù)據(jù)組成一個新的Pandas.DataFrame數(shù)據(jù)變量

wind_2012 = wind_data.loc[wind_data['Year']==2012, ['Month','Day','Hour', 'WD', 'WS']]

# 然后繪制畫箱式圖，其目的是查看各月份風(fēng)向和風(fēng)速數(shù)據(jù)的統(tǒng)計(jì)分布特征(中值、分布對稱性、奇異值、分布范圍、集中性)；使用了Seaborn繪圖函數(shù)庫的boxplot函數(shù)繪制該變量的統(tǒng)計(jì)箱式圖

fg = plt.figure()plt.subplot(121)wd_2012 = wind_data[['Month','WD']]sns.boxplot(x = 'WD', y = 'Month', data = wd_2012, orient = 'h', fliersize=0.5, linewidth = 0.7)plt.subplot(122)ws_2012 = wind_data[['Month','WS']]ax = sns.boxplot(x = 'WS', y = 'Month', data = ws_2012, orient = 'h', fliersize=0.5, linewidth = 0.7)ax.set_xlabel('WS (m/s)')plt.subplots_adjust(top=0.92, bottom=0.21, left=0.10, right=1.0, hspace=0.25, wspace=0.45)

# 數(shù)據(jù)統(tǒng)計(jì)分布曲線的繪制；使用Seaborn庫的distplot函數(shù)繪制2012年風(fēng)向的統(tǒng)計(jì)分布曲線

fg = plt.figure()plt.subplot(121)ax?=?sns.distplot(wind_2012.WD)plt.ylabel('Probability')plt.grid()plt.subplot(122)ax?=?sns.distplot(wind_2012.WS)plt.ylabel('Probability')plt.xlabel('WS?(m/s)')plt.grid()plt.subplots_adjust(top=0.92, bottom=0.21, left=0.10, right=1.0, hspace=0.25, wspace=0.45)

# 比較2012年1月和7月的風(fēng)場；分別提取1月和7月的風(fēng)向、風(fēng)速數(shù)據(jù)

fg = plt.figure()wd_Jan?=?wind_2012.loc[wind_2012['Month']?==?1,?['Month','WD']]wd_Jul = wind_2012.loc[wind_2012['Month'] == 7, ['Month','WD']]plt.subplot(121)sns.distplot(wd_Jan.WD)sns.distplot(wd_Jul.WD)plt.ylabel('Probability')plt.grid()plt.legend(['Jan','Jul'], frameon = False)ws_Jan?=?wind_2012.loc[wind_2012['Month']?==?1,?['Month','WS']]ws_Jul = wind_2012.loc[wind_2012['Month'] == 7, ['Month','WS']]plt.subplot(122)sns.distplot(ws_Jan.WS)sns.distplot(ws_Jul.WS)plt.ylabel('Probability')plt.xlabel('WS?(m/s)')plt.grid()plt.legend(['Jan','Jul'], frameon = False)plt.subplots_adjust(top=0.92,?bottom=0.21,?left=0.10,?right=1.0,?hspace=0.25,?wspace=0.45)

# 繪制風(fēng)玫瑰

# 為了繪制風(fēng)玫瑰，需要加載風(fēng)玫瑰繪圖函數(shù)庫

from?windrose?import?WindroseAxesimport matplotlib.cm as cm

# 繪制1月的風(fēng)玫瑰，圖上數(shù)值為百分比；該圖采用堆積柱狀圖表示不同風(fēng)速區(qū)間的百分比

ax?=?WindroseAxes.from_ax()ax.bar(wd_Jan.WD,?ws_Jan.WS,?normed=True,?opening=0.8,?edgecolor='black')ax.set_legend()

# 繪制1月的風(fēng)玫瑰，圖上數(shù)值為在各風(fēng)段區(qū)間出現(xiàn)的風(fēng)速的總數(shù)；采用等值線圖繪制風(fēng)玫瑰，風(fēng)速被劃分為9個等份，contourf為填充圖；采用等值線圖繪制風(fēng)玫瑰，風(fēng)速被劃分為9個等份，contour為線圖

ax?=?WindroseAxes.from_ax()ax.contourf(wd_Jan.WD,?ws_Jan.WS,?bins=np.arange(0,9,1),?cmap=cm.hot)ax.contour(wd_Jan.WD,?ws_Jan.WS,?bins=np.arange(0,9,1),?colors?=?'black',?lw=1)ax.set_title('January')ax.set_legend()

# 繪制7月的風(fēng)玫瑰，圖上數(shù)值為在各風(fēng)段區(qū)間出現(xiàn)的風(fēng)速的總數(shù)；演示另一種繪圖函數(shù)plot_windrose

df?=?pd.DataFrame({'speed':ws_Jul.WS.values,?'direction':wd_Jul.WD.values})ax?=?plot_windrose(df,?kind?=?'contour',?bins=np.arange(0,?9,?1),?cmap=cm.jet,?lw=1)ax.set_title('July')

# 繪制1月和7月堆積柱狀圖式的風(fēng)玫瑰組圖

fg?=?plt.figure()ax1?=?fg.add_subplot(1,2,1,?projection='windrose')ax1.bar(wd_Jan.WD,?ws_Jan.WS,?normed=True,?opening=0.8,?edgecolor='white')ax1.set_title('January', pad = 20)ax2?=?fg.add_subplot(1,2,2,?projection='windrose')ax2.bar(wd_Jul.WD,?ws_Jul.WS,?normed=True,?opening=0.8,?edgecolor='white')ax2.set_title('July', pad = 20)plt.subplots_adjust(top=1, bottom=0, left=0.0, right=1.0, hspace=0.25, wspace=0.45)

# 繪制1月和7月等值線圖式的風(fēng)玫瑰組圖

fg?=?plt.figure()ax_Jan?=?fg.add_subplot(1,2,1,?projection='windrose')ax_Jan.contourf(wd_Jan.WD,?ws_Jan.WS,?bins?=?np.arange(0,9,1.5),?cmap?=?cm.hot)ax_Jan.contour(wd_Jan.WD,?ws_Jan.WS,?bins?=?np.arange(0,9,1.5),?colors?=?'black',?lw?=?0.5)ax_Jan.set_title('January',?pad?=?20)ax_Jul?=?fg.add_subplot(1,2,2,?projection?=?'windrose')ax_Jul.contour(wd_Jul.WD,?ws_Jul.WS,?bins?=?np.arange(0,9,1.5),?colors?=?'black',?lw?=?0.5)ax_Jul.set_title('July', pad = 20)plt.subplots_adjust(top=1, bottom=0, left=0.0, right=1.0, hspace=0.25, wspace=0.45)

# 按16方位風(fēng)向繪制風(fēng)向頻次統(tǒng)計(jì)柱狀圖

ax_Jan.bar(wd_Jan.WD,?ws_Jan.WS,?normed=True,?nsector=16)table_Jan?=?ax_Jan._info['table']wd_freq_Jan?=?np.sum(table_Jan,?axis=0)direction_Jan = ax_Jan._info['dir']plt.figure(figsize=(10,4))plt.bar(np.arange(16),?wd_freq_Jan,?align='center',?facecolor?=?'green')plt.xticks(np.arange(16))xlabels?=?['N','NNE','NE','ENE','E','ESE','SE','SSE','S','SSW','SW','WSW','W','WNW','NW','NNW']_?=?plt.gca().set_xticklabels(xlabels)_?=?plt.text(x?=?14,?y?=?17.5,?s?=?'January')_?=?plt.ylabel('Frequency?(%)')

# 有風(fēng)速觀測數(shù)據(jù)，反推風(fēng)速統(tǒng)計(jì)分布—Weibull統(tǒng)計(jì)分布參數(shù)；params為輸出的Weibull統(tǒng)計(jì)分布參數(shù)，包括shape、loc、scale

from?windrose?import?WindAxesax?=?WindAxes.from_ax()bins?=?np.arange(0,?9+1,?0.5)bins = bins[1:]ax,?params?=?ax.pdf(ws_Jan.WS,?bins?=?bins)ax.set_xlabel('WS?(m/s)')ax.set_ylabel('Probability')print(params)(1, 2.6614301902227524, 0, 3.734483068109758)

# 采用scipy函數(shù)庫中的統(tǒng)計(jì)函數(shù)庫來擬合風(fēng)速的統(tǒng)計(jì)分布參數(shù)

from?scipy?import?statsfg?=?plt.figure()plt.subplot(121)

# 對觀測的風(fēng)速數(shù)據(jù)的統(tǒng)計(jì)分布進(jìn)行Weibull統(tǒng)計(jì)分布擬合，獲取統(tǒng)計(jì)分布參數(shù)

shape, loc, scale = stats.weibull_min.fit(ws_Jan.WS, floc=0)x = np.arange(0, 10, 0.01)

# 利用擬合獲得的參數(shù)，繪制Weibull統(tǒng)計(jì)分布曲線

sns.lineplot(x, stats.weibull_min.pdf(x, c = shape, scale = scale, loc = loc))

# 另一種方法是直接繪制觀測的風(fēng)速數(shù)據(jù)的Weibull統(tǒng)計(jì)分布擬合曲線，這個方法不返回統(tǒng)計(jì)分布參數(shù)

sns.distplot(ws_Jan.WS, bins = 30, fit = stats.weibull_min, kde = False, hist = True)plt.xlabel('WS?(m/s)')plt.ylabel('Probability')plt.title('January')plt.legend(['Weibull 1','Weibull 2'], frameon = False)plt.subplot(122)shape,?loc,?scale?=?stats.weibull_min.fit(ws_Jul.WS,?floc=0)x?=?np.arange(0,?10,?0.01)sns.lineplot(x,?stats.weibull_min.pdf(x,?c?=?shape,?scale?=?scale,?loc?=?loc))sns.distplot(ws_Jul.WS,?bins?=?30,?fit?=?stats.weibull_min,?kde?=?False,?hist?=?True)plt.xlabel('WS?(m/s)')plt.ylabel('Probability')plt.title('July')plt.legend(['Weibull?1','Weibull?2'],?frameon?=?False)plt.subplots_adjust(top=0.92, bottom=0.21, left=0.10, right=1.0, hspace=0.25, wspace=0.45)

# 風(fēng)向的矢量分解，調(diào)用了Ute_WRF的庫函數(shù)

# https://github.com/blaylockbk/Ute_WRF/tree/master/functions

# 庫函數(shù)wind_spddir_to_uv的定義

def?wind_spddir_to_uv(wspd,?wdir):????"""????calculated?the?u?and?v?wind?components?from?wind?speed?and?direction????Input:????????wspd:?wind?speed????????wdir:?wind?direction????Output:????????u:?u?wind?component????????v:?v?wind?component????"""????????rad?=?4.0*np.arctan(1)/180.????u?=?-wspd*np.sin(rad*wdir)????v?=?-wspd*np.cos(rad*wdir)return u,v

# 設(shè)置繪圖環(huán)境

%matplotlib

# 矢量分解函數(shù)的調(diào)用

u1, v1 = wind_spddir_to_uv(ws_Jan.WS, wd_Jan.WD)

# 繪制矢量動畫，即顯示不同時刻風(fēng)矢量的方向和大小

plt.ionfg?=?plt.figure()ax?=?fg.add_subplot(111)for?i?in?range(0,?len(u1)):????#?清楚當(dāng)前繪圖區(qū)????plt.clf()????#?矢量繪圖函數(shù)quiver????plt.quiver(-u1[i],?-v1[i],?scale=20.0,?color?=?'blue')????#?設(shè)置橫坐標(biāo)范圍????plt.xlim([-2,2])????#?設(shè)置縱坐標(biāo)范圍????plt.ylim([-2,2])????#?設(shè)置圖中文字信息，顯示年月日時????plt.text(x?=?0.6,?y?=?0.0,?s='{0}-{1}-{2}-{3}'.format(????????????ws_Jan.index[i].year,?ws_Jan.index[i].month,?ws_Jan.index[i].day,?????????????ws_Jan.index[i].hour))????#?設(shè)置圖中文字信息，顯示(U,V)????plt.text(x?=?0.6,?y?=?-1.5,?s='(U,V)?=?({0},{1})'.format(round(u1[i],2),?round(v1[i],2)))????#?設(shè)置圖中文字信息，顯示風(fēng)向(度)????plt.text(x?=?0.6,?y?=?-0.5,?s='Wind?direction:?{0}'.format(round(wd_Jan.WD[i],?2)))????#?設(shè)置圖中文字信息，顯示風(fēng)速(m/s)????plt.text(x?=?0.6,?y?=?-1.0,?s='Wind?speed:?{0}?(m/s)'.format(round(ws_Jan.WS[i],?2)))????#?保存每一個圖為一個png格式的圖像文件????plt.savefig('{0}.png'.format(i))????fg.canvas.draw()????#?刷新繪圖區(qū) fg.canvas.flush_events()

# 風(fēng)矢量的GIF動畫輸出

# 加載glob函數(shù)庫和imageio圖像輸入輸出函數(shù)庫

import globimport imageio

# 輸出GIF

graphs?=?[]for?i?in?range(0,?len(u1)):????#?filename記錄的是當(dāng)前目錄下所有風(fēng)矢量圖????filename?=?str('{0}.png'.format(i))????#?把所有風(fēng)矢量圖讀取到內(nèi)存，并串起來形成一個列表????graphs.append(imageio.imread(filename))#?定義輸出的GIF文件名output_file?=?'Gif_{0}.gif'.format(datetime.today())#?使用mimsave函數(shù)輸出動畫imageio.mimsave(output_file,?graphs,?duration?=?0.1,?loop?=?1)

上述代碼在Python 3.6以上測試通過，源程序可以訪問作者的GitHub主頁下的Meteo_data_analysis：https://github.com/yangsbin/Meteo_data_analysis

如果對本案例有什么補(bǔ)充和建議，歡迎留言?。謝謝！

往期回顧：

氣象招聘 | 中國民用航空華東地區(qū)空中交通管理局招聘氣象崗位宣講會

氣象招聘 | 天象科技2019秋冬季人才招聘計(jì)劃

行業(yè)動態(tài) | “遠(yuǎn)征杯”中國信用防雷宣傳活動開鑼了！

2020屆氣象本科考研與就業(yè)QQ群：9301072622020屆氣象研究生畢業(yè)就業(yè)QQ群：3465774272021屆氣象本科考研與就業(yè)QQ群：639522239

總結(jié)

以上是生活随笔為你收集整理的python subplot_气象编程 | 一个简单的风数据处理和分析案例（Python版）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。