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. 绘制效果
我们以深圳为测点,绘制效果如下。
与网上找到的深圳风玫瑰图(来源不太明确)做一下比较,可知结果整体上是合理的。