使用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 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 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()