登录/ 注册
 找回密码
 立即注册
关闭

CSDN热搜

程序园 精品问答 技术交流 资源下载
返回列表 发新帖
首页 › 业界区 › 业界 ›Python：风坡夹角/风效因子的计算

Python：风坡夹角/风效因子的计算

户烫擞 7 小时前
01 风坡夹角的定义

风向与坡向夹角的余弦值和坡度正弦值的乘积.
02 说明

计算风坡夹角需要使用到ERA5-Land再分析资料的u10(横向上风的分量)和v10(纵向上风的分量)(ps: 计算风速wind使用勾股定理u10 ^ 2 + v10 ^ 2再对结果开方即可)，坡向和坡度可以通过DEM进行计算(利用ArcGIS/QGIS等).
2.1 u10和v10的定义

u10全称为10 metre U wind component，v10全称是10 metre V wind component。
ERA5-Land官方文档的定义是:
1.png

2.png

附官方文档链接:
u10，v10
从上面的描述中可以知道：u10实际上就是东西向的风分量朝东移动的速度，v10是南北向的风分量朝北移动的速度。
而根据气象学上关于风的定义，例如东风表示风从东边吹过来，所以是朝正西方向移动的风，某种程度上可以说这个定义与上述的u10和v10的定义刚好相反，前者是从哪边吹过来，后者是朝哪边吹过去。因此后续进行计算时尤其需要注意正负号的关系。
关于风向，应该计算为0-360范围(类似地，东风的风向为90°)，具体如下：
3.png

2.2 u10和v10的下载

下载：ERA5-Land每小时产品下载
2.2.1 选择下载的变量

4.png

2.2.2 选择下载的时间

5.png

2.2.3 选择地理范围和下载格式

6.png

2.2.4 下载

7.png

8.png

2.3 DEM的下载

对于DEM的下载已经坡向和坡度的计算较为简单，此处略
03计算

3.1 风向的计算

关于风行的定义前文已经提及，这里着重说明计算的逻辑。
由于气象学上风向的逻辑与ERA5-Land的u10和v10的逻辑正好相反，譬如：
风u10v10东风 10m/s-10 m/s0 m/s北风 5 m/s0 m/s-5 m/s西风 8 m/s8 m/s0 m/s正西南风 2m/s$\sqrt{2}$ m/s$\sqrt{2}$ m/s正东南风 6 m/s$-3\sqrt{2}$ m/s$3\sqrt{2}$ m/s其他风就不一一列举。
因此为了方便数学上关于向量的计算，提前将u10和v10添上负号，其负号的两个横纵变量在笛卡尔xy坐标系上的合分量的指向就是正确的风向了。
也就是说，利用arctan($\frac{-v}{-u}$)反算夹角即可得到风向。
代码如下：
  1. import numpy as np
  3. def cal_wind_direction(u, v):
  4.     # 计算来向向量 (-u, -v) 的数学角度 -- (来向向量: 风吹来的方向, 例如北风是从北边吹来的,所以来向向量的方向是正北)
  5.     wind_dir_rad = np.arctan2(-v, -u)  # 参数顺序为 y=-v, x=-u
  6.     """
  7.     这里是对于np.arctan2(y, x)是传入一个向量(x, y), 计算其与向量(1, 0)(该向量与X轴正方向一致,此处坐标均为数学坐标系/二维笛卡尔坐标系)的夹角,
  8.     夹角范围为(-180, 180), 若为一二象限即是正,若为三四象限即为负,这与数学上的逆时针为正角顺时针为负角一致.
  9.     此外,
  10.     此处的u和v,首先解释u表示风吹向东边的速度;v表示风吹向北边的速度;(负号表示与定义方向相反)
  11.     由上面可以知道, u表示x轴方向上的分速度, v表示y轴方向上的分速度;
  12.     基于u和v可以计算出风运动的方向和速度,速度这里我们暂时不管,但是对于风运动的方向,显然与向量(u, v)的方向一致;
  13.     但是风运动的方向和向量(u, v)的方向都是表示风往哪个方向吹或者风要吹到哪里去;这与风向的定义正好相反,例如:
  14.     北风表示从北边吹来的风,所以是吹往南边的风,假定风速是1m/s,那么u=0,v=-1m/s(由v的定义知往南吹与定义方向相反所以添上负号).
  15.     如果直接传入arctan2(-1, 0),那么实际上得到的是向量(0, -1)的方向与向量(1, 0)方向上的夹角.
  16.     但是我们应该是需要得到向量(0, 1)的方向与向量(1, 0)方向上的夹角.
  17.     大家自行体会为什么需要添上负号,本质就是u和v定义正方向是风吹向的方向,而风向定义的方向是风来时或者从哪里吹来的方向--方向刚好相反
  18.     """
  20.     # 将弧度转换为度数 [-pi, pi] ==> [-180, 180]
  21.     wind_dir_deg = np.degrees(wind_dir_rad)
  23.     # 转换为地理角度(正北为0°，顺时针旋转)
  24.     wind_dir_geo = (90 - wind_dir_deg) % 360
  25.     """
  26.     根据arctan2的定义知,其输出是指风向向量(当你输入是-u, -v时符合风向定义)与向量(1, 0)方向或者说x轴正方向之间的夹角
  27.     但是风向的定义北风为0°,如此顺时针旋转到360度回到正北.两个夹角的定义不相同, 需要进行换算
  28.     """
  30.     return wind_dir_geo
