ERA5小时压力层级数据处理的Python实践（进阶篇）

一、需求分析

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

四、结果

作者：dairui1240