风玫瑰图整合(一)风速玫瑰图
文章目录
一、前言
全网中关于python绘制风玫瑰图的参考材料较少,基于第一期的风玫瑰图,我们再做一次拓展与整合解析
首先风玫瑰图分为:风向玫瑰图 和 风速玫瑰图。在多数论文中,提及的风玫瑰图一般指风向玫瑰图。其中风速玫瑰图常见有:
风向玫瑰图常见有:
本篇基于第一篇风玫瑰图的图形,主要解析风速玫瑰图。
二、风速玫瑰图怎么看
这是一张风速玫瑰图(a),采用的绘制方法对应极坐标系下的柱状图bar。而右图为,基于左图数据得到的无静风的各风向对应风速的频率。针对绘图的数据,其设置的分割区间为:
0.3—1.5 | 1.6—3.3 | 3.4—5.5 | 5.6—7.9 | 8.0—10.7 | 10.8—— |
---|---|---|---|---|---|
1级 | 2级 | 3级 | 4级 | 5级 | >5级 |
其中小于1级,其实就是低于0.3m/s以下,气象上定义为静风。
现在我们专注于风频最大、最突出的北北东风(NNE),可以看到,代表着1 级风的深蓝色轮辐,占比约6.07%;代表2 级风的蓝色轮辐,占比约7.80%,所处位置接近15%。
也由此得出,2级风的占比是由蓝色轮辐的上限减去下限得到的7.80%。只有对应颜色的部分,才是该方向上的风在某个风速下的风频占比。
所以针对风速玫瑰图,我个人的理解,需要配备一张对应的表格,才能更明确分析各方向上风速出现的普遍情况。因为风向的集中情况,风向玫瑰图便足够直观表现了。
而增加对风速的分析,单一张图只能粗略地得到,某方向上某一风级的突出情况,无法展开更细节的比较。
相较柱状图(a)和等值线图(c),前者更专注于突显大值内的分布变化,如NNE中风级频率最大的2级风与其他风级的对比,或者和相邻风向的对比,会很明了地分清在2级风上,NNE风的频率要比NE风要大。
而后者,则更注重方向上其最常出现的某风级、其频率随总体趋势变化而变化,能更容易明晰各风级总体上的占比面积。
三、各类风速玫瑰图怎么画
数据读取相关函数及类
(1)气象时制列表导出
不做解析,详见python设计风频统计
class meteorological_date2datetime:
def __init__(self, extract_date: datetime.date):
self.date = extract_date
@property
def day_hour(self):
_day_hour = []
_year = self.date.year
_month = self.date.month
_day = self.date.day
_hours_list = [21, 22, 23] + [i for i in range(21)]
# 是否为昨天的判断标识
pre_day = 1
for i in _hours_list:
if pre_day == 1 and datetime.time(i,0) != datetime.time(0,0):
_datetime = datetime.datetime(_year,_month,_day,i,0) - datetime.timedelta(days=1)
_day_hour.append(_datetime)
else:
pre_day = 0
_datetime = datetime.datetime(_year, _month, _day, i, 0)
_day_hour.append(_datetime)
return _day_hour
@property
def day_minute(self):
_day_minute = []
_year = self.date.year
_month = self.date.month
_day = self.date.day
# 预组合时分
_hours_list = [21, 22, 23] + [i for i in range(21)]
_minute_list = [i+1 for i in range(59)] + [0]
# 是否为昨天的判断标识,是0,否1
pre_day = 1
for i in _hours_list:
for j in _minute_list:
# 先判断是否为昨天,循环判断是否经过00时,昨天今天日的分界线。
if datetime.time(i,j) == datetime.time(0,0):
pre_day = 0
_datetime = datetime.datetime(_year, _month, _day, i, 0)
elif pre_day == 0:
# 今天日的时分
if datetime.time(i,j).minute == 0:
_datetime = datetime.datetime(_year, _month, _day, i, 0)
else:
_datetime = datetime.datetime(_year, _month, _day, i, j) - datetime.timedelta(hours=1)
else:
# 昨天日的时分
if datetime.time(i,j).minute == 0:
_datetime = datetime.datetime(_year, _month, _day, i, 0) - datetime.timedelta(days=1)
else:
_datetime = datetime.datetime(_year, _month, _day, i, j) - datetime.timedelta(hours=1,days=1)
_day_minute.append(_datetime)
return _day_minute
(2)excel数据导出
不做解析,详见python设计风频统计
def extract_data(_file_path,name):
df_list = []
date = []
for m in range(12):
year = 2023
month = str(m + 1) + '月'
data = pd.read_excel(_file_path, sheet_name=month, header=0, index_col=0)
day = data.columns.values
for i in day:
_date = datetime.date(year,m+1,i)
date.append(_date)
time_list = meteorological_date2datetime(_date).day_hour
day_data = data[i].replace(-999, float('nan')).values
day_df = pd.DataFrame(day_data, index=time_list)
df_list.append(day_df)
_df = pd.concat(df_list)
_df = _df.rename(columns=dict(zip(_df.columns, [name])))
return _df
(3)将数组转换为百分制
该函数用于表格中将频次转换成频率。
# 将数组转换为百分制
def array_str_percentage(arr):
_values = []
_shape = arr.shape
_n = len(_shape)
for i in range(_shape[0]):
value = arr[i]
# 第一层循环,判断出目前是一维,传入的就是二维
if len(value.shape) == 1:
_values.append(['{:.2f}%'.format(a * 100) for a in value])
else:
_values.append(array_str_percentage(arr[i]))
return np.array(_values)
绘图函数类设置
以下设计的wind_subplots
类,可以说是作为windrose
库的补丁。其中windrose
库本身是支持对数据静风的筛选的。如下windrose源代码:
if calm_limit is not None:
mask = var > calm_limit
self.calm_count = len(var) - np.count_nonzero(mask)
if normed:
self.calm_count = self.calm_count * 100 / len(var)
var = var[mask]
direction = direction[mask]
"""
var 是风速
direction 是风向
calm_limit 是判定静风的最大值
"""
根据以上内容,可以得知,要想使用它本身的静风筛出,要求传入的风向,必须要可以能通过var > calm_limit
这种条件索引的类。
所以可以用pd.Series
,pd.DataFrame
,np.ndarray
。
(0)源码类-风速玫瑰图设置
class wind_subplots:
def __init__(self,
_wd, _ws, _sum_freq,
_C:str=None,
_calm_limit:float=None,
sub_tup:tuple=(1,1),
pic_type:str=None):
"""
:param _wd: 风向数据
:param _ws: 风速数据
:param _sum_freq: 总频次
:param _C: 自算静风率
:param _calm_limit:静风限制,windrose会根据大于该值计算静风频次
:param sub_tup: 子图元组(1,1)或(1,2),默认图1为风玫瑰图,图2为图表(选设)
:param pic_type: 图1选择绘制图种样式
"""
# 基础参数
self.wd = _wd
self.ws = _ws
self.sum_freq = _sum_freq
self.C = _C
self.calm_limit = _calm_limit
if self.calm_limit:
if self.calm_limit <= 0.0:
self.calm_limit = None
if isinstance(self.wd,pd.Series) or isinstance(self.wd, pd.DataFrame) or isinstance(self.wd,np.ndarray):
pass
elif isinstance(self.ws,pd.Series) or isinstance(self.ws, pd.DataFrame) or isinstance(self.wd,np.ndarray):
pass
else:
raise TypeError("如果使用windrose自行判断静风,传入的数据必须是数组。仅支持pd.Series,pd.DataFrame,np.ndarray。")
else:
self.calm_limit = None
# 子图设置
if sub_tup[1] == 2:
self.fig = plt.figure(figsize=(12,6))
self.ax2 = None
else:
self.fig = plt.figure(figsize=(6, 6))
self.nrow = sub_tup[0]
self.ncol = sub_tup[1]
self.pic_type = pic_type
self.ax1 = self.fig.add_subplot(self.nrow,self.ncol,1,projection='windrose')
# 默认参数
self.dir = ['E', 'ENE', 'NE', 'NNE', 'N', 'NNW', 'NW', 'WNW', 'W', 'WSW', 'SW', 'SSW', 'S', 'SSE',
'SE', 'ESE']
self.levels = [' 1 级', ' 2 级', ' 3 级', ' 4 级', ' 5 级', '>5级']
self.colors = ['#0000bc', '#3a93f1', '#40ffe7', '#cbfc7c', '#fdd305', '#f46c06']
self.patches = [mpatches.Patch(color=self.colors[i], label=self.levels[i]) for i in range(len(self.colors))]
self.bins = [0, 1.6, 3.4, 5.5, 8.0, 10.8]
# 预设存放表格数据
self.table_data = []
def draw_windrose(self):
# 设置风玫瑰图偏移画布中心,多出空白放置图例
self.ax1.set_position([0.15, 0.15, 0.75, 0.75]) # [left, bottom, width, height]
# 设置成16方位
self.ax1.set_thetagrids(np.arange(0.0, 360.0, 22.5), labels=self.dir)
if self.pic_type == 'bar':
self.ax1.bar(self.wd, self.ws, opening=0.8, colors=self.colors, bins=self.bins,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
elif self.pic_type == 'contour':
self.ax1.contour(self.wd, self.ws, bins=self.bins, colors=self.colors,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
elif self.pic_type == 'contourf':
self.ax1.contourf(self.wd, self.ws, bins=self.bins, colors=self.colors,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
# 换算每一圈的比例
self.ax1.set_yticks([self.sum_freq * 0.05, self.sum_freq * 0.10, self.sum_freq * 0.15, self.sum_freq * 0.20])
self.ax1.set_yticklabels(['5%', '10%', '15%', '20%'])
# 因为采用原理为极坐标,在极坐标上,90°朝上。因此我们希望将每一圈的比例坐标轴指向北,则选择90°方位
self.ax1.set_rlabel_position(90)
def set_colorbar(self):
# 上设颜色划分为6,色标共6格,7个坐标,从0开始,6最大
norm = mpl.colors.Normalize(vmin=0, vmax=len(self.levels))
cmap = mpl.colors.LinearSegmentedColormap.from_list('How2matplotlib_custom', self.colors, N=len(self.levels))
sm = cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([]) # 防止显示数值
# 创建色标
cbar = plt.colorbar(sm, ax=self.ax1, orientation='horizontal')
cbar.ax.set_xticks(np.arange(0,len(self.levels)+1,1)) # 设置新的刻度位置
cbar.ax.set_xticklabels(['']+self.levels)
def set_legend(self):
# Windroses设置图例方法
self.ax1.legend_ = legend.Legend(parent=self.ax1, handles=self.patches, labels=self.levels,
loc='lower left', bbox_to_anchor=(-0.15, -0.15), fontsize='small')
# 表格显示频率
def add_table(self):
if self.nrow == 2 or self.ncol == 2:
self.ax2 = self.fig.add_subplot(self.nrow,self.ncol,2)
plt.subplots_adjust(wspace=0.2)
rowLabel = self.dir[5:] + self.dir[:5]
_table_data = array_str_percentage(self.table_data / self.sum_freq)
_table_data = np.vstack((self.levels,_table_data))
table = self.ax2.table(cellText=_table_data, loc='center left', rowLoc='center',
rowLabels=[''] + rowLabel[::-1],edges='horizontal')
table.set_fontsize(14)
table.scale(1.2, 1.2)
self.ax2.spines['bottom'].set_visible(False)
self.ax2.spines['left'].set_visible(False)
self.ax2.spines['right'].set_visible(False)
self.ax2.spines['top'].set_visible(False)
self.ax2.get_xaxis().set_ticks([])
self.ax2.get_yaxis().set_ticks([])
# 表格显示频次
def add_table1(self):
if self.nrow == 2 or self.ncol == 2:
self.ax2 = self.fig.add_subplot(self.nrow,self.ncol,2)
plt.subplots_adjust(wspace=0.2)
rowLabel = self.dir[5:] + self.dir[:5]
_table_data = self.table_data.astype(int)
_table_data = np.vstack((self.levels,_table_data))
table = self.ax2.table(cellText=_table_data, loc='center left', rowLoc='center',
rowLabels=[''] + rowLabel[::-1],edges='horizontal')
table.set_fontsize(14)
table.scale(1.2, 1.2)
self.ax2.spines['bottom'].set_visible(False)
self.ax2.spines['left'].set_visible(False)
self.ax2.spines['right'].set_visible(False)
self.ax2.spines['top'].set_visible(False)
self.ax2.get_xaxis().set_ticks([])
self.ax2.get_yaxis().set_ticks([])
def add_calm_circle(self):
self.ax1.set_rorigin(-120)
# 在圆圈中添加文本标注
if self.calm_limit is None:
self.ax1.text(0, -10, self.C,
ha='right', va='center', fontsize=8, color='black')
else:
self.ax1.text(0, -10, f'{round(self.ax1.calm_count * 100 / self.sum_freq, 1)}%',
ha='right', va='center',fontsize=8, color='black')
def add_calm_text(self):
# ax1.transAxes 转换成轴坐标。在轴坐标系中,坐标的原点(0,0)位于图表的左下角,而(1,1)则代表图表的右上角,无论图表的实际数据范围如何。
if self.C is None:
self.C = str(np.round(self.ax1.calm_count * 100 / self.sum_freq,3)) + '%'
self.ax1.text(0.5, -0.1, '静风频率:' + self.C,
horizontalalignment='center', verticalalignment='center', transform=self.ax1.transAxes)
(1)构造函数__init__
基础参数
# 风向
self.wd = _wd
# 风速
self.ws = _ws
# 总频次
self.sum_freq = _sum_freq
# 自算静风率
self.C = _C
# windrose的静风大小设置
self.calm_limit = _calm_limit
子图设置
# 子图设置
# ax1:风玫瑰图;ax2:表格
if sub_tup[1] == 2:
self.fig = plt.figure(figsize=(12,6))
self.ax2 = None
else:
self.fig = plt.figure(figsize=(6, 6))
self.nrow = sub_tup[0]
self.ncol = sub_tup[1]
self.pic_type = pic_type
self.ax1 = self.fig.add_subplot(self.nrow,self.ncol,1,projection='windrose')
默认参数
# 坐标风向列表
self.dir = ['E', 'ENE', 'NE', 'NNE', 'N', 'NNW', 'NW', 'WNW', 'W', 'WSW', 'SW', 'SSW', 'S', 'SSE',
'SE', 'ESE']
# 风级标签列表
self.levels = ['≤1级', ' 2 级', ' 3 级', ' 4 级', ' 5 级', '>5级']
# 颜色映射列表
self.colors = ['#0000bc', '#3a93f1', '#40ffe7', '#cbfc7c', '#fdd305', '#f46c06']
# 用于图例设置,不设置可以删除
self.patches = [mpatches.Patch(color=self.colors[i], label=self.levels[i]) for i in range(len(self.colors))]
# 风级下限列表
self.bins = [0, 1.6, 3.4, 5.5, 8.0, 10.8]
# 预设存放表格数据
self.table_data = []
(2)绘图方法draw_windrose
【-】调整画布布局
# [left, bottom, width, height]
self.ax1.set_position([0.15, 0.15, 0.75, 0.75])
【-】设置坐标16方位
self.ax1.set_thetagrids(np.arange(0.0, 360.0, 22.5), labels=self.dir)
【-】图种选择设置+表格数据读取
if self.pic_type == 'bar':
self.ax1.bar(self.wd, self.ws, opening=0.8, colors=self.colors, bins=self.bins,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
elif self.pic_type == 'contour':
self.ax1.contour(self.wd, self.ws, bins=self.bins, colors=self.colors,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
elif self.pic_type == 'contourf':
self.ax1.contourf(self.wd, self.ws, bins=self.bins, colors=self.colors,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
【1】图种分析
柱状图(bar
)和等值线图(contour
或contourf
),这两个方法Windrose库自设了对角度的转换,所以输入的风向数据不需要进行弧度制的转换。
同时,这两种图种的设置,调用bar
或contour
,并不会像笛卡尔坐标系或者极坐标系那样,会返回一个对象。如下:
>>>ax = plt.subplot(111, projection='polar')
>>>CS1 = ax.contourf(theta,r,Z)
>>>CS1
<matplotlib.contour.QuadContourSet object at ······>
【2】表格读取(不需要的跳过)
Windrose
添加图种后是没有对象返回的,是直接将图像录入到了子图(ax
)里面。所以,表格数据的获取,就需要用到私有属性ax1._info['table']
(强烈建议。参考书籍:python气象资料处理与可视化[李开元 豆京华 著 ](只是高亮,实体书没链接))。
或者如果你介意、使用私有属性不符合编程习惯,会出现这样波浪线的提示警告:
可以读取子图里的patches
或patches_list
属性。
(1)ax1._info['table']
方法(强烈建议)
通过ax._info['table']
,会得到bins x dir
的一个各风向上风速的频次的二维数组,len(dir)
列,len(bins)
行。
(2)patches
或patches_list
方法
**表格数据读取:**不同图种使用不一样的填充图形对象。采用ax1.patches_list
方法,获取到:
bar
:对应matplotlib.patches.Rectangle
对象,patches_list
中包含6个元素。contour
:对应matplotlib.lines.Line2D
对象,patches_list
中包含6个元素。contourf
:对应matplotlib.patches.Polygon
对象,patches_list
中包含6个元素。
所以我们需要根据不同对象的属性获取其中各个风向上的风级量。而在柱状图(bar
)中,用patches_list
只能得到北风方向上的数据,要获取全部分量上的内容,需要使用ax.patches
方法。
各图种采用ax.patches
方法得到的内容:
bar
:<Axes.ArtistList of 96 patches>
,即16方位 × 风速分级数contour
:<Axes.ArtistList of 0 patches>
contourf
:<Axes.ArtistList of 6 patches>
,而contourf
在patches
和patches_list
的使用上并无区别。
所以bar
通过patches[i].get_height()
方法,循环遍历得到风向各风级列表,再重组成数组:
for i in range(len(self.ax1.patches)):
self.table_data.append(self.ax1.patches[i].get_height())
_table_data = np.array(self.table_data).reshape(16,6)
self.table_data = _table_data
contour
通过ax.patches_list[i].get_ydata()
方法,得到每级风的各风向对应的累积频次,注意只是累积频次。
每次获取到的列表共17个元素,第0个和-1个是相同方向,需要剔除最后一个。
for i in range(len(self.ax1.patches_list)):
self.table_data.append(self.ax1.patches_list[i].get_data()[1])
# 或者
# self.table_data.append(self.ax1.patches_list[i].get_ydata())
_table_data = np.array(self.table_data)
_table_data = np.delete(_table_data, -1, axis=1)
self.table_data = _table_data.T
self.table_data = column_transform(1, self.table_data)
>>>_draw.ax1.patches_list[0].get_ydata()
array([377., 531., 539., 237., 187., 258., 175., 155., 179., 204., 238.,226., 203., 203., 213., 305., 377.])
contourf
并没有get_ydata()
方法,通过ax.patches_list[i].xy
方法,得到(19, 2)
的xy
数组。
该数组倒数第二行为原点,倒数第一行为初始点,即等于第0行。
通过遍历获取前16行的y轴,得到每级风的各风向对应的累积频次,注意同样也只是累积频次。
for i in range(len(self.ax1.patches_list)):
self.table_data.append(self.ax1.patches_list[i].get_data()[1])
# 或者
# self.table_data.append(self.ax1.patches_list[i].get_ydata())
_table_data = np.array(self.table_data)
_table_data = np.delete(_table_data, -1, axis=1)
self.table_data = _table_data.T
self.table_data = column_transform(1, self.table_data)
>>>_draw.ax1.patches_list[0].xy.shape
(19, 2)
针对累积频次,需要除了第一列外所有列,都要减去前一列,得到的数才是每个风级对应频次。具体数组列变换的函数设计如下:
def column_transform(n,arr):
change_arr_list = []
if n == 0:
raise ValueError('起始列不能为0!')
# 计算差值列的最终值
for i in range(arr[:,n:].shape[1]):
change_arr = arr[:,n+i] - arr[:,n+i-1]
change_arr_list.append(change_arr)
# 赋予对应列
for j in range(arr[:,n:].shape[1]):
arr[:,n+j] = change_arr_list[j]
return arr
【-】设置每一圈的坐标
主要是将频次改成比例。
self.ax1.set_yticks([self.sum_freq * 0.05, self.sum_freq * 0.10, self.sum_freq * 0.15, self.sum_freq * 0.20])
self.ax1.set_yticklabels(['5%', '10%', '15%', '20%'])
【-】设置每一圈的坐标朝北
self.ax1.set_rlabel_position(90)
(3)色标设置
def set_colorbar(self):
# 设置颜色归一化,上设颜色划分为6,色标共6格,7个坐标,从0开始,6最大
norm = mpl.colors.Normalize(vmin=0, vmax=len(self.levels))
# 创建一个线性分段的颜色映射,使用自定义的颜色列表
cmap = mpl.colors.LinearSegmentedColormap.from_list('How2matplotlib_custom', self.colors, N=len(self.levels))
# 创建一个ScalarMappable对象,用于映射数据到颜色
sm = cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([]) # 防止显示数值
# 创建色标
cbar = plt.colorbar(sm, ax=self.ax1, orientation='horizontal')
cbar.ax.set_xticks(np.arange(0,len(self.levels)+1,1)) # 设置新的刻度位置
cbar.ax.set_xticklabels(['']+self.levels)
Windrose
库是不具备色标的设置的,如果想要设置色标,只能通过matplotlib
库生成自定义色标。
- 首先设置颜色归一化,添加划分的层次,如设置了6层风级,则是0到6。
- 颜色映射:
mpl.colors.LinearSegmentedColormap.from_list(风格名称,颜色列表,颜色分段数)
- ScalarMappable对象:映射数据到颜色
- 创建色标
不过一般不画色标,纯属顺手带上。
(4)图例设置
在windrose
库中是具备set_legend
方法,但是这个设置图例的方法,是不允许修改它的标签值的。它只有默认的形式:
根据其源代码:
def set_legend(self, **pyplot_arguments):
if "borderaxespad" not in pyplot_arguments:
pyplot_arguments["borderaxespad"] = -0.10
legend = self.legend(**pyplot_arguments)
plt.setp(legend.get_texts(), fontsize=8)
return legend
def legend(self, loc="lower left", decimal_places=1, units=None, **kwargs):
……
……
def get_labels(decimal_places=1, units=None):
digits = np.copy(self._info["bins"]).tolist()
if not digits:
return ""
digits[-1] = digits[-2]
digits = [f"{label:.{decimal_places}f}" for label in digits]
……
……
return labels
kwargs.pop("labels", None)
kwargs.pop("handles", None)
handles = get_handles()
labels = get_labels(decimal_places, units)
self.legend_ = mpl.legend.Legend(self, handles, labels, loc=loc, **kwargs)
return self.legend_
注意看它的get_labels
,它是直接copy
了self._info["bins"]
的数值,你即便传入的labels不是空值,它也不会将其传给函数get_labels
上,所以你的图例无论如何也动不了它的标签值。
只能采用matplotlib
的legend
方法,跟平常自定义图例方法一样。
# 颜色映射列表
self.colors = ['#0000bc', '#3a93f1', '#40ffe7', '#cbfc7c', '#fdd305', '#f46c06']
# 用于色标设置,不设置可以删除
self.patches = [mpatches.Patch(color=self.colors[i], label=self.levels[i]) for i in range(len(self.colors))]
def set_legend(self):
# Windroses设置图例方法
self.ax1.legend_ = legend.Legend(parent=self.ax1, handles=self.patches, labels=self.levels,
loc='lower left', bbox_to_anchor=(-0.15, -0.15), fontsize='small')
先将颜色列表做成patch
对象,作为handles`参数,其余的基本就是按部就班。
(5)添加图表
【-】频次表
def add_table1(self):
# 设置判定,只有子图元组某列等于2时,才会绘图
if self.nrow == 2 or self.ncol == 2:
self.ax2 = self.fig.add_subplot(self.nrow,self.ncol,2)
# 设置子图间间距
plt.subplots_adjust(wspace=0.2)
# 定义表行名,为各风向
rowLabel = self.dir[5:] + self.dir[:5]
# 由于是频次统计,所以只有整数
_table_data = self.table_data.astype(int)
# 将预设的风级标签并入数组中充当列名,之所以这么做是希望行列名的表格线统一
_table_data = np.vstack((self.levels,_table_data))
# 添加表格
table = self.ax2.table(cellText=_table_data, loc='center left', rowLoc='center',
rowLabels=[''] + rowLabel[::-1],edges='horizontal')
# 设置字体大小
table.set_fontsize(14)
# 表格放大1.2倍
table.scale(1.2, 1.2)
# 关闭轴线,取消刻度标签
self.ax2.spines['bottom'].set_visible(False)
self.ax2.spines['left'].set_visible(False)
self.ax2.spines['right'].set_visible(False)
self.ax2.spines['top'].set_visible(False)
self.ax2.get_xaxis().set_ticks([])
self.ax2.get_yaxis().set_ticks([])
【-】频率表
基本与频次表一致,只是在_table_data
上增加了除以总频次和所有元素转换为字符串。详见完整代码
_table_data = array_str_percentage(self.table_data / self.sum_freq)
(6)静风率文本
def add_calm_text(self):
# ax1.transAxes 转换成轴坐标。在轴坐标系中,坐标的原点(0,0)位于图表的左下角,而(1,1)则代表图表的右上角,无论图表的实际数据范围如何。
if self.C is None:
self.C = str(np.round(self.ax1.calm_count * 100 / self.sum_freq,3)) + '%'
self.ax1.text(0.5, -0.1, '静风频率:' + self.C,horizontalalignment='center', verticalalignment='center', transform=self.ax1.transAxes)
可选择手动计算静风率带入进去,主要方便传入的不是可以索引的数据。也可以使用windrose
本身算得的静风率。
(7)静风圈
实际上,调用windrose
绘图时,它本身就将静风的分布绘制进去,只是静风占比很小,所以看着并不明显。
而全网上有些中间放置静风率的圆圈该怎么画呢?
其实也是用它源代码里的方法:
将其内的圆圈半径放大,然后增加静风率文本:
def add_calm_circle(self):
self.ax1.set_rorigin(-120)
# 在圆圈中添加文本标注
if self.calm_limit is None:
self.ax1.text(0, -10, self.C,ha='right', va='center', fontsize=8, color='black')
else:
self.ax1.text(0, -10, f'{round(self.ax1.calm_count * 100 / self.sum_freq, 1)}%',
ha='right', va='center',fontsize=8, color='black')
四、完整代码
import os
import calendar
import matplotlib.cm as cm
import numpy as np
import pandas as pd
import datetime
from matplotlib import pyplot as plt
import matplotlib as mpl
from windrose import WindroseAxes
import matplotlib.patches as mpatches
import matplotlib.legend as legend
mpl.rcParams['font.sans-serif'] = ['SimHei'] #解决中文显示的问题
plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
plt.rcParams['axes.unicode_minus'] = False
class meteorological_date2datetime:
def __init__(self, extract_date: datetime.date):
self.date = extract_date
@property
def day_hour(self):
_day_hour = []
_year = self.date.year
_month = self.date.month
_day = self.date.day
_hours_list = [21, 22, 23] + [i for i in range(21)]
# 是否为昨天的判断标识
pre_day = 1
for i in _hours_list:
if pre_day == 1 and datetime.time(i,0) != datetime.time(0,0):
_datetime = datetime.datetime(_year,_month,_day,i,0) - datetime.timedelta(days=1)
_day_hour.append(_datetime)
else:
pre_day = 0
_datetime = datetime.datetime(_year, _month, _day, i, 0)
_day_hour.append(_datetime)
return _day_hour
@property
def day_minute(self):
_day_minute = []
_year = self.date.year
_month = self.date.month
_day = self.date.day
# 预组合时分
_hours_list = [21, 22, 23] + [i for i in range(21)]
_minute_list = [i+1 for i in range(59)] + [0]
# 是否为昨天的判断标识,是0,否1
pre_day = 1
for i in _hours_list:
for j in _minute_list:
# 先判断是否为昨天,循环判断是否经过00时,昨天今天日的分界线。
if datetime.time(i,j) == datetime.time(0,0):
pre_day = 0
_datetime = datetime.datetime(_year, _month, _day, i, 0)
elif pre_day == 0:
# 今天日的时分
if datetime.time(i,j).minute == 0:
_datetime = datetime.datetime(_year, _month, _day, i, 0)
else:
_datetime = datetime.datetime(_year, _month, _day, i, j) - datetime.timedelta(hours=1)
else:
# 昨天日的时分
if datetime.time(i,j).minute == 0:
_datetime = datetime.datetime(_year, _month, _day, i, 0) - datetime.timedelta(days=1)
else:
_datetime = datetime.datetime(_year, _month, _day, i, j) - datetime.timedelta(hours=1,days=1)
_day_minute.append(_datetime)
return _day_minute
def extract_data(_file_path,name):
df_list = []
date = []
for m in range(12):
year = 2023
month = str(m + 1) + '月'
data = pd.read_excel(_file_path, sheet_name=month, header=0, index_col=0)
day = data.columns.values
for i in day:
_date = datetime.date(year,m+1,i)
date.append(_date)
time_list = meteorological_date2datetime(_date).day_hour
day_data = data[i].replace(-999, float('nan')).values
day_df = pd.DataFrame(day_data, index=time_list)
df_list.append(day_df)
_df = pd.concat(df_list)
_df = _df.rename(columns=dict(zip(_df.columns, [name])))
return _df
# 将数组转换为百分制
def array_str_percentage(arr):
_values = []
_shape = arr.shape
_n = len(_shape)
for i in range(_shape[0]):
value = arr[i]
# 第一层循环,判断出目前是一维,传入的就是二维
if len(value.shape) == 1:
_values.append(['{:.2f}%'.format(a * 100) for a in value])
else:
_values.append(array_str_percentage(arr[i]))
return np.array(_values)
class wind_subplots:
def __init__(self,
_wd, _ws, _sum_freq,
_C:str=None,
_calm_limit:float=None,
sub_tup:tuple=(1,1),
pic_type:str=None):
"""
:param _wd: 风向数据
:param _ws: 风速数据
:param _sum_freq: 总频次
:param _C: 自算静风率
:param _calm_limit:静风限制,windrose会根据大于该值计算静风频次
:param sub_tup: 子图元组(1,1)或(1,2),默认图1为风玫瑰图,图2为图表(选设)
:param pic_type: 图1选择绘制图种样式
"""
# 基础参数
self.wd = _wd
self.ws = _ws
self.sum_freq = _sum_freq
self.C = _C
self.calm_limit = _calm_limit
if self.calm_limit:
if self.calm_limit <= 0.0:
self.calm_limit = None
if isinstance(self.wd,pd.Series) or isinstance(self.wd, pd.DataFrame) or isinstance(self.wd,np.ndarray):
pass
elif isinstance(self.ws,pd.Series) or isinstance(self.ws, pd.DataFrame) or isinstance(wd,np.ndarray):
pass
else:
raise TypeError("如果使用windrose自行判断静风,传入的数据必须是数组。仅支持pd.Series,pd.DataFrame,np.ndarray。")
if self.calm_limit is None:
raise ValueError(
"错误的静风限制值,如果使用windrose自行判断静风,calm_limit必须大于0。")
else:
self.calm_limit = None
# 子图设置
if sub_tup[1] == 2:
self.fig = plt.figure(figsize=(12,6))
self.ax2 = None
else:
self.fig = plt.figure(figsize=(6, 6))
self.nrow = sub_tup[0]
self.ncol = sub_tup[1]
self.pic_type = pic_type
self.ax1 = self.fig.add_subplot(self.nrow,self.ncol,1,projection='windrose')
# 默认参数
self.dir = ['E', 'ENE', 'NE', 'NNE', 'N', 'NNW', 'NW', 'WNW', 'W', 'WSW', 'SW', 'SSW', 'S', 'SSE',
'SE', 'ESE']
self.levels = [' 1 级', ' 2 级', ' 3 级', ' 4 级', ' 5 级', '>5级']
self.colors = ['#0000bc', '#3a93f1', '#40ffe7', '#cbfc7c', '#fdd305', '#f46c06']
self.patches = [mpatches.Patch(color=self.colors[i], label=self.levels[i]) for i in range(len(self.colors))]
self.bins = [0, 1.6, 3.4, 5.5, 8.0, 10.8]
# 预设存放表格数据
self.table_data = []
def draw_windrose(self):
# 设置风玫瑰图偏移画布中心,多出空白放置图例
self.ax1.set_position([0.15, 0.15, 0.75, 0.75]) # [left, bottom, width, height]
# 设置成16方位
self.ax1.set_thetagrids(np.arange(0.0, 360.0, 22.5), labels=self.dir)
if self.pic_type == 'bar':
self.ax1.bar(self.wd, self.ws, opening=0.8, colors=self.colors, bins=self.bins,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
elif self.pic_type == 'contour':
self.ax1.contour(self.wd, self.ws, bins=self.bins, colors=self.colors,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
elif self.pic_type == 'contourf':
self.ax1.contourf(self.wd, self.ws, bins=self.bins, colors=self.colors,
calm_limit=self.calm_limit)
# 读取表格数据,只画一图可删除
self.table_data = self.ax1._info['table'].T
# 换算每一圈的比例
self.ax1.set_yticks([self.sum_freq * 0.05, self.sum_freq * 0.10, self.sum_freq * 0.15, self.sum_freq * 0.20])
self.ax1.set_yticklabels(['5%', '10%', '15%', '20%'])
# 因为采用原理为极坐标,在极坐标上,90°朝上。因此我们希望将每一圈的比例坐标轴指向北,则选择90°方位
self.ax1.set_rlabel_position(90)
def set_colorbar(self):
# 上设颜色划分为6,色标共6格,7个坐标,从0开始,6最大
norm = mpl.colors.Normalize(vmin=0, vmax=len(self.levels))
cmap = mpl.colors.LinearSegmentedColormap.from_list('How2matplotlib_custom', self.colors, N=len(self.levels))
sm = cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([]) # 防止显示数值
# 创建色标
cbar = plt.colorbar(sm, ax=self.ax1, orientation='horizontal')
cbar.ax.set_xticks(np.arange(0,len(self.levels)+1,1)) # 设置新的刻度位置
cbar.ax.set_xticklabels(['']+self.levels)
def set_legend(self):
# Windroses设置图例方法
self.ax1.legend_ = legend.Legend(parent=self.ax1, handles=self.patches, labels=self.levels,
loc='lower left', bbox_to_anchor=(-0.15, -0.15), fontsize='small')
# 表格显示频率
def add_table(self):
if self.nrow == 2 or self.ncol == 2:
self.ax2 = self.fig.add_subplot(self.nrow,self.ncol,2)
plt.subplots_adjust(wspace=0.2)
rowLabel = self.dir[5:] + self.dir[:5]
_table_data = array_str_percentage(self.table_data / self.sum_freq)
_table_data = np.vstack((self.levels,_table_data))
table = self.ax2.table(cellText=_table_data, loc='center left', rowLoc='center',
rowLabels=[''] + rowLabel[::-1],edges='horizontal')
table.set_fontsize(14)
table.scale(1.2, 1.2)
self.ax2.spines['bottom'].set_visible(False)
self.ax2.spines['left'].set_visible(False)
self.ax2.spines['right'].set_visible(False)
self.ax2.spines['top'].set_visible(False)
self.ax2.get_xaxis().set_ticks([])
self.ax2.get_yaxis().set_ticks([])
# 表格显示频次
def add_table1(self):
if self.nrow == 2 or self.ncol == 2:
self.ax2 = self.fig.add_subplot(self.nrow,self.ncol,2)
plt.subplots_adjust(wspace=0.2)
rowLabel = self.dir[5:] + self.dir[:5]
_table_data = self.table_data.astype(int)
_table_data = np.vstack((self.levels,_table_data))
table = self.ax2.table(cellText=_table_data, loc='center left', rowLoc='center',
rowLabels=[''] + rowLabel[::-1],edges='horizontal')
table.set_fontsize(14)
table.scale(1.2, 1.2)
self.ax2.spines['bottom'].set_visible(False)
self.ax2.spines['left'].set_visible(False)
self.ax2.spines['right'].set_visible(False)
self.ax2.spines['top'].set_visible(False)
self.ax2.get_xaxis().set_ticks([])
self.ax2.get_yaxis().set_ticks([])
def add_calm_circle(self):
self.ax1.set_rorigin(-120)
# 在圆圈中添加文本标注
if self.calm_limit is None:
self.ax1.text(0, -10, self.C,
ha='right', va='center', fontsize=8, color='black')
else:
self.ax1.text(0, -10, f'{round(self.ax1.calm_count * 100 / self.sum_freq, 1)}%',
ha='right', va='center',fontsize=8, color='black')
def add_calm_text(self):
# ax1.transAxes 转换成轴坐标。在轴坐标系中,坐标的原点(0,0)位于图表的左下角,而(1,1)则代表图表的右上角,无论图表的实际数据范围如何。
if self.C is None:
self.C = str(np.round(self.ax1.calm_count * 100 / self.sum_freq,3)) + '%'
self.ax1.text(0.5, -0.1, '静风频率:' + self.C,
horizontalalignment='center', verticalalignment='center', transform=self.ax1.transAxes)
if __name__ == '__main__':
dir_name = r'.\2winddir旧.xlsx'
spd_name = r'.\2windspd旧.xlsx'
dir_df = extract_data(dir_name,'风向')
spd_df = extract_data(spd_name,'风速')
df = pd.concat([dir_df, spd_df],axis=1).dropna()
# 计算静风率
df_C = df[df['风向'] == 'PPC']
C = np.round(len(df_C.index)/len(df.index), 3)
C = str(np.round(C * 100, 1)) + '%'
# 风向风速
wd = df['风向']
ws = df['风速']
# df = df[df['风向'] != 'PPC']
# wd = list(df['风向'])
# ws = list(df['风速'])
pic = 'bar'
_draw = wind_subplots(wd,ws,len(df.index),_C=C,_calm_limit=0.2,sub_tup=(1,1),pic_type=pic)
_draw.draw_windrose()
# _draw.set_colorbar()
_draw.set_legend()
_draw.add_table()
# _draw.add_calm_text()
_draw.add_calm_circle()
# plt.show()
plt.savefig('{}.png'.format(pic))
plt.close()
# plt.savefig('table-{}静风法.png'.format(pic))
# plt.close()