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

2024-08-21 08:32:45 浏览数 (2)

如何使用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环境的配置

代码语言：javascript复制
import pyart
import pydda
import matplotlib.pyplot as plt
import numpy as np
import pooch
代码语言：javascript复制
## 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

代码语言：javascript复制
## 下载数据
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)
代码语言：javascript复制
/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): 生成多部多普勒风场所使用的参数。
代码语言：javascript复制
# 设置参数开始反演
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",
)
代码语言：javascript复制
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格式，让我看看里边有啥变量

代码语言：javascript复制
ds = new_grids[0].to_xarray()
ds
代码语言：javascript复制
<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变量，反演成功

代码语言：javascript复制
ds.u
代码语言：javascript复制
<xarray.DataArray 'u' (time: 1, z: 21, y: 356, x: 151)> Size: 9MB
代码语言：javascript复制
我们得到了一个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进行多部雷达数据的风场反演， 不仅能够克服单一雷达观测的局限性， 还能在复杂天气条件下提供更准确、更详细的风场结构。 如对您有用点个赞吧

python nan 函数 数据 优化

0 人点赞