一、需求分析
1.1 ERA5数据的U、V风合成计算风向、风速。
#在(一)的基础上加入风向、风速计算
二、程序设计
2.1(引自原文链接:u, v风和风速风向的相互转换_uv分量的风怎么换算成风速-CSDN博客)
博主“气海同途”的文章中详细介绍了u, v风和风速风向的相互转换。
气象上定义正北方向为0°(即风从北吹向南,是数学坐标系中的 -90°), 顺时针转动角度增大。
风向Dir=0°(或360°), u=0, v<0, 正北风;
风向Dir=90°, u<0, v=0,正东风;
风向Dir=180°, u=0, v>0,正南风;
风向Dir=270°, u>0, v=0,正西风。
三角函数所使用的极坐标系,其0°对应X轴正方向,逆时针为正;而气象学中0°对应的是Y轴正方向,顺时针为正。
这里需要通过u, v风计算风速wspd和wdir:
1.风速计算:
风速是U和V分量的矢量和的模。这里使用np.sqrt(U**2 + V**2)
来计算风速。
2.风向计算:
风向是U和V分量合成后的方向,结果以度为单位。
deg = 180.0/np.pi #(角度制与弧度制的相互转换)
wdir1 = 180.0 + np.arctan2(u, v)*deg
wdir2 = 270.0 - np.arctan2(v, u)*deg
#上面两行代码均可将计算出的风向转到0-360°之间。
三、完整代码
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 8 15:47:26 2024
@author: DR
"""
import netCDF4 as nc
file_path = "D:/test/202405.nc"
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
# 加载ERA5数据
file_path = 'D:/test/202405.nc' # 替换为你的ERA5数据文件路径
dataset = xr.open_dataset(file_path)
# 筛选时间点
timestamp = np.datetime64('2024-05-03T14:00:00')
dataset_time_filtered = dataset.sel(time=timestamp, method='nearest')
# 筛选地点
latitude1, longitude1 = 31.8667, 117.2833 # 北纬31度52分、东经117度17分
latitude2, longitude2 = 31.0, 117.3 # 北纬31度、东经117.3度
# 提取两个地点的温度廓线、U-component of wind和V-component of wind
temperature_profile1 = dataset_time_filtered.sel(latitude=latitude1, longitude=longitude1, method='nearest').t
temperature_profile2 = dataset_time_filtered.sel(latitude=latitude2, longitude=longitude2, method='nearest').t
u_wind_profile1 = dataset_time_filtered.sel(latitude=latitude1, longitude=longitude1, method='nearest').u
u_wind_profile2 = dataset_time_filtered.sel(latitude=latitude2, longitude=longitude2, method='nearest').u
v_wind_profile1 = dataset_time_filtered.sel(latitude=latitude1, longitude=longitude1, method='nearest').v
v_wind_profile2 = dataset_time_filtered.sel(latitude=latitude2, longitude=longitude2, method='nearest').v
#计算风速、风向
deg = 180.0/np.pi
wind_speed1 = np.sqrt(u_wind_profile1**2 + v_wind_profile1**2)
wind_direction1 = 180.0 + np.arctan2(v_wind_profile1, u_wind_profile1)*deg
wind_speed2 = np.sqrt(u_wind_profile2**2 + v_wind_profile2**2)
wind_direction2 = 180.0 + np.arctan2(v_wind_profile2, u_wind_profile2)*deg
# 创建一个图表和两个子图
fig, (plt1, plt2, plt3) = plt.subplots(1, 3, figsize=(14, 6))
# 温度廓线
temperature_profile1.plot(y='level', marker='o', label=f'Temperature 31.8667N, 117.2833E', ax=plt1)
temperature_profile2.plot(y='level', marker='x', label=f'Temperature 31.0N, 117.3E', ax=plt1)
plt1.set_title('Temperature Profiles')
plt1.set_xlabel('Temperature (K)')
plt1.set_ylabel('Pressure Level (hPa)')
plt1.invert_yaxis() # 压力随高度增加而减小,因此反转y轴
plt1.legend()
plt1.grid(True)
# 注意:风速可能需要转换为其他单位,这里假设单位是m/s
wind_speed1.plot(y='level', marker='o', linestyle='--', label=f'wind_speed 31.8667N, 117.2833E', ax=plt2)
wind_speed2.plot(y='level', marker='x', linestyle='--', label=f'wind_speed 31.0N, 117.3E', ax=plt2)
plt2.set_title('Wind Components')
plt2.set_xlabel('Wind Speed (m/s)')
plt2.set_ylabel('Pressure Level (hPa)')
plt2.invert_yaxis() # 压力随高度增加而减小,因此反转y轴
plt2.legend()
plt2.grid(True)
wind_direction1.plot(y='level', marker='o', linestyle='-.', label=f'wind_direction 31.8667N, 117.2833E', ax=plt3)
wind_direction2.plot(y='level', marker='x', linestyle='-.', label=f'wind_direction 31.0N, 117.3E', ax=plt3)
plt3.set_title('Wind2 Components')
plt3.set_xlabel('wind_direction (°)')
plt3.set_ylabel('Pressure Level (hPa)')
plt3.invert_yaxis() # 压力随高度增加而减小,因此反转y轴
plt3.legend(loc='upper center')#将图例放置在顶部
plt3.grid(True)
# 调整布局并显示图表
plt.tight_layout()
plt.show()
四、结果