Example of a wind retrieval in a tornado over Sydney#

This shows an example of how to retrieve winds from 4 radars over Sydney.

We use smoothing to decrease the magnitude of the updraft in the region of the mesocyclone. The reduction of noise also helps the solution converge much faster since the cost function is smoother and therefore less susecptible to find a local minimum that is in noise.

The observational constraint is reduced to 0.01 from the usual 1 because we are factoring in many more data points as we are using 4 radars instead of the two in the Darwin example.

This example uses pooch to download the data files.

import pydda
import matplotlib.pyplot as plt
import numpy as np


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 = pydda.io.read_grid(grid1_path)
grid2 = pydda.io.read_grid(grid2_path)
grid3 = pydda.io.read_grid(grid3_path)
grid4 = pydda.io.read_grid(grid4_path)

# Set initialization and do retrieval
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=10,
    Cy=10,
    Cz=10,
    vel_name="VRADH_corr",
    refl_field="DBZH",
    mask_outside_opt=True,
    wind_tol=0.5,
    max_iterations=200,
    engine="scipy",
)
# Make a neat plot
fig = plt.figure(figsize=(10, 7))
ax = pydda.vis.plot_horiz_xsection_quiver_map(
    new_grids,
    background_field="DBZH",
    level=3,
    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.0, 151.35])
ax.set_ylim([-34.15, -33.9])
plt.show()
## 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 2.1.0
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!
False
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
Total number of model points: 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.1228988818772181
Nfeval | Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Jpoint  | Max w  
     10| 117.2266|  23.2756|   0.0000|   0.0000|   0.0000|   0.0000|   0.0000|  18.4482
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[1], line 17
     15 # Set initialization and do retrieval
     16 grid1 = pydda.initialization.make_constant_wind_field(grid1, vel_field="VRADH_corr")
---> 17 new_grids, _ = pydda.retrieval.get_dd_wind_field(
     18     [grid1, grid2, grid3, grid4],
     19     Co=1e-2,
     20     Cm=256.0,
     21     Cx=10,
     22     Cy=10,
     23     Cz=10,
     24     vel_name="VRADH_corr",
     25     refl_field="DBZH",
     26     mask_outside_opt=True,
     27     wind_tol=0.5,
     28     max_iterations=200,
     29     engine="scipy",
     30 )
     31 # Make a neat plot
     32 fig = plt.figure(figsize=(10, 7))

File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1474, in get_dd_wind_field(Grids, u_init, v_init, w_init, engine, **kwargs)
   1467     w_init = new_grids[0]["w"].values.squeeze()
   1469 if (
   1470     engine.lower() == "scipy"
   1471     or engine.lower() == "jax"
   1472     or engine.lower() == "auglag"
   1473 ):
-> 1474     return _get_dd_wind_field_scipy(
   1475         new_grids, u_init, v_init, w_init, engine, **kwargs
   1476     )
   1477 elif engine.lower() == "tensorflow":
   1478     return _get_dd_wind_field_tensorflow(
   1479         new_grids, u_init, v_init, w_init, **kwargs
   1480     )

File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:602, in _get_dd_wind_field_scipy(Grids, u_init, v_init, w_init, engine, points, vel_name, refl_field, u_back, v_back, z_back, frz, Co, Cm, Cx, Cy, Cz, Cb, Cv, Cmod, Cpoint, cvtol, gtol, Jveltol, Ut, Vt, low_pass_filter, mask_outside_opt, weights_obs, weights_model, weights_bg, max_iterations, mask_w_outside_opt, filter_window, filter_order, min_bca, max_bca, upper_bc, model_fields, output_cost_functions, roi, wind_tol, tolerance, const_boundary_cond, max_wind_mag)
    600 parameters.print_out = False
    601 if engine.lower() == "scipy":
--> 602     winds = fmin_l_bfgs_b(
    603         J_function,
    604         winds,
    605         args=(parameters,),
    606         maxiter=max_iterations,
    607         pgtol=tolerance,
    608         bounds=bounds,
    609         fprime=grad_J,
    610         disp=0,
    611         iprint=-1,
    612         callback=_vert_velocity_callback,
    613     )
    614 else:
    616     def loss_and_gradient(x):

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_lbfgsb_py.py:237, in fmin_l_bfgs_b(func, x0, fprime, args, approx_grad, bounds, m, factr, pgtol, epsilon, iprint, maxfun, maxiter, disp, callback, maxls)
    225 callback = _wrap_callback(callback)
    226 opts = {'disp': disp,
    227         'iprint': iprint,
    228         'maxcor': m,
   (...)
    234         'callback': callback,
    235         'maxls': maxls}
