Python中坡风夹角及风效因子的测算

01 风坡夹角的界定

风坡夹角指的是风向和坡向夹角的余弦值与坡度正弦值的乘积。

02 相关说明

要计算风坡夹角，需用到ERA5 – Land再分析资料中的u10（横向风的分量）和v10（纵向风的分量）。（注：计算风速wind可通过勾股定理，对u10平方加上v10平方的结果开方得到），而坡向和坡度可借助DEM来计算（可利用ArcGIS、QGIS等软件）。

2.1 u10与v10的定义

u10全称为10 metre U wind component，v10全称是10 metre V wind component。ERA5 – Land官方文档对其定义如下：

附官方文档链接：
u10，v10

从上述描述可知：u10实际是东西向风分量朝东移动的速度，v10是南北向风分量朝北移动的速度。但依据气象学中风的定义，例如东风表示风从东边吹来，是朝正西方向移动的风，这与上述u10和v10的定义刚好相反，前者是风从哪边吹来，后者是风朝哪边吹去。所以后续计算时要尤其留意正负号的关系。

关于风向，应计算为0 – 360范围（类似东风的风向为90°），具体如下：

2.2 u10和v10的获取

获取方式：ERA5 – Land每小时产品下载

2.2.1 选择下载的变量

2.2.2 选择下载的时间

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

2.2.4 进行下载

2.3 DEM的获取

对于DEM的获取以及坡向和坡度的计算相对简单，此处不再详细阐述

03 计算过程

3.1 风向的计算

前文已提及风的相关定义，此处重点说明计算的逻辑。

由于气象学上风向的逻辑和ERA5 – Land的u10、v10的逻辑恰好相反，例如：

其他风不再一一列举。

为了便于数学上向量的计算，提前给u10和v10添上负号，其负号对应的横纵两个变量在笛卡尔xy坐标系上的合分量指向就是正确的风向。也就是说，利用arctan(\frac{-v}{-u})反算夹角就能得到风向。代码如下：

import numpy as np

def cal_wind_direction(u, v):
    # 计算来向向量 (-u, -v) 的数学角度 -- (来向向量: 风吹来的方向, 例如北风是从北边吹来的,所以来向向量的方向是正北)
    wind_dir_rad = np.arctan2(-v, -u)  # 参数顺序为 y=-v, x=-u
    """
    这里是对于np.arctan2(y, x)是传入一个向量(x, y), 计算其与向量(1, 0)(该向量与X轴正方向一致,此处坐标均为数学坐标系/二维笛卡尔坐标系)的夹角,
    夹角范围为(-180, 180), 若为一二象限即是正,若为三四象限即为负,这与数学上的逆时针为正角顺时针为负角一致.
    此外,
    此处的u和v,首先解释u表示风吹向东边的速度;v表示风吹向北边的速度;(负号表示与定义方向相反)
    由上面可以知道, u表示x轴方向上的分速度, v表示y轴方向上的分速度;
    基于u和v可以计算出风运动的方向和速度,速度这里我们暂时不管,但是对于风运动的方向,显然与向量(u, v)的方向一致;
    但是风运动的方向和向量(u, v)的方向都是表示风往哪个方向吹或者风要吹到哪里去;这与风向的定义正好相反,例如:
    北风表示从北边吹来的风,所以是吹往南边的风,假定风速是1m/s,那么u=0,v=-1m/s(由v的定义知往南吹与定义方向相反所以添上负号).
    如果直接传入arctan2(-1, 0),那么实际上得到的是向量(0, -1)的方向与向量(1, 0)方向上的夹角.
    但是我们应该是需要得到向量(0, 1)的方向与向量(1, 0)方向上的夹角.
    大家自行体会为什么需要添上负号,本质就是u和v定义正方向是风吹向的方向,而风向定义的方向是风来时或者从哪里吹来的方向--方向刚好相反
    """

    # 将弧度转换为度数 [-pi, pi] ==> [-180, 180]
    wind_dir_deg = np.degrees(wind_dir_rad)

    # 转换为地理角度(正北为0°，顺时针旋转)
    wind_dir_geo = (90 - wind_dir_deg) % 360
    """
    根据arctan2的定义知,其输出是指风向向量(当你输入是-u, -v时符合风向定义)与向量(1, 0)方向或者说x轴正方向之间的夹角
    但是风向的定义北风为0°,如此顺时针旋转到360度回到正北.两个夹角的定义不相同, 需要进行换算
    """

    return wind_dir_geo

注：函数中有两点需要注意，一是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计算得到）：

3.2 风向与坡向夹角的计算

坡向的计算可通过GIS软件如ArcGIS来实现，此处不再赘述。

风向和坡向的夹角范围是0 – 180，计算方式为：

# 计算风向和坡向的夹角
theta_diff = np.abs(wind_dir - aspect)
theta = np.minimum(theta_diff, 360 - theta_diff)

3.3 风坡夹角/风效因子的计算

