如何使用python进行多部雷达数据反演风
前言
之前在公众号气python风雨发的素材募集中，有读者询问如何使用多部雷达反演风场
这个好办，前人早已写了库
雷达反演风一直是难题，今天我们介绍一个利用雷达数据反演风的库：pydda
需要注意的是PyDDA 现在支持使用 Jax 和 TensorFlow 来求解三维风场。PyDDA 需要启用 TensorFlow 2.6 和 TensorFlow 的 tensorflow-probability 包。此外，这两个软件包都可以利用支持 CUDA 的 GPU 来加快处理速度。这两个依赖项是可选的，因为用户仍然可以在 SciPy 生态系统中使用 PyDDA。Jax 优化器使用与 SciPy 的 （L-BFGS-B） 相同的优化器。
以下是它的官方示例之一
这显示了如何从悉尼上空的 4 个雷达中检索风的示例。
我们使用平滑来降低中气旋区域的上升气流的幅度。噪声的降低还有助于解更快地收敛，因为成本函数更平滑，因此更难找到噪声中的局部最小值。
观测约束从通常的 1 减少到 0.01，因为我们使用了 4 个雷达，因此我们考虑了更多的数据点。
此示例使用 pooch 下载数据文件。
使用镜像：气象分析3.9
资源：4核16g
ps：由于基于cpu计算，运行可能需较长时间
如需在线运行可点击原文链接
导入库
此处感谢龙清大佬对pydda环境的配置
import pyart
import pydda
import matplotlib.pyplot as plt
import numpy as np
import pooch## You are using the Python ARM Radar Toolkit (Py-ART), an open source
## library for working with weather radar data. Py-ART is partly
## supported by the U.S. Department of Energy as part of the Atmospheric
## Radiation Measurement (ARM) Climate Research Facility, an Office of
## Science user facility.
##
## If you use this software to prepare a publication, please cite:
##
## JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119
Welcome to PyDDA 1.5.1
If you are using PyDDA in your publications, please cite:
Jackson et al. (2020) Journal of Open Research Science
Detecting Jax...
Jax/JaxOpt are not installed on your system, unable to use Jax engine.
Detecting TensorFlow...
TensorFlow detected. Checking for tensorflow-probability...
TensorFlow-probability detected. TensorFlow engine enabled!数据下载与读取
pydda是与pyart库搭配计算的，如有不同格式的雷达数据需要转为pyart可以读取的格式
小编才知道悉尼英文是sydney
## 下载数据
grid1_path = pydda.tests.get_sample_file("grid1_sydney.nc")
grid2_path = pydda.tests.get_sample_file("grid2_sydney.nc")
grid3_path = pydda.tests.get_sample_file("grid3_sydney.nc")
grid4_path = pydda.tests.get_sample_file("grid4_sydney.nc")
## 读取四个雷达
grid1 = pyart.io.read_grid(grid1_path)
grid2 = pyart.io.read_grid(grid2_path)
grid3 = pyart.io.read_grid(grid3_path)
grid4 = pyart.io.read_grid(grid4_path)/opt/conda/lib/python3.9/site-packages/pyart/io/cfradial.py:412: UserWarning: WARNING: valid_min not used since it
cannot be safely cast to variable data type
data = self.ncvar[:]
/opt/conda/lib/python3.9/site-packages/pyart/io/cfradial.py:412: UserWarning: WARNING: valid_max not used since it
cannot be safely cast to variable data type
data = self.ncvar[:]反演部分
可通过调整不同的参数来优化反演结果，另外计算时间较长，小编还不知如何调用gpu计算，还请知道的大佬指点指点pydda.retrieval.get_dd_wind_field() 函数接受一系列 Py-ART 格网对象作为输入，从中推导出风场。该函数要求所有输入的 Py-ART 格网必须具有相同的格网规格，即形状、X 坐标、Y 坐标和 Z 坐标都相同。
参数说明：
• Grids (list of Py-ART/DDA Grids): 一系列 Py-ART 或 PyDDA 格网，对应每一部雷达。所有格网必须有相同的形状和坐标。
• u_init (3D ndarray): 经向风场的初始猜测，作为与 Grids 中字段相同形状的 3D 数组输入。如果为 None，则 PyDDA 将使用第一个 Grid 中的 u 字段作为初始化。
• v_init (3D ndarray): 纬向风场的初始猜测，作为与 Grids 中字段相同形状的 3D 数组输入。如果为 None，则 PyDDA 将使用第一个 Grid 中的 v 字段作为初始化。
• w_init (3D ndarray): 垂直风场的初始猜测，作为与 Grids 中字段相同形状的 3D 数组输入。如果为 None，则 PyDDA 将使用第一个 Grid 中的 w 字段作为初始化。
• engine (str): 设定此标志将使用基于 SciPy、TensorFlow 或 Jax 的求解器。使用 TensorFlow 或 Jax 能够让 PyDDA 利用基于 GPU 的系统。此外，这两种实现使用自动微分计算代价函数的梯度以优化梯度计算。
• points (None or list of dicts): 由
pydda.constraints.get_iem_obs()返回的点观测。设置为 None 禁用点观测。• vel_name (string): 径向速度场名称。设置为 None 将使 PyDDA 自动检测速度场名称。
• refl_field (string): 反射率场名称。设置为 None 将使 PyDDA 自动检测反射率场名称。
• u_back (1D array): 随高度变化的背景经向风场，从探空数据获取。
• v_back (1D array): 随高度变化的背景纬向风场，从探空数据获取。
• z_back (1D array): 对应背景风场级别的高度，单位米。
• frz (float): 冻结层高度，用于计算下落速度，单位米。
• Co (float): 与观测径向速度相关的代价函数权重。
• Cm (float): 与质量连续性方程相关的代价函数权重。
• Cx (float): 与 X 方向平滑度相关的代价函数权重。
• Cy (float): 与 Y 方向平滑度相关的代价函数权重。
• Cz (float): 与 Z 方向平滑度相关的代价函数权重。
• Cv (float): 与垂直涡度方程相关的代价函数权重。
• Cmod (float): 与自定义约束相关的代价函数权重。
• Cpoint (float): 与点观测相关的代价函数权重。
• weights_obs (list of floating point arrays or None): 来自 Grids 中每部雷达的网格中每一点的权重列表。设为 None 让 PyDDA 自动确定。
• weights_model (list of floating point arrays or None): 来自 model_fields 中每项自定义场的网格中每一点的权重列表。设为 None 让 PyDDA 自动确定。
• weights_bg (list of floating point arrays or None): 来自探空数据的网格中每一点的权重列表。设为 None 让 PyDDA 自动确定。
• Ut (float): 规定的风暴运动经向分量。如果 Cv 不为零则需要此参数。
• Vt (float): 规定的风暴运动纬向分量。如果 Cv 不为零则需要此参数。
• filter_winds (bool): 若为 True，PyDDA 将对检索到的风场运行低通滤波器。设为 False 禁用低通滤波器。
• mask_outside_opt (bool): 若设为 True，风值在多部多普勒波瓣外将被屏蔽，即如果少于 2 部雷达覆盖某一点。
• max_iterations (int): 优化循环的最大迭代次数。
• mask_w_outside_opt (bool): 若设为 True，垂直风在多部多普勒波瓣外将被屏蔽，即如果少于 2 部雷达覆盖某一点。
• filter_window (int): 低通滤波器的窗口大小。较大的窗口将增加多项式拟合中考虑的点数，从而增加平滑度。
• filter_order (int): 低通滤波器使用的多项式阶数。高阶多项式允许保留更小尺度特征但可能也无法去除足够噪声。
• min_bca (float): 两部雷达之间的最小波束交叉角，单位度。30.0 是文献中常用的值。
• max_bca (float): 两部雷达之间的最大波束交叉角，单位度。150.0 是文献中常用的值。
• upper_bc (bool): 设为 True 强制顶层大气 w = 0，即所谓的不可渗透条件。
• model_fields (list of strings): 第一个 Grid 中包含自定义数据的字段列表，这些数据已插值到 Grid 的格网规格。PyDDA 将寻找形如 U_(model field name), V_(model field name), 和 W_(model field name) 的字段。
• output_cost_functions (bool): 设为 True 在每 10 次迭代后输出每个代价函数的值。
• roi (float): 点观测的影响半径。点观测在该半径之外将不具有任何权重。
• parallel_iterations (int): 优化循环中并行运行的迭代次数。仅对 TensorFlow 基础引擎有效。
• wind_tol (float): 最大风速变化小于此值时停止迭代。
• tolerance (float): L2 范数的梯度容忍度，在此之前停止。
• max_wind_magnitude (float): 约束优化使得 |u|, |w|, 和 |w| < x m/s。
返回值：
• new_grid_list (list): 包含推导出风场的 Py-ART 格网列表。这些字段可由可视化模块显示。
• parameters (struct): 生成多部多普勒风场所使用的参数。
# 设置参数开始反演
grid1 = pydda.initialization.make_constant_wind_field(grid1, vel_field="VRADH_corr")
new_grids, _ = pydda.retrieval.get_dd_wind_field(
[grid1, grid2, grid3, grid4],
Co=1e-2,
Cm=256.0,
Cx=1e-4,
Cy=1e-4,
Cz=1e-4,
vel_name="VRADH_corr",
refl_field="DBZH",
mask_outside_opt=True,
wind_tol=0.1,
max_iterations=200,
engine="scipy",
)
Calculating weights for radars 0 and 1
Calculating weights for radars 0 and 2
Calculating weights for radars 0 and 3
Calculating weights for radars 1 and 0
Calculating weights for radars 1 and 2
Calculating weights for radars 1 and 3
Calculating weights for radars 2 and 0
Calculating weights for radars 2 and 1
Calculating weights for radars 2 and 3
Calculating weights for radars 3 and 0
Calculating weights for radars 3 and 1
Calculating weights for radars 3 and 2
Calculating weights for models...
Starting solver
rmsVR = 7.242250167336596
Total points: 281276
The max of w_init is 0.0
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
0|2864.7900| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000
The gradient of the cost functions is 0.10783124122961905
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
10| 57.7027| 17.0172| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 51.2499
/opt/conda/lib/python3.9/site-packages/pydda/retrieval/wind_retrieve.py:600: RuntimeWarning: All-NaN slice encountered
delta = np.nanmax(diff)
Max change in w: nan
The gradient of the cost functions is 0.13319518639944616
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
20| 54.9807| 18.7575| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 75.2207
The gradient of the cost functions is 0.06265389287213236
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
30| 62.5696| 40.3943| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 100.0000
The gradient of the cost functions is 0.05768538154971628
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
40| 54.2658| 19.1128| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 78.5232
The gradient of the cost functions is 0.05505137139375907
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
50| 53.5801| 19.9986| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 82.7234
The gradient of the cost functions is 0.05439831132705333
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
60| 53.5883| 19.6409| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 82.0780
The gradient of the cost functions is 0.05374405443091169
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
Done! Time = 855.1得到的结果依然是pyart的grid格式，让我看看里边有啥变量
ds = new_grids[0].to_xarray()
ds<xarray.Dataset> Size: 60MB
Dimensions: (time: 1, z: 21, y: 356, x: 151)
Coordinates:
• time (time) object 8B 2015-12-15 23:36:32
• z (z) float64 168B 0.0 1e+03 2e+03 3e+03 ... 1.8e+04 1.9e+04 2e+04
lat (y, x) float64 430kB -36.11 -36.11 -36.11 ... -32.93 -32.93
lon (y, x) float64 430kB 150.6 150.6 150.6 ... 152.2 152.2 152.2• y (y) float64 3kB -5e+04 -4.901e+04 -4.803e+04 ... 2.99e+05 3e+05
• x (x) float64 1kB 1e+05 1.01e+05 1.02e+05 ... 2.49e+05 2.5e+05
Data variables:
DBZH (time, z, y, x) float32 5MB nan nan nan nan ... nan nan nan nan
VRADH_corr (time, z, y, x) float32 5MB nan nan nan nan ... nan nan nan nan
ROI (time, z, y, x) float32 5MB 1.171e+03 1.18e+03 ... 4.09e+03
AZ (time, z, y, x) float64 9MB 116.6 116.3 116.1 ... 39.69 39.81
EL (time, z, y, x) float64 9MB -1.086 -1.083 -1.081 ... 1.416 1.409
u (time, z, y, x) float64 9MB nan nan nan nan ... nan nan nan nan
v (time, z, y, x) float64 9MB nan nan nan nan ... nan nan nan nan
w (time, z, y, x) float64 9MB nan nan nan nan ... nan nan nan nan
可以看到文件中多了u和v变量，反演成功
ds.u<xarray.DataArray 'u' (time: 1, z: 21, y: 356, x: 151)> Size: 9MB我们得到了一个3d的风场
绘制反演结果
使用matplotlib绘制反演后的水平切面风矢图，背景展示反射率因子，同时标注风速等值线。
平面图
# Make a neat plot
fig = plt.figure(figsize=(20, 12))
ax = pydda.vis.plot_horiz_xsection_quiver_map(
new_grids,
background_field="DBZH",
level=5,
show_lobes=False,
bg_grid_no=3,
vmin=0,
vmax=60,
quiverkey_len=20.0,
w_vel_contours=[1.0, 3.0, 5.0, 10.0, 20.0],
quiver_spacing_x_km=2.0,
quiver_spacing_y_km=2.0,
quiverkey_loc="top",
colorbar_contour_flag=True,
cmap="ChaseSpectral",
)
ax.set_xticks(np.arange(150.5, 153, 0.1))
ax.set_yticks(np.arange(-36, -32.0, 0.1))
ax.set_xlim([151, 152.3])
ax.set_ylim([-34.15, -32.8])
plt.show()剖面图
plt.figure(figsize=(9, 9))
pydda.vis.plot_xz_xsection_barbs(
new_grids,
background_field="DBZH",
level=40,
w_vel_contours=[5, 10, 15],
barb_spacing_x_km=10.0,
barb_spacing_z_km=2.0,
)
plt.show()
# Plot a vertical Y-Z cross section
plt.figure(figsize=(9, 9))
pydda.vis.plot_yz_xsection_barbs(
new_grids,
background_field="DBZH",
level=40,
barb_spacing_y_km=10.0,
barb_spacing_z_km=2.0,
)<Axes: title={'center': 'PyDDA retreived winds @140.0 km east of origin.'}, xlabel='Y [km]', ylabel='Z [km]'>小结
通过pydda进行多部雷达数据的风场反演，
不仅能够克服单一雷达观测的局限性，
还能在复杂天气条件下提供更准确、更详细的风场结构。
如对您有用点个赞吧
