Optimizing your wind retrieval#
In the Retrieving your first wind field section, we showed how to perform an example wind retrieval with PyDDA. However, there were some issues to be resolved in the wind retrieval, including artificial updrafts at the boundaries of the Dual Doppler lobes caused by discontinuities in the horizontal winds from changing data sources from the radar network to the model constraint. In this section, we will show how to adjust the parameters of your wind retrieval in order to minimize such artifacts. First, we will show the wind retrieval that we did in Retrieving your first wind field.
grids_out, _ = pydda.retrieval.get_dd_wind_field([grid_kict, grid_ktlx],
Cm=256.0, Co=1e-2, Cx=1, Cy=1,
Cz=1, Cmod=1e-5, model_fields=["hrrr"],
refl_field='DBZ', wind_tol=0.5,
max_iterations=100, filter_window=15,
filter_order=3, engine='scipy')
pydda.vis.plot_horiz_xsection_quiver(grids_out, level=15, cmap='ChaseSpectral', vmin=-10, vmax=80,
quiverkey_len=20.0, background_field='DBZ', bg_grid_no=1,
w_vel_contours=[1, 2, 5, 10], quiver_spacing_x_km=10.0,
quiver_spacing_y_km=10.0, quiverkey_loc='bottom_right')
(Source code, png, hires.png, pdf)
We can see several potential issues with the wind retrieval. First, there are artifacts at the
Dual Doppler lobe edges where updrafts are being produced by the optimization code simply
because of a discontinuity in the horizontal winds at the edges of the lobes. In addition,
there are other discontinuities in the horizontal winds that should be addressed. One thing
we can do to mitigate these discontinuities is to increase the weight of the horizontal
smoothness constraints. Therefore, let’s prescribe Cx = 100. and Cy = 100
to the above retrieval.
grids_out, _ = pydda.retrieval.get_dd_wind_field([grid_kict, grid_ktlx],
Cm=256.0, Co=1e-2, Cx=100, Cy=100,
Cz=1, Cmod=1e-5, model_fields=["hrrr"],
refl_field='DBZ', wind_tol=0.5,
max_iterations=100, filter_window=15,
filter_order=3, engine='scipy')
pydda.vis.plot_horiz_xsection_quiver(grids_out, level=15, cmap='ChaseSpectral', vmin=-10, vmax=80,
quiverkey_len=20.0, background_field='DBZ', bg_grid_no=1,
w_vel_contours=[1, 2, 5, 10], quiver_spacing_x_km=10.0,
quiver_spacing_y_km=10.0, quiverkey_loc='bottom_right')
(Source code, png, hires.png, pdf)
As we can see, the artifact at the edge of the Dual Doppler lobe has reduced in size. However, we also have lost some detail on the updraft structure at this level because the wind field has been smoothed out. This therefore coarsens the effective resolution of the retrieval. Let’s see what happens when we increase the level of smoothing.
grids_out, _ = pydda.retrieval.get_dd_wind_field([grid_kict, grid_ktlx],
Cm=256.0, Co=1e-2, Cx=250., Cy=250.,
Cz=250.0, Cmod=1e-5, model_fields=["hrrr"],
refl_field='DBZ', wind_tol=0.5,
max_iterations=100, filter_window=15,
filter_order=3, engine='scipy')
pydda.vis.plot_horiz_xsection_quiver(grids_out, level=15, cmap='ChaseSpectral', vmin=-10, vmax=80,
quiverkey_len=20.0, background_field='DBZ', bg_grid_no=1,
w_vel_contours=[1, 2, 5, 10], quiver_spacing_x_km=10.0,
quiver_spacing_y_km=10.0, quiverkey_loc='bottom_right')
(Source code, png, hires.png, pdf)
In the above retrieval, the updrafts appear to be smoothed out. To help the optimization loop resolve the updrafts, we recommend, from here, decreasing the tolerance required for the optimization loop to converge. In addition, decreasing the smoothness will allow more details of the resolved wind field to appear. In the below example, we observe this, though part of the artifact near the edge of the Dual Doppler lobe re-appears.
grids_out, _ = pydda.retrieval.get_dd_wind_field([grid_kict, grid_ktlx],
Cm=256.0, Co=1e-2, Cx=150., Cy=150.,
Cz=150.0, Cmod=1e-5, model_fields=["hrrr"],
refl_field='DBZ', wind_tol=0.1,
max_iterations=400, filter_window=15,
filter_order=3, engine='scipy')
pydda.vis.plot_horiz_xsection_quiver(grids_out, level=15, cmap='ChaseSpectral', vmin=-10, vmax=80,
quiverkey_len=20.0, background_field='DBZ', bg_grid_no=1,
w_vel_contours=[1, 2, 5, 10], quiver_spacing_x_km=10.0,
quiver_spacing_y_km=10.0, quiverkey_loc='bottom_right')
(Source code, png, hires.png, pdf)
Generally, these parameters need to be tuned for your particular radar configuration in order to obtain the most optimal wind retrieval for your situation. If you are placing more importance on horizontal winds compared to updraft velocities, then you may be willing to tolerate more errors in the vertical velocity field so that finer details of the horizontal wind field can be generated. The above parameters are examples that apply to a 1 km resolution grid from two NEXRADs and vary for given radar configurations and storm coverages.
Choosing an upper boundary condition#
PyDDA constrains the vertical velocity at the boundaries of the analysis domain by
zeroing the gradient of the cost function with respect to \(w\) there, so that
\(w\) is held at its first guess. At the surface this impermeability condition is
always applied. At the top of the domain it is controlled by the upper_bc
keyword of pydda.retrieval.get_dd_wind_field(), which takes one of three values:
|
Condition applied |
|---|---|
|
No condition at the top of the domain. |
|
\(w = 0\) at the top vertical level of the grid. |
|
\(w = 0\) above the echo top, i.e. wherever no radar reports an observation and the point is higher than above km in the grid’s vertical coordinate. |
The classic choice is upper_bc=1, which imposes \(w = 0\) at the top of the
analysis domain. That is only physically defensible when the domain top is genuinely
above the storm. If the grid is truncated through the middle of deep convection, the
condition forces the mass continuity constraint to close the divergence profile at an
arbitrary height and pushes a spurious compensating signal down into the levels you
actually care about.
Setting upper_bc=2 instead applies the impermeability condition at the echo
top: every grid point at which none of the radars report a valid radial velocity is
treated as being outside the storm, and \(w\) is held at its first guess there.
Because the first guess for \(w\) is normally zero, this is equivalent to requiring
that no mass crosses the top of the observed echo. The method is described in Thompson
et al. (2026), https://doi.org/10.5194/egusphere-2026-4631.
The above keyword sets the lowest altitude, in km, at which the echo top
condition may be applied. It exists so that clear air at low levels – gaps between
cells, the cone of silence, the far edge of the Dual Doppler lobes – does not pin
\(w\) to zero close to the surface, which would suppress the very updrafts you are
trying to retrieve. The default of 2 km is a reasonable starting point; raise it if
your radars have poor low-level coverage.
grids_out, _ = pydda.retrieval.get_dd_wind_field([grid_kict, grid_ktlx],
Cm=256.0, Co=1e-2, Cx=1, Cy=1,
Cz=1, Cmod=1e-5, model_fields=["hrrr"],
refl_field='DBZ', wind_tol=0.5,
max_iterations=50, filter_window=15,
filter_order=3, engine='scipy',
upper_bc=2, above=2.0)
Both upper_bc and above are supported by all of PyDDA’s engines
("scipy", "jax", "tensorflow" and "auglag"). The set of
points at which the condition is applied is calculated once at the start of the
retrieval by pydda.cost_functions.calculate_echo_top_mask() and is returned on the
upper_bc_mask attribute of the parameters object, so it can be inspected
afterwards:
grids_out, parameters = pydda.retrieval.get_dd_wind_field(
[grid_kict, grid_ktlx], upper_bc=2, above=2.0, ...)
print("Fraction of the domain held impermeable: %.2f"
% parameters.upper_bc_mask.mean())
Note
For backwards compatibility upper_bc=True and upper_bc=False are
still accepted and mean the same as 1 and 0.