最新热点文章：行业动态 | 2020中国气象现代化建设科技博览会
行业动态 | “远征杯”中国信用防雷宣传活动开锣了！
气象招聘 | 中南空管局气象中心2020年校园招聘公告
今天我们来讲讲风数据的简单处理和分析。通过这个案例，你将知道一些与风数据处理和分析有关的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
