python subplot_气象编程 | 一个简单的风数据处理和分析案例（Python版）
今天我們來講講風數據的簡單處理和分析。通過這個案例，你將知道一些與風數據處理和分析有關的Python函數庫和繪圖方法。具體內容如下：
1)如何從網絡上獲取數據并讀取格式化數據(reading data)
2)如何處理時間變量(datetime variable)
3)如何判別異常值和去除異常值(Outlier)
4)如何使用Pandas.DataFrame來篩選、整合、選擇數據
5)如何繪制風玫瑰圖
6)如何繪制數據的統計分布曲線
7)如何從數據中獲取Weibull分布參數
8)如何繪制16位風向統計柱狀圖
9)如何繪制單點風矢量動畫和GIF輸出
10)如何繪制箱式統計圖(boxplot)
11)風的矢量分解函數的使用
Python編程語言已經成為世界上最流行和受歡迎的編程語言。它的靈活性和易用性已經強大的函數庫吸引了數據分析愛好者。今天我們一起來看一例使用Python語言來處理和分析風數據的一個案例。希望這個案例能夠幫助到大家。
風數據是我們應用氣象專業的同行們經常要接觸的。相比溫度、降水量、太陽輻射等氣象要素，風要素的重要特征是它是一個矢量要素，有方向有大小。風方向從正北方向為0°開始，順時針轉向360°，因此正東方向就對應風向90°，正南為180°，正西為270°，正北為0°。在氣象上通常將風向劃分為16個方向，用符號N，NNE，NE，...來表示。風速則根據風速大小劃分等級，比如國際通用的風力等級：蒲福風級(https://baike.baidu.com/item/蒲福風級)。因此有了風向、風速的分類分級，我們就可以對風要素進行分析，獲取不同地方的風場特征。其中最重要的風場特征表現形式是風玫瑰圖。下面我們就風數據的讀取、預處理和分析來做一個案例。
# 加載一些常用的Python函數庫
# numpy是支持多維數據存儲、運算、分析函數庫
# pandas是結構化數據的存儲、運算、分析函數庫
# matplotlib是比用繪圖函數庫??
# seaborn是繪圖函數庫??
# warnings是運行警告函數庫
import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as snsimport warnings# 運行環境設置：在Notebook中抑制警告輸出? ? ? ? ? ? ? ??
warnings.filterwarnings('ignore')# 在notebook中顯示繪制的圖像
%matplotlib# 加載風數據。風數據存放在網絡上。注：本案例的風數據是在原始數據基礎上增加了噪聲
url?=?'https://raw.githubusercontent.com/yangsbin/Meteo_data_analysis/master/data/wind_data.txt'# 通過Pandas的read_csv函數讀取格式化數據，參數中sep是讀取的原數據列的分隔符參數定義， header是頭行(None就是無頭行)，names定義了每個數據列的名稱；WD定義為風向；WS定義為風速
wind_data?=?pd.read_csv(url,?sep?=?r'\s+',?header?=?None,?names?=?['Year','Month','Day','Hour',?'WD','WS'])# 原始數據中風速是10倍值，因此將讀取后的數據中的風速列數據除10轉化為正常值
wind_data.WS = wind_data.WS/10.0# 瀏覽wind_data數據變量的前5行的數值
wind_data.head()# 加載日期和時間函數庫，處理與日期和時間有關的變量或數值
from datetime import datetime# 定義一個空列表變量，用于存儲時間變量
datetimes = []# 定義year_v變量存放原數據中的Year數據列
year_v = wind_data['Year'].values# 定義month_v變量存放原數據中的Month數據列
month_v = wind_data['Month'].values# 定義day_v變量存放原數據中的Day數據列
day_v = wind_data['Day'].values# 定義hour_v變量存放原數據中的Hour數據列? ??
hour_v = wind_data['Hour'].values# 下面的循環是整合各時間要素變量，每個時間值都是“年月日時”
# List數據類型的append方法就是把每個時間值串起來形成一個列表
for?i?in?range(0,?len(year_v)): datetimes.append(datetime(year_v[i], month_v[i], day_v[i], hour_v[i]))# 把整合后的時間變量作為導入數據的行索引，在氣象數據的處理上很有用，因為氣象數據多是時間序列數據
wind_data.index = pd.Series(datetimes)# 給這個剛定義的索引命名
wind_data.index.name = 'Datetime'# 瀏覽前5行的數值
wind_data.head()# 根據年和月進行分組統計，WD代表風向，WS代表風速
wind_data[['Year','Month','WD', 'WS']].groupby(['Year', 'Month']).describe().head()# 下面繪制圖形，看看有沒有奇異值(outliers)
# 首先創建一個新的繪圖區
fg = plt.figure()# 創建組圖中的第一個圖
plt.subplot(121)plt.plot(wind_data.WD)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)# 從上圖可以看到有很多異常值。定義異常值：第一個是風向角度超過360°的，第二個是風速異常大的超過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# 對剔除了異常值后的數據進行再次統計，可以看出，風向和風速統計值正常了
wind_data[['Year','Month','WD', 'WS']].groupby(['Year', 'Month']).describe().head()# 再畫圖看看，這回數值都正常了
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年數據為例，我們做一些分析和統計
# 首先，選擇2012年的5列數據組成一個新的Pandas.DataFrame數據變量
wind_2012 = wind_data.loc[wind_data['Year']==2012, ['Month','Day','Hour', 'WD', 'WS']]# 然后繪制畫箱式圖，其目的是查看各月份風向和風速數據的統計分布特征(中值、分布對稱性、奇異值、分布范圍、集中性)；使用了Seaborn繪圖函數庫的boxplot函數繪制該變量的統計箱式圖
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)# 數據統計分布曲線的繪制；使用Seaborn庫的distplot函數繪制2012年風向的統計分布曲線
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月的風場；分別提取1月和7月的風向、風速數據
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)# 繪制風玫瑰
# 為了繪制風玫瑰，需要加載風玫瑰繪圖函數庫
from?windrose?import?WindroseAxesimport matplotlib.cm as cm# 繪制1月的風玫瑰，圖上數值為百分比；該圖采用堆積柱狀圖表示不同風速區間的百分比
ax?=?WindroseAxes.from_ax()ax.bar(wd_Jan.WD,?ws_Jan.WS,?normed=True,?opening=0.8,?edgecolor='black')ax.set_legend()# 繪制1月的風玫瑰，圖上數值為在各風段區間出現的風速的總數；采用等值線圖繪制風玫瑰，風速被劃分為9個等份，contourf為填充圖；采用等值線圖繪制風玫瑰，風速被劃分為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月的風玫瑰，圖上數值為在各風段區間出現的風速的總數；演示另一種繪圖函數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月堆積柱狀圖式的風玫瑰組圖
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月等值線圖式的風玫瑰組圖
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方位風向繪制風向頻次統計柱狀圖
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?(%)')# 有風速觀測數據，反推風速統計分布—Weibull統計分布參數；params為輸出的Weibull統計分布參數，包括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函數庫中的統計函數庫來擬合風速的統計分布參數
from?scipy?import?statsfg?=?plt.figure()plt.subplot(121)# 對觀測的風速數據的統計分布進行Weibull統計分布擬合，獲取統計分布參數
shape, loc, scale = stats.weibull_min.fit(ws_Jan.WS, floc=0)x = np.arange(0, 10, 0.01)# 利用擬合獲得的參數，繪制Weibull統計分布曲線
sns.lineplot(x, stats.weibull_min.pdf(x, c = shape, scale = scale, loc = loc))# 另一種方法是直接繪制觀測的風速數據的Weibull統計分布擬合曲線，這個方法不返回統計分布參數
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)# 風向的矢量分解，調用了Ute_WRF的庫函數
# https://github.com/blaylockbk/Ute_WRF/tree/master/functions
# 庫函數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# 設置繪圖環境
%matplotlib# 矢量分解函數的調用
u1, v1 = wind_spddir_to_uv(ws_Jan.WS, wd_Jan.WD)# 繪制矢量動畫，即顯示不同時刻風矢量的方向和大小
plt.ionfg?=?plt.figure()ax?=?fg.add_subplot(111)for?i?in?range(0,?len(u1)):????#?清楚當前繪圖區????plt.clf()????#?矢量繪圖函數quiver????plt.quiver(-u1[i],?-v1[i],?scale=20.0,?color?=?'blue')????#?設置橫坐標范圍????plt.xlim([-2,2])????#?設置縱坐標范圍????plt.ylim([-2,2])????#?設置圖中文字信息，顯示年月日時????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))????#?設置圖中文字信息，顯示(U,V)????plt.text(x?=?0.6,?y?=?-1.5,?s='(U,V)?=?({0},{1})'.format(round(u1[i],2),?round(v1[i],2)))????#?設置圖中文字信息，顯示風向(度)????plt.text(x?=?0.6,?y?=?-0.5,?s='Wind?direction:?{0}'.format(round(wd_Jan.WD[i],?2)))????#?設置圖中文字信息，顯示風速(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()????#?刷新繪圖區 fg.canvas.flush_events()# 風矢量的GIF動畫輸出
# 加載glob函數庫和imageio圖像輸入輸出函數庫
import globimport imageio# 輸出GIF
graphs?=?[]for?i?in?range(0,?len(u1)):????#?filename記錄的是當前目錄下所有風矢量圖????filename?=?str('{0}.png'.format(i))????#?把所有風矢量圖讀取到內存，并串起來形成一個列表????graphs.append(imageio.imread(filename))#?定義輸出的GIF文件名output_file?=?'Gif_{0}.gif'.format(datetime.today())#?使用mimsave函數輸出動畫imageio.mimsave(output_file,?graphs,?duration?=?0.1,?loop?=?1)上述代碼在Python 3.6以上測試通過，源程序可以訪問作者的GitHub主頁下的Meteo_data_analysis：https://github.com/yangsbin/Meteo_data_analysis
如果對本案例有什么補充和建議，歡迎留言?。謝謝！
往期回顧：
氣象招聘 | 中國民用航空華東地區空中交通管理局招聘氣象崗位宣講會
氣象招聘 | 天象科技2019秋冬季人才招聘計劃
行業動態 | “遠征杯”中國信用防雷宣傳活動開鑼了！
2020屆氣象本科考研與就業QQ群：9301072622020屆氣象研究生畢業就業QQ群：3465774272021屆氣象本科考研與就業QQ群：639522239總結
以上是生活随笔為你收集整理的python subplot_气象编程 | 一个简单的风数据处理和分析案例（Python版）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： python自学行_python自学行吗
- 下一篇： 筒灯里的模块脱落了会导致筒灯更暗吗？