个人备案域名可以做企业网站吗,登录百度账号,电商 网站建设,苏州哪家网站公司做的好的使用ERA5数据绘制风向玫瑰图的简易流程
今天需要做一个2017年-2023年的平均风向的统计,做一个风向玫瑰图#xff0c;想到的还是高分辨率的ERA5land的数据#xff08;0.1分辨率#xff0c;逐小时分辨率#xff0c;1950年至今#xff09;。
风向#xff0c;我分为了16个想到的还是高分辨率的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, titleWind 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-, linewidth2, colorpurple)ax.fill(angles, frequencies, alpha0.25, colorpurple)# 设置方向标签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(fDirection {i*22.5:6.1f}°: {freq:6.3f})# 绘制风向玫瑰图
fig plot_wind_rose(wind_data, fWind Rose at ({x}, {y}))# 保存图片
plt.savefig(wind_rose.png, dpi300, bbox_inchestight)
plt.show()在代码中填入需要生成的风玫瑰图的经纬度即可获得2017-2023年的该地区风向情况。 参考
Claude sonnet 3.5
