Python | WRF任意剖面风矢量的三种投影方法及绘图

用户11172986

发布于 2025-01-17 19:08:08

11000

代码可运行

文章被收录于专栏：气python风雨气python风雨

运行总次数：0

代码可运行

Python | WRF任意剖面风矢量的三种投影方法及绘图

环境设置

安装依赖

!pip install geopy -i https://pypi.mirrors.ustc.edu.cn/simple/

Looking in indexes: https://pypi.mirrors.ustc.edu.cn/simple/
Collecting geopy
  Downloading https://mirrors.ustc.edu.cn/pypi/packages/e5/15/cf2a69ade4b194aa524ac75112d5caac37414b20a3a03e6865dfe0bd1539/geopy-2.4.1-py3-none-any.whl (125 kB)
[2K     [90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━[0m [32m125.4/125.4 kB[0m [31m3.1 MB/s[0m eta [36m0:00:00[0m00:01[0m
[?25hCollecting geographiclib<3,>=1.52
  Downloading https://mirrors.ustc.edu.cn/pypi/packages/9f/5a/a26132406f1f40cf51ea349a5f11b0a46cec02a2031ff82e391c2537247a/geographiclib-2.0-py3-none-any.whl (40 kB)
[2K     [90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━[0m [32m40.3/40.3 kB[0m [31m4.8 MB/s[0m eta [36m0:00:00[0m
[?25hInstalling collected packages: geographiclib, geopy
Successfully installed geographiclib-2.0 geopy-2.4.1

摘要

本文介绍了三种将WRF模型输出的u和v风分量投影到垂直剖面上的方法：

向量点积法：通过计算风向量与剖面方向向量的点积，直接得到沿剖面的风速分量。

角度差法：通过计算风向量与剖面方向的角度差，利用余弦函数得到沿剖面的风速分量。

旋转法：通过旋转风向量，使其与剖面方向对齐，从而得到沿剖面和垂直于剖面的风速分量。

我们将通过Python代码实现这三种方法，并使用Matplotlib绘制风矢量剖面图，直观比较不同方法的结果。

项目目标

理解风场投影的基本原理。
掌握三种风场投影方法的实现。
使用Python对WRF模型输出的风场数据进行处理和分析。
通过可视化对比不同方法的结果，评估其优缺点。

投影原理

详情访问https://zhuanlan.zhihu.com/p/380626604

方法一旋转法

from zemega https://github.com/NCAR/wrf-python/issues/143

旋转法通过旋转风向量，使其与剖面方向对齐，从而得到沿剖面和垂直于剖面的风速分量。以下是实现代码

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from cartopy import crs
from cartopy.feature import NaturalEarthFeature, COLORS
from netCDF4 import Dataset
from math import acos, degrees, pi, cos, sin
from geopy.distance import geodesic
from wrf import (getvar, to_np, get_cartopy, latlon_coords, vertcross,
                 cartopy_xlim, cartopy_ylim, interpline, CoordPair)
# Rotate u and v components to align with the cross-section
def rotate_ua_va_vert_cross(ua_vertcross, va_vertcross):
    '''
    Takes u and v wind component vertcross, rotates them to align with
    transect meteorological direction, from start_point to end_point. 
    '''

    coord_pairs_1 = to_np(ua_vertcross.coords["xy_loc"])
    coord_pairs_2 = to_np(va_vertcross.coords["xy_loc"])
    if (any(coord_pairs_1 != coord_pairs_2)):
        print("u-component and v component does not match")
        return
    coord_pairs = coord_pairs_1
    main_lat = [(x.lat) for x in coord_pairs]
    main_lon = [(x.lon) for x in coord_pairs]
    # Create an emptry transect
    met_dir_transect = []

    point_a = main_lat[0], main_lon[0]
    point_b = main_lat[1], main_lon[1]
    point_c = main_lat[0], main_lon[1]
    A = geodesic(point_a, point_b).km
    B = geodesic(point_b, point_c).km
    C = geodesic(point_c, point_a).km
    degrees_A_C = 90 + degrees(acos((A * A + C * C - B * B)/(2.0 * A * C)))
    met_dir_transect.append(degrees_A_C)

    for point in range(len(main_lat)):
        if point == 0:
            continue
        point_a = main_lat[point-1], main_lon[point-1]
        point_b = main_lat[point], main_lon[point]
        point_c = main_lat[point-1], main_lon[point]
        A = geodesic(point_a, point_b).km
        B = geodesic(point_b, point_c).km
        C = geodesic(point_c, point_a).km
        degrees_A_B = 180 - \
            degrees(acos((A * A + B * B - C * C)/(2.0 * A * B)))
        met_dir_transect.append(degrees_A_B)

    met_dir_transect_2 = np.array(met_dir_transect, ndmin=1)
    # if deg == True:
    a = met_dir_transect_2/180
    # else:
    #    a = met_dir_transect_2/pi
    c = [cos(pi*X) for X in a]
    s = [sin(pi*X) for X in a]
    c_tile = np.tile(c, (len(ua_vertcross.vertical), 1))
    s_tile = np.tile(s, (len(ua_vertcross.vertical), 1))
    # if clockwise == True:
    un = ua_vertcross * c_tile - va_vertcross * s_tile
    vn = ua_vertcross * s_tile + va_vertcross * c_tile
    # else:
    #    un = ua_cross * c_tile + va_cross * s_tile
    #    vn = va_cross * c_tile - ua_cross * s_tile
    return (vn, un)


# 读取数据
path = '/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_20_00_00'
wrf_file = Dataset(path)
cross_start = CoordPair(lat=32, lon=110)
cross_end = CoordPair(lat=28, lon=115)
ter = getvar(wrf_file, "ter", timeidx=0)
ht = getvar(wrf_file, "z", timeidx=0)
w = getvar(wrf_file, "wa", timeidx=0)
u = getvar(wrf_file, "ua", timeidx=0)
v = getvar(wrf_file, "va", timeidx=0)

# 计算垂直剖面
w_cross = vertcross(w, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)
u_cross = vertcross(u, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)
v_cross = vertcross(v, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)

# 旋转风分量
h_wind_c, p_wind_c = rotate_ua_va_vert_cross(u_cross, v_cross)

# 填充垂直速度数据
w_cross_filled = np.ma.copy(to_np(w_cross))
for i in range(w_cross_filled.shape[-1]):
    column_vals = w_cross_filled[:, i]
    first_idx = int(np.transpose((column_vals > -200).nonzero())[0])
    w_cross_filled[0:first_idx, i] = w_cross_filled[first_idx, i]

# 获取地形高度
ter_line = interpline(ter, wrfin=wrf_file, start_point=cross_start, end_point=cross_end)
lats, lons = latlon_coords(w)
cart_proj = get_cartopy(w)

# 创建图形和轴对象
fig, ax_cross = plt.subplots(figsize=(8, 6),dpi=200)

# 绘制垂直速度剖面
w_levels = np.arange(-4E-1, +4E-1, 5E-2)
xs = np.arange(0, w_cross.shape[-1], 1)
ys = to_np(w_cross.coords["vertical"])
w_contours = ax_cross.contourf(xs, ys, to_np(w_cross_filled), levels=w_levels, cmap='seismic', extend="both")

# 添加颜色条
cbar = fig.colorbar(w_contours, ax=ax_cross)
cbar.ax.tick_params(labelsize=12)
cbar.set_label('Vertical Velocity (m.s-1)', rotation=-270, fontsize=12, labelpad=15)

# 填充地形区域
ht_fill = ax_cross.fill_between(xs, 0, to_np(ter_line), facecolor="black")

# 绘制风矢量
ax_cross.quiver(xs[::5], ys[::5], to_np(h_wind_c[::5, ::5]), to_np(w_cross[::5, ::5] * 100), scale=500)

# 设置x轴刻度
coord_pairs = to_np(u_cross.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape[0])
x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]
num_ticks = 5
thin = int((len(x_ticks) / num_ticks) + .5)
ax_cross.set_xticks(x_ticks[::thin])
ax_cross.set_xticklabels(x_labels[::thin], rotation=90, fontsize=8)

# 设置轴标签和范围
ax_cross.set_xlabel("Latitude, Longitude", fontsize=12)
ax_cross.set_ylabel("Height (m)", fontsize=12)
ax_cross.set_ylim([0, 10000])

# 添加网格
ax_cross.grid(True, linestyle='--', alpha=0.6)

# 显示图形
plt.tight_layout()
plt.show()

/tmp/ipykernel_139/1344016512.py:93: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  first_idx = int(np.transpose((column_vals > -200).nonzero())[0])

方法二角度差余弦法

from 爱垂钓的猫

https://www.heywhale.com/mw/project/618e8ca486227300178d33df

角度差法通过计算风向量与剖面方向的角度差，利用余弦函数得到沿剖面的风速分量。以下是实现代码

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from cartopy import crs
from cartopy.feature import NaturalEarthFeature, COLORS
from netCDF4 import Dataset
from math import acos, degrees, pi, cos, sin
from wrf import (getvar, to_np, get_cartopy, latlon_coords, vertcross,
                 cartopy_xlim, cartopy_ylim, interpline, CoordPair)


def project_wind_to_transect_angle_method(ua_vert, va_vert, start_point, end_point):
    """
    Projects u and v wind components onto a transect using the angle difference method.
    Returns the projected wind component along the transect.
    """
    # Calculate wind speed
    ws_vert = np.sqrt(ua_vert**2 + va_vert**2)
    
    # Calculate wind direction in degrees
    wdir_vert = np.arctan2(va_vert, ua_vert) * 180 / np.pi
    
    # Calculate transect direction in degrees
    dy = end_point.lat - start_point.lat
    dx = end_point.lon - start_point.lon
    line_angle = np.arctan2(dy, dx) * 180 / np.pi
    
    # Calculate angle difference between wind and transect
    angle_diff = wdir_vert - line_angle
    
    # Convert angle difference to radians and calculate projection factor
    projection_factor = np.cos(np.deg2rad(angle_diff))
    
    # Project wind onto the transect
    ws_projected = ws_vert * projection_factor
    
    return ws_projected

# 读取数据
path = '/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_20_00_00'
wrf_file = Dataset(path)

cross_start = CoordPair(lat=32, lon=110)
cross_end = CoordPair(lat=28, lon=115)
ter = getvar(wrf_file, "ter", timeidx=0)
ht = getvar(wrf_file, "z", timeidx=0)
w = getvar(wrf_file, "wa", timeidx=0)
u = getvar(wrf_file, "ua", timeidx=0)
v = getvar(wrf_file, "va", timeidx=0)

# 计算垂直剖面
w_cross = vertcross(w, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)
u_cross = vertcross(u, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)
v_cross = vertcross(v, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)

# 旋转风分量
ws_projected= project_wind_to_transect_angle_method(u_cross , v_cross, cross_start , cross_end)

# 填充垂直速度数据
w_cross_filled = np.ma.copy(to_np(w_cross))
for i in range(w_cross_filled.shape[-1]):
    column_vals = w_cross_filled[:, i]
    first_idx = int(np.transpose((column_vals > -200).nonzero())[0])
    w_cross_filled[0:first_idx, i] = w_cross_filled[first_idx, i]

# 获取地形高度
ter_line = interpline(ter, wrfin=wrf_file, start_point=cross_start, end_point=cross_end)
lats, lons = latlon_coords(w)
cart_proj = get_cartopy(w)

# 创建图形和轴对象
fig, ax_cross = plt.subplots(figsize=(8, 6),dpi=200)

# 绘制垂直速度剖面
w_levels = np.arange(-4E-1, +4E-1, 5E-2)
xs = np.arange(0, w_cross.shape[-1], 1)
ys = to_np(w_cross.coords["vertical"])
w_contours = ax_cross.contourf(xs, ys, to_np(w_cross_filled), levels=w_levels, cmap='seismic', extend="both")

# 添加颜色条
cbar = fig.colorbar(w_contours, ax=ax_cross)
cbar.ax.tick_params(labelsize=12)
cbar.set_label('Vertical Velocity (m.s-1)', rotation=-270, fontsize=12, labelpad=15)

# 填充地形区域
ht_fill = ax_cross.fill_between(xs, 0, to_np(ter_line), facecolor="black")

# 绘制风矢量
ax_cross.quiver(xs[::5], ys[::5], to_np(ws_projected[::5, ::5]), to_np(w_cross[::5, ::5] * 100), scale=500)

# 设置x轴刻度
coord_pairs = to_np(u_cross.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape[0])
x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]
num_ticks = 5
thin = int((len(x_ticks) / num_ticks) + .5)
ax_cross.set_xticks(x_ticks[::thin])
ax_cross.set_xticklabels(x_labels[::thin], rotation=90, fontsize=8)

# 设置轴标签和范围
ax_cross.set_xlabel("Latitude, Longitude", fontsize=12)
ax_cross.set_ylabel("Height (m)", fontsize=12)
ax_cross.set_ylim([0, 10000])

# 添加网格
ax_cross.grid(True, linestyle='--', alpha=0.6)

# 显示图形
plt.tight_layout()
plt.show()

/tmp/ipykernel_180/1442593371.py:64: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  first_idx = int(np.transpose((column_vals > -200).nonzero())[0])

方法三向量点积法

向量点积法通过计算风向量与剖面方向向量的点积，直接得到沿剖面的风速分量。说罢了是博主基于方法二的简化版本以下是实现代码

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from cartopy import crs
from cartopy.feature import NaturalEarthFeature, COLORS
from netCDF4 import Dataset
from math import acos, degrees, pi, cos, sin
from wrf import (getvar, to_np, get_cartopy, latlon_coords, vertcross,
                 cartopy_xlim, cartopy_ylim, interpline, CoordPair)

def project_wind_to_transect(ua_vert, va_vert, start_point, end_point):
    """
    Projects u and v wind components onto a transect using the vector dot product method.
    Returns the projected wind component along the transect and perpendicular to it.
    """
    # Calculate the direction vector of the transect
    dx = end_point.lon - start_point.lon
    dy = end_point.lat - start_point.lat
    # Calculate the length of the direction vector
    length = np.sqrt(dx**2 + dy**2)
    
    # Normalize the direction vector
    dx_norm = dx / length
    dy_norm = dy / length
    
    # Project the wind components onto the transect
    ws_along = ua_vert * dx_norm + va_vert * dy_norm  # Wind along the transect
    ws_perp = -ua_vert * dy_norm + va_vert * dx_norm  # Wind perpendicular to the transect
    
    return ws_along, ws_perp

# 读取数据
path = '/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_20_00_00'
wrf_file = Dataset(path)

cross_start = CoordPair(lat=32, lon=110)
cross_end = CoordPair(lat=28, lon=115)
ter = getvar(wrf_file, "ter", timeidx=0)
ht = getvar(wrf_file, "z", timeidx=0)
w = getvar(wrf_file, "wa", timeidx=0)
u = getvar(wrf_file, "ua", timeidx=0)
v = getvar(wrf_file, "va", timeidx=0)

# 计算垂直剖面
w_cross = vertcross(w, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)
u_cross = vertcross(u, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)
v_cross = vertcross(v, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True)

# 旋转风分量
ws_projected,_= project_wind_to_transect(u_cross , v_cross, cross_start , cross_end)

# 填充垂直速度数据
w_cross_filled = np.ma.copy(to_np(w_cross))
for i in range(w_cross_filled.shape[-1]):
    column_vals = w_cross_filled[:, i]
    first_idx = int(np.transpose((column_vals > -200).nonzero())[0])
    w_cross_filled[0:first_idx, i] = w_cross_filled[first_idx, i]

# 获取地形高度
ter_line = interpline(ter, wrfin=wrf_file, start_point=cross_start, end_point=cross_end)
lats, lons = latlon_coords(w)
cart_proj = get_cartopy(w)

# 创建图形和轴对象
fig, ax_cross = plt.subplots(figsize=(8, 6),dpi=200)

# 绘制垂直速度剖面
w_levels = np.arange(-4E-1, +4E-1, 5E-2)
xs = np.arange(0, w_cross.shape[-1], 1)
ys = to_np(w_cross.coords["vertical"])
w_contours = ax_cross.contourf(xs, ys, to_np(w_cross_filled), levels=w_levels, cmap='seismic', extend="both")

# 添加颜色条
cbar = fig.colorbar(w_contours, ax=ax_cross)
cbar.ax.tick_params(labelsize=12)
cbar.set_label('Vertical Velocity (m.s-1)', rotation=-270, fontsize=12, labelpad=15)

# 填充地形区域
ht_fill = ax_cross.fill_between(xs, 0, to_np(ter_line), facecolor="black")

# 绘制风矢量
ax_cross.quiver(xs[::5], ys[::5], to_np(ws_projected[::5, ::5]), to_np(w_cross[::5, ::5] * 100), scale=500)

# 设置x轴刻度
coord_pairs = to_np(u_cross.coords["xy_loc"])
x_ticks = np.arange(coord_pairs.shape[0])
x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]
num_ticks = 5
thin = int((len(x_ticks) / num_ticks) + .5)
ax_cross.set_xticks(x_ticks[::thin])
ax_cross.set_xticklabels(x_labels[::thin], rotation=90, fontsize=8)

# 设置轴标签和范围
ax_cross.set_xlabel("Latitude, Longitude", fontsize=12)
ax_cross.set_ylabel("Height (m)", fontsize=12)
ax_cross.set_ylim([0, 10000])

# 添加网格
ax_cross.grid(True, linestyle='--', alpha=0.6)

# 显示图形
plt.tight_layout()
plt.show()

/tmp/ipykernel_221/491975736.py:57: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  first_idx = int(np.transpose((column_vals > -200).nonzero())[0])