--> 237 res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
    238                        **opts)
    239 d = {'grad': res['jac'],
    240      'task': res['message'],
    241      'funcalls': res['nfev'],
    242      'nit': res['nit'],
    243      'warnflag': res['status']}
    244 f = res['fun']

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_lbfgsb_py.py:407, in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)
    401 task_str = task.tobytes()
    402 if task_str.startswith(b'FG'):
    403     # The minimization routine wants f and g at the current x.
    404     # Note that interruptions due to maxfun are postponed
    405     # until the completion of the current minimization iteration.
    406     # Overwrite f and g:
--> 407     f, g = func_and_grad(x)
    408 elif task_str.startswith(b'NEW_X'):
    409     # new iteration
    410     n_iterations += 1

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:344, in ScalarFunction.fun_and_grad(self, x)
    342     self._update_x(x)
    343 self._update_fun()
--> 344 self._update_grad()
    345 return self.f, self.g

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:306, in ScalarFunction._update_grad(self)
    304 if self._orig_grad in FD_METHODS:
    305     self._update_fun()
--> 306 self.g = self._wrapped_grad(self.x, f0=self.f)
    307 self.g_updated = True

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:41, in _wrapper_grad.<locals>.wrapped(x, **kwds)
     38 def wrapped(x, **kwds):
     39     # kwds present to give function same signature as numdiff variant
     40     ncalls[0] += 1
---> 41     return np.atleast_1d(grad(np.copy(x), *args))

File ~/work/PyDDA/PyDDA/pydda/cost_functions/cost_functions.py:541, in grad_J(winds, parameters)
    528     grad += _cost_functions_numpy.calculate_mass_continuity_gradient(
    529         winds[0],
    530         winds[1],
   (...)
    537         upper_bc=parameters.upper_bc,
    538     )
    540 if parameters.Cx > 0 or parameters.Cy > 0 or parameters.Cz > 0:
--> 541     grad += _cost_functions_numpy.calculate_smoothness_gradient(
    542         winds[0],
    543         winds[1],
    544         winds[2],
    545         parameters.dx,
    546         parameters.dy,
    547         parameters.dz,
    548         Cx=parameters.Cx,
    549         Cy=parameters.Cy,
    550         Cz=parameters.Cz,
    551         upper_bc=parameters.upper_bc,
    552     )
    554 if parameters.Cb > 0:
    555     grad += _cost_functions_numpy.calculate_background_gradient(
    556         winds[0],
    557         winds[1],
   (...)
    562         parameters.Cb,
    563     )

File ~/work/PyDDA/PyDDA/pydda/cost_functions/_cost_functions_numpy.py:246, in calculate_smoothness_gradient(u, v, w, dx, dy, dz, Cx, Cy, Cz, upper_bc)
    244 grad_w = np.zeros(w.shape)
    245 scipy.ndimage.laplace(u, du, mode="wrap")
--> 246 scipy.ndimage.laplace(v, dv, mode="wrap")
    247 scipy.ndimage.laplace(w, dw, mode="wrap")
    248 du = du / dx

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/ndimage/_filters.py:593, in laplace(input, output, mode, cval)
    591 def derivative2(input, axis, output, mode, cval):
    592     return correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0)
--> 593 return generic_laplace(input, derivative2, output, mode, cval)

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/ndimage/_filters.py:550, in generic_laplace(input, derivative2, output, mode, cval, extra_arguments, extra_keywords)
    548 if len(axes) > 0:
    549     modes = _ni_support._normalize_sequence(mode, len(axes))
--> 550     derivative2(input, axes[0], output, modes[0], cval,
    551                 *extra_arguments, **extra_keywords)
    552     for ii in range(1, len(axes)):
    553         tmp = derivative2(input, axes[ii], output.dtype, modes[ii], cval,
    554                           *extra_arguments, **extra_keywords)

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/ndimage/_filters.py:592, in laplace.<locals>.derivative2(input, axis, output, mode, cval)
    591 def derivative2(input, axis, output, mode, cval):
--> 592     return correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0)

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/ndimage/_filters.py:140, in correlate1d(input, weights, axis, output, mode, cval, origin)
    136     raise ValueError('Invalid origin; origin must satisfy '
    137                      '-(len(weights) // 2) <= origin <= '
    138                      '(len(weights)-1) // 2')
    139 mode = _ni_support._extend_mode_to_code(mode)
--> 140 _nd_image.correlate1d(input, weights, axis, output, mode, cval,
    141                       origin)
    142 return output

KeyboardInterrupt: