使用ERA5数据绘制风向玫瑰图的简易流程

使用ERA5数据绘制风向玫瑰图的简易流程

今天需要做一个2017年-2023年的平均风向的统计,做一个风向玫瑰图，想到的还是高分辨率的ERA5land的数据（0.1°分辨率，逐小时分辨率，1950年至今）。

风向，我分为了16个（0-360°，北方向为0），统计该时间段内的16个风向频率。

下载

使用Google earth engine快速统计风向频率：

var ROI = 你的区域;
var startDate = '2023-1-01';
var endDate = '2023-01-30';var dataset = ee.ImageCollection('ECMWF/ERA5_LAND/HOURLY').select(['u_component_of_wind_10m', 'v_component_of_wind_10m']).filterDate(startDate, endDate).filter(ee.Filter.calendarRange(11, 4, 'month'));var calculateWindDirection = function(image) {var direction = image.select('u_component_of_wind_10m', 'v_component_of_wind_10m').expression('atan2(v, u) * 180 / PI + 180',{'u': image.select('u_component_of_wind_10m'),'v': image.select('v_component_of_wind_10m'),'PI': Math.PI});return direction.rename('wind_direction');
};// 计算16个方向的频率
var directions = ee.List.sequence(0, 15);
var binSize = 360/16;var directionalFrequency = directions.map(function(dir) {var start = ee.Number(dir).multiply(binSize);var end = start.add(binSize);var directionMask = dataset.map(calculateWindDirection).map(function(img) {return img.gte(start).and(img.lt(end));});    return directionMask.mean().rename(ee.String('dir_').cat(ee.Number(dir).format('%02d')));
});// 将List转换为Image Collection，然后合并为一个多波段图像
var allDirections = ee.ImageCollection.fromImages(directionalFrequency).toBands();// 修改波段名称
var newBandNames = directions.map(function(dir) {return ee.String('dir_').cat(ee.Number(dir).format('%02d'));
}).getInfo();// 重命名波段
allDirections = allDirections.rename(newBandNames);// 导出数据
Export.image.toDrive({image: allDirections,description: 'Wind_Direction_Frequency_16dirs',scale: 10000,region: ROI,fileFormat: 'GeoTIFF',maxPixels: 1e9
});

下载下来后，放到qgis里面看看，一共16个波段，每个波段都代表着一个方向上的频率，16个波段加起来是1：

制图

使用python3实现：

import numpy as np
import matplotlib.pyplot as pltdef plot_wind_rose(data, title='Wind Rose'):"""绘制风向玫瑰图data: 包含16个方向频率的数组"""# 创建图形fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': 'polar'})# 设置方向角度（16个方向，每个22.5度）angles = np.arange(0, 360, 22.5) * np.pi/180# 确保数据是闭合的（首尾相连）frequencies = np.append(data, data[0])angles = np.append(angles, angles[0])# 绘制极坐标图ax.plot(angles, frequencies, 'o-', linewidth=2, color='purple')ax.fill(angles, frequencies, alpha=0.25, color='purple')# 设置方向标签ax.set_xticks(angles[:-1])direction_labels = ['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE','S', 'SSW', 'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW']ax.set_xticklabels(direction_labels)# 设置网格和刻度ax.grid(True)# 设置频率刻度范围max_freq = np.max(frequencies)ax.set_ylim(0, max_freq * 1.1)# 设置标题ax.set_title(title)return figdef read_wind_data(tiff_path, x, y):"""读取特定位置的风向数据"""with rasterio.open(tiff_path) as src:# 将经纬度转换为像素坐标row, col = src.index(x, y)# 读取所有波段在该位置的值data = []for i in range(1, src.count + 1):value = src.read(i)[row, col]data.append(float(value))return np.array(data)# 使用示例
import rasterio
#输入tif路径
tiff_path = r'\风向数据\Wind_Direction_Frequency_16dirs.tif'
#输入经纬度
x, y = 99, 25.312# 读取数据
wind_data = read_wind_data(tiff_path, x, y)# 打印数据检查
print("Wind direction frequencies:")
for i, freq in enumerate(wind_data):print(f"Direction {i*22.5:>6.1f}°: {freq:>6.3f}")# 绘制风向玫瑰图
fig = plot_wind_rose(wind_data, f'Wind Rose at ({x}, {y})')# 保存图片
plt.savefig('wind_rose.png', dpi=300, bbox_inches='tight')
plt.show()

在代码中填入需要生成的风玫瑰图的经纬度，即可获得2017-2023年的该地区风向情况。

参考

Claude sonnet 3.5