基于python和era再分析风场的的风玫瑰图制作

1. 简介

当今风电新能源规划如火如荼，对于风电建设，确定规划区域的风玫瑰图是基础的步骤，对于缺少实测资料的地区，可以通过再分析数据获取该区域的风场历时数据，并基于此绘制风玫瑰图。
近年来，era5的再分析气象数据得到了越来越广泛的认可，本文基于python和era5数据绘制某区域的风玫瑰图。

2. 数据源

我们下载了南海区域2022年的风场数据，时间间隔为6小时，数据范围如下图所示。

3. 源代码

#导入需要的模块
import matplotlib.pyplot as plt
import numpy as np
from netCDF4 import Dataset
import datetime as dt
import windrose as wd
from windrose import WindroseAxes, WindAxes, plot_windrose
import matplotlib.cm as cm
plt.rcParams['font.sans-serif']=['Times New Roman']
#读取nc文件
file = Dataset('D:\软件\浏览器下载的文件\\2022nanhaiwind.nc')
print(file)
lon = file['longitude'][:]
lat = file['latitude'][:]
u10 = file['u10'][:]
v10 = file['v10'][:]
time =file['time'][:]
#给定测点的位置
xlon= 114
xlat = 22.51#运用双线性插值计算测点的u、v风速时间序列
num_lon = len(lon)
num_lat = len(lat)
for i in range(num_lon-1):
        if xlon<lon[i+1] and lon[i]<=xlon:
                gnumx=i
                break
for i in range(num_lat-1):
        if xlat>lat[i+1] and lat[i]>=xlat:
                gnumy=i
                break
s=0.25*0.25
a1 = np.abs((xlat-lat[gnumy])*(xlon-lon[gnumx]))/s
a2 = np.abs((xlat-lat[gnumy+1])*(xlon-lon[gnumx]))/s
a3 = np.abs((xlat-lat[gnumy])*(xlon-lon[gnumx+1]))/s
a4 = np.abs((xlat-lat[gnumy+1])*(xlon-lon[gnumx+1]))/s
print(a1,a2,a3,a4)
u_=[]
v_=[]
for i in range(len(time)):
     	u_.append(a1*u10[i,gnumy+1,gnumx+1]+a2*u10[i,gnumy,gnumx+1]+a3*u10[i,gnumy+1,gnumx]+a4*u10[i,gnumy,gnumx])
        v_.append(a1 * v10[i, gnumy + 1, gnumx + 1] + a2 * v10[i, gnumy, gnumx + 1] + a3 * v10[i, gnumy + 1, gnumx] + a4 *v10[i, gnumy, gnumx])
u_ = np.array(u_)
v_ = np.array(v_)
#计算和风速
vel = np.sqrt(u_**2+v_**2)  
#计算风向
dir = []
for i in range(len(u_)):
                if u_[i]>0 and v_[i]>0:
                dir.append(1.5*np.pi-np.arctan(abs(v_[i]/u_[i])))
        elif u_[i]<0 and v_[i]>0:
                dir.append(0.5*np.pi+np.arctan(abs(v_[i]/u_[i])))
        elif u_[i] < 0 and v_[i] < 0:
                dir.append(0.5*np.pi - np.arctan(abs(v_[i] / u_[i])))
        elif u_[i] > 0 and v_[i] < 0:
                dir.append(1.5*np.pi + np.arctan(abs(v_[i] / u_[i])))
#风向转换为角度单位
dir =np.array(dir)/np.pi*180
font1 = {'family': 'Times New Roman',
         'weight': 'normal',
         'size': 20,
         }
#绘制风玫瑰图程序
ax = WindroseAxes.from_ax(figsize=(16, 10))  #自北向南为0度，顺时针旋
ax.bar(dir, vel, bins=np.arange(0.,16,2), cmap=cm.rainbow, lw=5)
ax.set_legend(fontsize=20)
plt.legend(loc=(2.09,0.6),bbox_to_anchor=(1.07,0.015),fontsize=20)
ax.tick_params(labelsize=20)
ax.tick_params(labelleft=False)
plt.show()

4. 绘制效果

我们以深圳为测点，绘制效果如下。

与网上找到的深圳风玫瑰图（来源不太明确）做一下比较，可知结果整体上是合理的。