风玫瑰图整合(一)风速玫瑰图

风玫瑰图整合(一)风速玫瑰图

一、前言

全网中关于python绘制风玫瑰图的参考材料较少,基于第一期的风玫瑰图,我们再做一次拓展与整合解析

基于气象A文件绘制风玫瑰图

首先风玫瑰图分为:风向玫瑰图 和 风速玫瑰图。在多数论文中,提及的风玫瑰图一般指风向玫瑰图。其中风速玫瑰图常见有:
风速玫瑰图
风向玫瑰图常见有:
风向玫瑰图

本篇基于第一篇风玫瑰图的图形,主要解析风速玫瑰图。

二、风速玫瑰图怎么看

在这里插入图片描述

这是一张风速玫瑰图(a),采用的绘制方法对应极坐标系下的柱状图bar。而右图为,基于左图数据得到的无静风的各风向对应风速的频率。针对绘图的数据,其设置的分割区间为:

0.3—1.51.6—3.33.4—5.55.6—7.98.0—10.710.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.Seriespd.DataFramenp.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)和等值线图(contourcontourf),这两个方法Windrose库自设了对角度的转换,所以输入的风向数据不需要进行弧度制的转换。

同时,这两种图种的设置,调用barcontour,并不会像笛卡尔坐标系或者极坐标系那样,会返回一个对象。如下:

>>>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气象资料处理与可视化[李开元 豆京华 著 ](只是高亮,实体书没链接))。

或者如果你介意、使用私有属性不符合编程习惯,会出现这样波浪线的提示警告:
在这里插入图片描述

可以读取子图里的patchespatches_list属性。

(1)ax1._info['table']方法(强烈建议)

通过ax._info['table'],会得到bins x dir的一个各风向上风速的频次的二维数组,len(dir)列,len(bins)行。

(2)patchespatches_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>,而 contourfpatchespatches_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,它是直接copyself._info["bins"]的数值,你即便传入的labels不是空值,它也不会将其传给函数get_labels上,所以你的图例无论如何也动不了它的标签值。

只能采用matplotliblegend方法,跟平常自定义图例方法一样。

# 颜色映射列表
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()

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值