# 计算风坡夹角
ws_angle = np.cos(np.radians(theta)) * np.sin(np.radians(slope))

具体函数（函数内部的cal_wind_direction为自定义函数，具体见3.1内容）：

import rasterio as rio
import numpy as np

def cal_wind_slope_angle(out_path, u10_path, v10_path, aspect_path, slope_path):
    """
    基于u10和v10计算风向, 基于风向和坡向、坡度计算风坡夹角
    风坡夹角: 风向与坡向夹角的余弦值和坡度正弦值的乘积
    :param out_path: 风坡夹角的输出路径
    :param u10_path: u10的输入路径
    :param v10_path: v10的输入路径
    :param aspect_path: 坡向的输入路径
    :param slope_path: 坡度的输入路径
    :return: None
    """

    # 读取经纬向风速,坡向坡度
    with rio.open(u10_path) as f:
        u10 = f.read(1, masked=True)  # 读取第一个波段
        meta = f.meta
        """
        .meta返回当前tiff文件的元数据, 包括格格式(GTiff)、数据类型(dtype)、无效值(nodata)、行列数和波段数(width,height,count),
        坐标参考系(crs)、仿射参数(transform)
        示例:
        {'driver': 'GTiff',
         'dtype': 'float32',
         'nodata': nan,
         'width': 129,
         'height': 133,
         'count': 1,
         '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]]'),
         'transform': Affine(0.1, 0.0, 97.30000000000001,
                0.0, -0.1, 34.4)}
        """
    with rio.open(v10_path) as f:
        v10 = f.read(1, masked=True)
    with rio.open(aspect_path) as f:
        aspect = f.read(1, masked=True)
    with rio.open(slope_path) as f:
        slope = f.read(1, masked=True)

    # 计算风向
    wind_dir = cal_wind_direction(u10, v10)
    # 计算风向和坡向的夹角
    theta_diff = np.abs(wind_dir - aspect)
    theta = np.minimum(theta_diff, 360 - theta_diff)
    # 计算风坡夹角
    ws_angle = np.cos(np.radians(theta)) * np.sin(np.radians(slope))

    with rio.open(out_path, 'w', **meta) as f:
        f.write(ws_angle, 1)

使用示例：

# @Author  : ChaoQiezi
# @Time    : 2025/4/21 下午10:12
# @Email   : chaoqiezi.one@qq.com
# @Wechat  : GIS茄子
# @FileName: wind_slope_angle_cal

"""
This script is used to 基于ERA5的u和v风向计算风向角度(0~360),再利用风向和坡度坡向计算风坡夹角.

风坡夹角: 风向与坡向夹角的余弦值和坡度正弦值的乘积.
"""
import os.path

import numpy as np
from numpy import arctan2, pi
from datetime import datetime
from dateutil.relativedelta import relativedelta
from rasterio.plot import show

from Src.utils import cal_wind_direction, cal_wind_slope_angle

# 准备
out_dir = r'E:\Datasets\Objects\PrecipitationDownscaling\wind_slope_angle'
slope_dir = r"E:\Datasets\Objects\PrecipitationDownscaling\Slope"
aspect_dir = r"E:\Datasets\Objects\PrecipitationDownscaling\Aspect"
u10_dir = r'E:\Datasets\Objects\PrecipitationDownscaling\u10'
v10_dir = r'E:\Datasets\Objects\PrecipitationDownscaling\v10'
res_folder_name = '1km'
start_date = datetime(2019, 1, 1)
end_date = datetime(2023, 12, 31)
out_dir = os.path.join(out_dir, res_folder_name)
if not os.path.exists(out_dir):
    os.makedirs(out_dir)

rd = relativedelta(end_date, start_date)
month_range = rd.years * 12 + rd.months + 1
for cur_month in range(month_range):
    cur_date = start_date + relativedelta(months=cur_month)
    cur_slope_path = os.path.join(slope_dir, 'slope_{}.tif'.format(res_folder_name))
    cur_aspect_path = os.path.join(aspect_dir, 'aspect_{}.tif'.format(res_folder_name))
    cur_u10_path = os.path.join(u10_dir, res_folder_name, 'u10_{}{:02}_{}.tif'.format(
        cur_date.year, cur_date.month, res_folder_name))
    cur_v10_path = os.path.join(v10_dir, res_folder_name, 'v10_{}{:02}_{}.tif'.format(
        cur_date.year, cur_date.month, res_folder_name))
    cur_out_filename = 'ws_angle_{}{:02}_{}.tif'.format(cur_date.year, cur_date.month, res_folder_name)
    cur_out_path = os.path.join(out_dir, cur_out_filename)
    cal_wind_slope_angle(cur_out_path, cur_u10_path, cur_v10_path, cur_aspect_path, cur_slope_path)
    print('处理: {}'.format(cur_out_filename))
print('处理完成.')

需注意，原文中最后的广告相关内容已被过滤。