复制代码
ps：函数内有两个需要注意的就是np.arctan2的输入和输出，由于我们计算出来的夹角是需要在0-360°，普通的np.arctan函数只能计算一个周期Π即范围是[$\frac{-Π}{2}$, $\frac{Π}{2}$]，因此需要使用np.arctan2：

  • 其输入是先y(对应-v10)后x(对应-u10),注意参数的顺序
  • 其输出是[-Π，Π]，夹角(即为下方的src_angle)表示为向量(x, y)(即(-u10, -v10))与x轴正方向上的(1, 0)之间的夹角，按照数学逻辑逆时针为正（即在x轴上方的向量为正），顺时针为负（x轴下方的向量为负），换算为气象学上的风向为dst_angle = (90 - src_angle) % 360，其中dst_angle为气象学的风向角度值(0-360,0表示北风)。注意此处的%不同编程语言的运算逻辑在负数的取余上存在区别，因此注意使用其他语言计算时的区别（目前仅适用于python）。
我们尝试计算一些示例（风向完全依据cal_wind_direction函数输入下方u10和v10计算得到）：
风u10v10风向(°)东风 10m/s-10 m/s0 m/s90北风 5 m/s0 m/s-5 m/s0西风 8 m/s8 m/s0 m/s270正西南风 2m/s$\sqrt{2}$ m/s$\sqrt{2}$ m/s225正东南风 6 m/s$-3\sqrt{2}$ m/s$3\sqrt{2}$ m/s1353.2 风向和坡向夹角的计算

坡向的计算通过GIS软件例如ArcGIS即可实现，此处不再赘述.
风向和坡向的夹角范围为0-180，因此计算为：
  1. # 计算风向和坡向的夹角
  2. theta_diff = np.abs(wind_dir - aspect)
  3. theta = np.minimum(theta_diff, 360 - theta_diff)
复制代码
3.3 计算风坡夹角/风效因子
  1. # 计算风坡夹角
  2. ws_angle = np.cos(np.radians(theta)) * np.sin(np.radians(slope))
复制代码
这是具体函数(函数内部的cal_wind_direction为自定义函数,具体见3.1内容)：
  1. import rasterio as rio
  2. import numpy as np
  4. def cal_wind_slope_angle(out_path, u10_path, v10_path, aspect_path, slope_path):
  5.     """
  6.     基于u10和v10计算风向, 基于风向和坡向、坡度计算风坡夹角
  7.     风坡夹角: 风向与坡向夹角的余弦值和坡度正弦值的乘积
  8.     :param out_path: 风坡夹角的输出路径
  9.     :param u10_path: u10的输入路径
  10.     :param v10_path: v10的输入路径
  11.     :param aspect_path: 坡向的输入路径
  12.     :param slope_path: 坡度的输入路径
  13.     :return: None
  14.     """
  16.     # 读取经纬向风速,坡向坡度
  17.     with rio.open(u10_path) as f:
  18.         u10 = f.read(1, masked=True)  # 读取第一个波段
  19.         meta = f.meta
  20.         """
  21.         .meta返回当前tiff文件的元数据, 包括格格式(GTiff)、数据类型(dtype)、无效值(nodata)、行列数和波段数(width,height,count),
  22.         坐标参考系(crs)、仿射参数(transform)
  23.         示例:
  24.         {'driver': 'GTiff',
  25.          'dtype': 'float32',
  26.          'nodata': nan,
  27.          'width': 129,
  28.          'height': 133,
  29.          'count': 1,
  30.          'crs': CRS.from_wkt('GEOGCS["unknown",DATUM["unknown",SPHEROID["unknown",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]'),
  31.          'transform': Affine(0.1, 0.0, 97.30000000000001,
  32.                 0.0, -0.1, 34.4)}
  33.         """
  34.     with rio.open(v10_path) as f:
  35.         v10 = f.read(1, masked=True)
  36.     with rio.open(aspect_path) as f:
  37.         aspect = f.read(1, masked=True)
  38.     with rio.open(slope_path) as f:
  39.         slope = f.read(1, masked=True)
  41.     # 计算风向
  42.     wind_dir = cal_wind_direction(u10, v10)
  43.     # 计算风向和坡向的夹角
  44.     theta_diff = np.abs(wind_dir - aspect)
  45.     theta = np.minimum(theta_diff, 360 - theta_diff)
  46.     # 计算风坡夹角
  47.     ws_angle = np.cos(np.radians(theta)) * np.sin(np.radians(slope))
  49.     with rio.open(out_path, 'w', **meta) as f:
  50.         f.write(ws_angle, 1)
