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

<frozen importlib._bootstrap>:283: DeprecationWarning: the load_module() method is deprecated and slated for removal in Python 3.12; use exec_module() instead

/usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/pyart/testing/example_data.py:1: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  import pkg_resources
/usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/pkg_resources/__init__.py:2832: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('google')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  declare_namespace(pkg)
/usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/pkg_resources/__init__.py:2832: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('sphinxcontrib')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  declare_namespace(pkg)
/usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/tensorflow_probability/python/__init__.py:57: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead.
  if (distutils.version.LooseVersion(tf.__version__) <
/usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/tensorflow_probability/python/__init__.py:58: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead.
  distutils.version.LooseVersion(required_tensorflow_version)):

Welcome to PyDDA 2.0.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!

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
# Set initialization and do retrieval
grid1 = pydda.initialization.make_constant_wind_field(grid1, vel_field="VRADH_corr")
---> 17 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))

File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1474, in get_dd_wind_field(Grids, u_init, v_init, w_init, engine, **kwargs)
   w_init = new_grids[0]["w"].values.squeeze()
if (
   engine.lower() == "scipy"
   or engine.lower() == "jax"
   or engine.lower() == "auglag"
):
-> 1474     return _get_dd_wind_field_scipy(
       new_grids, u_init, v_init, w_init, engine, **kwargs
   )
elif engine.lower() == "tensorflow":
   return _get_dd_wind_field_tensorflow(
       new_grids, u_init, v_init, w_init, **kwargs
   )

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)
parameters.print_out = False
if engine.lower() == "scipy":
--> 602     winds = fmin_l_bfgs_b(
       J_function,
       winds,
       args=(parameters,),
       maxiter=max_iterations,
       pgtol=tolerance,
       bounds=bounds,
       fprime=grad_J,
       disp=0,
       iprint=-1,
       callback=_vert_velocity_callback,
   )
else:
   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)
callback = _wrap_callback(callback)
opts = {'disp': disp,
       'iprint': iprint,
       'maxcor': m,
   (...)
       'callback': callback,
       'maxls': maxls}
--> 237 res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
                      **opts)
d = {'grad': res['jac'],
    'task': res['message'],
    'funcalls': res['nfev'],
    'nit': res['nit'],
    'warnflag': res['status']}
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)
task_str = task.tobytes()
if task_str.startswith(b'FG'):
   # The minimization routine wants f and g at the current x.
   # Note that interruptions due to maxfun are postponed
   # until the completion of the current minimization iteration.
   # Overwrite f and g:
--> 407     f, g = func_and_grad(x)
elif task_str.startswith(b'NEW_X'):
   # new iteration
   n_iterations += 1

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

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:267, in ScalarFunction._update_grad(self)
def _update_grad(self):
   if not self.g_updated:
--> 267         self._update_grad_impl()
       self.g_updated = True

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:175, in ScalarFunction.__init__.<locals>.update_grad()
def update_grad():
--> 175     self.g = grad_wrapped(self.x)

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:172, in ScalarFunction.__init__.<locals>.grad_wrapped(x)
def grad_wrapped(x):
   self.ngev += 1
--> 172     return np.atleast_1d(grad(np.copy(x), *args))

File ~/work/PyDDA/PyDDA/pydda/cost_functions/cost_functions.py:528, in grad_J(winds, parameters)
grad = _cost_functions_numpy.calculate_grad_radial_vel(
   parameters.vrs,
   parameters.els,
   (...)
   upper_bc=parameters.upper_bc,
)
if parameters.Cm > 0:
--> 528     grad += _cost_functions_numpy.calculate_mass_continuity_gradient(
       winds[0],
       winds[1],
       winds[2],
       parameters.z,
       parameters.dx,
       parameters.dy,
       parameters.dz,
       coeff=parameters.Cm,
       upper_bc=parameters.upper_bc,
   )
if parameters.Cx > 0 or parameters.Cy > 0 or parameters.Cz > 0:
   grad += _cost_functions_numpy.calculate_smoothness_gradient(
       winds[0],
       winds[1],
   (...)
       upper_bc=parameters.upper_bc,
   )

File ~/work/PyDDA/PyDDA/pydda/cost_functions/_cost_functions_numpy.py:465, in calculate_mass_continuity_gradient(u, v, w, z, dx, dy, dz, coeff, anel, upper_bc)
grad_u = -np.gradient(div, dx, axis=2) * coeff
grad_v = -np.gradient(div, dy, axis=1) * coeff
--> 465 grad_w = -np.gradient(div, dz, axis=0) * coeff
# Impermeability condition
grad_w[0, :, :] = 0

File <__array_function__ internals>:180, in gradient(*args, **kwargs)

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/numpy/lib/function_base.py:1256, in gradient(f, axis, edge_order, *varargs)
dx_0 = ax_dx if uniform_spacing else ax_dx[0]
# 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
-> 1256 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
slice1[axis] = -1
slice2[axis] = -1

KeyboardInterrupt: 