雷达系列：如何使用python进行多部雷达数据反演风

用户11172986

发布于 2024-08-21 08:32:45

76000

代码可运行

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

运行总次数：0

代码可运行

如何使用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, )