复制代码
使用示例：
  1. # @Author  : ChaoQiezi
  2. # @Time    : 2025/4/21 下午10:12
  3. # @Email   : chaoqiezi.one@qq.com
  4. # @Wechat  : GIS茄子
  5. # @FileName: wind_slope_angle_cal
  7. """
  8. This script is used to 基于ERA5的u和v风向计算风向角度(0~360),再利用风向和坡度坡向计算风坡夹角.
  10. 风坡夹角: 风向与坡向夹角的余弦值和坡度正弦值的乘积.
  11. """
  12. import os.path
  14. import numpy as np
  15. from numpy import arctan2, pi
  16. from datetime import datetime
  17. from dateutil.relativedelta import relativedelta
  18. from rasterio.plot import show
  20. from Src.utils import cal_wind_direction, cal_wind_slope_angle
  22. # 准备
  23. out_dir = r'E:\Datasets\Objects\PrecipitationDownscaling\wind_slope_angle'
  24. slope_dir = r"E:\Datasets\Objects\PrecipitationDownscaling\Slope"
  25. aspect_dir = r"E:\Datasets\Objects\PrecipitationDownscaling\Aspect"
  26. u10_dir = r'E:\Datasets\Objects\PrecipitationDownscaling\u10'
  27. v10_dir = r'E:\Datasets\Objects\PrecipitationDownscaling\v10'
  28. res_folder_name = '1km'
  29. start_date = datetime(2019, 1, 1)
  30. end_date = datetime(2023, 12, 31)
  31. out_dir = os.path.join(out_dir, res_folder_name)
  32. if not os.path.exists(out_dir):
  33.     os.makedirs(out_dir)
  35. rd = relativedelta(end_date, start_date)
  36. month_range = rd.years * 12 + rd.months + 1
  37. for cur_month in range(month_range):
  38.     cur_date = start_date + relativedelta(months=cur_month)
  39.     cur_slope_path = os.path.join(slope_dir, 'slope_{}.tif'.format(res_folder_name))
  40.     cur_aspect_path = os.path.join(aspect_dir, 'aspect_{}.tif'.format(res_folder_name))
  41.     cur_u10_path = os.path.join(u10_dir, res_folder_name, 'u10_{}{:02}_{}.tif'.format(
  42.         cur_date.year, cur_date.month, res_folder_name))
  43.     cur_v10_path = os.path.join(v10_dir, res_folder_name, 'v10_{}{:02}_{}.tif'.format(
  44.         cur_date.year, cur_date.month, res_folder_name))
  45.     cur_out_filename = 'ws_angle_{}{:02}_{}.tif'.format(cur_date.year, cur_date.month, res_folder_name)
  46.     cur_out_path = os.path.join(out_dir, cur_out_filename)
  47.     cal_wind_slope_angle(cur_out_path, cur_u10_path, cur_v10_path, cur_aspect_path, cur_slope_path)
  48.     print('处理: {}'.format(cur_out_filename))
  49. print('处理完成.')
复制代码
若对博客存在疑虑或存在问题瑕疵，交流欢迎关注微信公众号：GIS茄子，或者邮箱联系(推荐 | chaoqiezi.one@qq.com).
本文由博客一文多发平台 OpenWrite 发布！

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！

使用道具 举报
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册
发帖
户烫擞
7 小时前

0

粉丝关注

13

主题发布

板块介绍填写区域，请于后台编辑
微信扫一扫