基于Python的三维风场风玫瑰图绘制方法研究与应用

周积强，魏 巍，姚肖萌，王雪妮，李龙艳，杨 勇，柳佳俊

(1.宁夏回族自治区气象灾害防御技术中心，银川 750002；2.银川市气象局，银川 750002；3.宁夏气象防灾减灾重点实验室，银川 750002)

0 引言

空气的水平运动称为风，不仅有数值的大小，还具有方向。地面风向一般用16方位表示，每相邻方位间的角差为22.5°。风玫瑰图是统计和显示风的专业方法，是指极坐标底图上绘制出某一地区在某一时间段内各风向出现的频率或各风向的平均风速的统计图，相应分为风向玫瑰图和风速玫瑰图。在真实大气中风还包括垂直运动，即上升或下沉气流。三维风速仪能够实时观测水平、垂直风速及风向，广泛应用于测量风场、湍流等观测研究中。

风速数据处理和玫瑰图的绘制有很多方式。常见的有使用专用的风统计软件、Excel[1]、VB[2]、VC[3]、C#[4]、Extjs[5]、Python[6]等进行绘制。Python具有强大的数据处理和分析、可视化、丰富的标准库及第三方库的特点，目前已在气象领域得到广泛应用[7]。但使用Python对大气三维风绘制玫瑰图仍然较少，包括水平风和垂直气流风速风向玫瑰图。

根据风速风向玫瑰图可视化展示特点，文章分析了三维风数据和数据结构存在的问题，使用三倍标准差等方法对数据进行了质量控制。基于Python平台使用第三方库Windrose库，设计了一种大气三维风包括水平风和垂直气流风速风向玫瑰图的绘制方法。使用实测三维风资料进行处理和应用，结果表明绘制的水平风和垂直气流风玫瑰图能够清晰地显示出风的三维风统计特征，为近地层大气三维风场研究提供了一种统计方法。

1 资料来源和方法

文章使用的三维风速仪是超声波风速风向仪，能够记录水平风速、垂直气流、风向等信息。数据保存格式为文本格式。采样数据存在一些问题，主要包括存储数据包含一些将风速风向数据记录为NAN记录。有些数据记录为奇异值，还有重复记录，如同一时间会出现两条以上的数据记录。极少情况下水平风数据出现负号。最后数据在使用前需进行粗大值筛选和剔除。使用3倍标准差法将粗大值剔除和插补。查找重复数据记录，仅保留首次记录的数据。将出现的个别空值使用邻近法进行插补，从而得到经数据清洗后的数据集。

2 数据分析流程和关键技术

文章使用Python第三方库Windrose来绘制三维风(水平风和垂直风)风速风向玫瑰图。仅安装Windrose第三方库并不能顺利使用其函数进行风玫瑰图绘制，同时还需要安装相关的依赖库，如Pandas，Numpy，Matplotlib，Scipy等。

2.1 数据分析流程

三维风数据存储格式为文本格式，在多个文档中存储。处理之前需选择待分析的多个文档合并、数据清洗。然后将需要绘制的风向、水平风、垂直气流信息提取出，使用Windrose第三方库分别完成水平风和垂直气流速度的风速风向玫瑰图。数据分析流程如图1所示。

2.2 关键技术

2.2.1 数据预处理

使用os和tkinter库函数选择需要处理的数据文件。读取数据使用函数pd.read_csv()。将循环读取的多个DataFrame表分别追加暂存入列表，使用concat()函数将列表内的DataFrame表进行垂直合并。按照数据质量控制方法使用Pandas库中的replace()函数将数据表中的NAN信息替换为程序能够识别的np.nan。使用astypes()函数强制转换数据类型。待数据空值和类型全部确定后，使用3倍标准差方式分别筛选出水平风和垂直气流速度的粗大值并将其定义为空值，将出现水平风为负值的数据定义为空值。最后将筛选出的空值使用就近填充法对其进行填充。

图1 数据处理流程

2.2.2 绘制方法

直接使用Windrose中WindroseAxes()函数绘制风速风向玫瑰图时，其显示内容不符合使用习惯和需求，要将风向信息由度数调整为风向表示方式。绘制出的0°方向指向右侧，需要将其旋转设置为正上方为N方向或0°方向。数据旋转方式为逆时针方向，需将风向旋转方向调整为顺时针方向。最重要的是直接使用显示数据结果存在错误，如将0°方向绘制在90°方向，将90°方向绘制在0°方向。使用ax.set_theta_zero_location('N') 函数调整方位，使用ax.set_theta_direction('clockwise')将数据旋转方式调整为顺时针方式显示。对风向数据进行转换处理，经调整后数据显示为正确结果，如图2所示。

