使用ERA5数据绘制风向玫瑰图的简易流程
今天我需要做一个2017年-2023年的
平均风向的统计
,生成一个风向玫瑰图。
还是一如既往地选择ERA5land的数据（0.1°分辨率，逐小时分辨率，1950年至今），当然你也可以用ERA5数据，我只是觉得分辨率比较低就没用。
做的时候，我把风向分为了16个（0-360°，北方向为0），即统计该时间段内的16个风向频率，然后生成tif，基于tif统计数据获取指定区域的任意风向玫瑰图。
下载
使用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 plt
def 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 fig
def read_wind_data(tiff_path, x, y):
"""
读取特定位置的风向数据
"""
with rasterio.open(tiff_path) as