图2 调整后风速风向玫瑰图

水平风速有明确的风力等级划分标准，因此在图例设置时按照风力等级进行分档和显示。但垂直气流速度没有标准可以参照划分，并且垂直气流运动上升和下沉时常交替出现，研究人员更关注垂直气流是上升还是下沉，以及其发展趋势。另外垂直气流速度一般情况下较小，经常为小于1.0 m/s，因此根据实际工作需要，将垂直气流上升或下沉分开显示，并分别具有多个分档。另外为了便于分析，程序还设计了8方位、16方位和带有中文或英文两种不同显示风向方式的玫瑰图。绘制方法程序详细如下：

风向数据转换程序：

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

风速风向玫瑰图绘制程序：

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='风速(单位：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 应用个例

文章使用数据的时间段为2021-07-01—2021-07-10。按照设计的程序，对个例数据进行测试应用。

3.1 水平风速风向玫瑰图

从水平风16方位风速风向玫瑰图统计结果可得出：7月上旬六盘山东坡山腰位置水平风主要风向为南风、西南偏南、西南风、东南偏南，风频占比分别为14.0%，9.4%，8.8%，9.0%左右。其中南风风向中风速在0.3～1.6 m/s时，风频占比为2.9%；风速在1.6～3.4 m/s时，风频占比为7.1%；风速在3.4～5.5 m/s时，风频占比为4.0%。从8方位风速风向玫瑰图中可以看出水平风主要风向为南风、东南、西南，风频占比分别为24.0%，13.0%，16.0%左右。其中南风风向中风速在0.3～1.6 m/s时，风频占比为6.0%；风速在1.6～3.4 m/s时，风频占比为12.0%；风速在3.4～5.5 m/s时，风频占比为6.0%。图3为16方位水平风速风向玫瑰图。

3.2 垂直气流风玫瑰图

从垂直气流16方位风速风向玫瑰图统计结果可得出：7月上旬在六盘山东坡山腰位置，当风向为东南风时，垂直气流表现为上升运动，当风向为西南、西风、西北时，垂直气流表现为下沉和上升运动均存在。其中当西风时，下沉气流占比大于下沉气流。其中当风向为南风时，下沉气流风速在-2.7～0.0 m/s时，风频占比为5.8%；上升气流风速区间在0.0～1.1 m/s时，风频占比为8.2%；上升气流风速在1.1～2.1 m/s时，风频占比为0.3%。当风向为东南风时，垂直气流大部分为上升运动，下沉气流风速在-2.7～0.0 m/s时，风频占比为1.0%；上升气流风速在0.0～1.1 m/s，风频占比为4.7%；上升气流风速在1.1～2.1 m/s，风频占比为6.2%。从8方位风速风向玫瑰图中可更加明显地看出大气垂直运动上升气流与风向的关系。当风向为南风时，垂直气流下沉运动风频占比为10.1%，上升运动风频占比为14.0%。当风向为东南风时，垂直气流下沉运动风频占比为2.0%，上升运动风频占比为12.0%。当风向为东风时，垂直气流下沉运动风频占比小于1.0%，上升运动风频占比为9.0%。因六盘山山脉总体呈南北走向，当近地层风从西向东吹时，近地层大气翻越山体后在山体东侧下沉，观测得到的垂直气流表现为下沉运动。当近地层风从东向西吹时，近地层大气由于山体强迫抬升，观测得到的垂直气流表现为上升运动。16方位垂直气流风速风向玫瑰图如图4所示。

图3 16方位水平风速风向玫瑰图

图4 16方位垂直气流风速风向玫瑰图

4 结束语

大气除了水平运动以外还有大气的垂直运动，三维风速仪能够观测大气三维运动情况。绘制三维风玫瑰图对分析研究大气三维运动具有重要作用。文章分析了三维风速仪资料特点，基于Python平台使用第三方Windrose函数设计了三维风玫瑰图绘制程序，并使用该程序绘制宁夏六盘山地区东坡山腰实际观测三维风玫瑰图，包括水平风风速风向玫瑰图和垂直气流风速风向玫瑰图。文章解决了直接使用Windrose库函数绘制风速风向玫瑰图存在风向0°方向没有指向北方、风向旋转方式为逆时针，显示数据存在错误等问题。绘制得到水平风和垂直气流风速风向玫瑰图能够清晰地反应出当地近地层大气运动情况及与地形之间的关系，为相关研究和玫瑰图自动绘制提供了一种方法。