Example on retrieving and plotting winds#

This is a simple example for how to retrieve and plot winds from 2 radars using PyDDA.

Author: Robert C. Jackson

import pydda
from matplotlib import pyplot as plt

berr_grid = pydda.io.read_grid(pydda.tests.EXAMPLE_RADAR0)
cpol_grid = pydda.io.read_grid(pydda.tests.EXAMPLE_RADAR1)

# Load sounding data and insert as an intialization
berr_grid = pydda.initialization.make_constant_wind_field(
    berr_grid, (0.0, 0.0, 0.0), vel_field="corrected_velocity"
)

# Start the wind retrieval. This example only uses the mass continuity
# and data weighting constraints.
Grids, _ = pydda.retrieval.get_dd_wind_field(
    [berr_grid, cpol_grid],
    Co=1.0,
    Cm=256.0,
    Cx=0.0,
    Cy=0.0,
    Cz=0.0,
    Cb=0.0,
    frz=5000.0,
    filter_window=5,
    mask_outside_opt=True,
    upper_bc=1,
    wind_tol=0.5,
    engine="scipy",
)

# Plot a horizontal cross section
plt.figure(figsize=(9, 9))
pydda.vis.plot_horiz_xsection_barbs(
    Grids,
    background_field="reflectivity",
    level=6,
    w_vel_contours=[5, 10, 15],
    barb_spacing_x_km=5.0,
    barb_spacing_y_km=15.0,
    vmin=0,
    vmax=70,
)
plt.show()

# Plot a vertical X-Z cross section
plt.figure(figsize=(9, 9))
pydda.vis.plot_xz_xsection_barbs(
    Grids,
    background_field="reflectivity",
    level=40,
    w_vel_contours=[5, 10, 15],
    barb_spacing_x_km=10.0,
    barb_spacing_z_km=2.0,
    vmin=0,
    vmax=70,
)
plt.show()

# Plot a vertical Y-Z cross section
plt.figure(figsize=(9, 9))
pydda.vis.plot_yz_xsection_barbs(
    Grids,
    background_field="reflectivity",
    level=40,
    barb_spacing_y_km=10.0,
    barb_spacing_z_km=2.0,
    vmin=0,
    vmax=70,
)
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.0.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!
/home/runner/work/PyDDA/PyDDA/pydda/retrieval/angles.py:25: RuntimeWarning: invalid value encountered in arccos
  elev = np.arccos((Re**2 + slantrsq - rh**2) / (2 * Re * slantr))
False
Calculating weights for radars 0 and 1
Calculating weights for radars 1 and 0
/home/runner/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1523: RuntimeWarning: invalid value encountered in divide
  inp_array1 = x / a
/home/runner/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1529: RuntimeWarning: invalid value encountered in divide
  inp_array3 = (a * a + b * b - c * c) / (2 * a * b)
/home/runner/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1526: RuntimeWarning: divide by zero encountered in divide
  inp_array2 = (x - rad2[1]) / b
/home/runner/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1529: RuntimeWarning: invalid value encountered in divide
  inp_array3 = (a * a + b * b - c * c) / (2 * a * b)
Calculating weights for models...
Starting solver 
rmsVR = 6.827303971100175
Total points: 81194
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|83859.8222|   0.0000|   0.0000|   0.0000|   0.0000|   0.0000|   0.0000|   0.0000
The gradient of the cost functions is 0.6631357968996711
Nfeval | Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Jpoint  | Max w  
     10|   1.8782|  41.0864|   0.0000|   0.0000|   0.0000|   0.0000|   0.0000|  11.3085
Max change in w: 10.561
The gradient of the cost functions is 0.13750649771706588
Nfeval | Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Jpoint  | Max w  
     20|   0.3524|  18.6065|   0.0000|   0.0000|   0.0000|   0.0000|   0.0000|  11.8823
Max change in w: 6.008
The gradient of the cost functions is 0.09828068193831543
Nfeval | Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Jpoint  | Max w  
     30|   0.1953|  10.8039|   0.0000|   0.0000|   0.0000|   0.0000|   0.0000|  12.4880
Max change in w: 4.003
The gradient of the cost functions is 0.09730430370400128
Nfeval | Jvel    | Jmass   | Jsmooth |   Jbg   | Jvort   | Jmodel  | Jpoint  | Max w  
     40|   0.1245|   7.0774|   0.0000|   0.0000|   0.0000|   0.0000|   0.0000|  16.4854
Max change in w: 5.381
The gradient of the cost functions is 0.05299677215060275
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[1], line 14
      8 berr_grid = pydda.initialization.make_constant_wind_field(
      9     berr_grid, (0.0, 0.0, 0.0), vel_field="corrected_velocity"
     10 )
     12 # Start the wind retrieval. This example only uses the mass continuity
     13 # and data weighting constraints.
---> 14 Grids, _ = pydda.retrieval.get_dd_wind_field(
     15     [berr_grid, cpol_grid],
     16     Co=1.0,
     17     Cm=256.0,
     18     Cx=0.0,
     19     Cy=0.0,
     20     Cz=0.0,
     21     Cb=0.0,
     22     frz=5000.0,
     23     filter_window=5,
     24     mask_outside_opt=True,
     25     upper_bc=1,
     26     wind_tol=0.5,
     27     engine="scipy",
     28 )
     30 # Plot a horizontal cross section
     31 plt.figure(figsize=(9, 9))

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:199, in fmin_l_bfgs_b(func, x0, fprime, args, approx_grad, bounds, m, factr, pgtol, epsilon, iprint, maxfun, maxiter, disp, callback, maxls)
    187 callback = _wrap_callback(callback)
    188 opts = {'disp': disp,
    189         'iprint': iprint,
    190         'maxcor': m,
   (...)
    196         'callback': callback,
    197         'maxls': maxls}
--> 199 res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
    200                        **opts)
    201 d = {'grad': res['jac'],
    202      'task': res['message'],
    203      'funcalls': res['nfev'],
    204      'nit': res['nit'],
    205      'warnflag': res['status']}
    206 f = res['fun']

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

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:285, in ScalarFunction.fun_and_grad(self, x)
    283 if not np.array_equal(x, self.x):
    284     self._update_x_impl(x)
--> 285 self._update_fun()
    286 self._update_grad()
    287 return self.f, self.g

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:251, in ScalarFunction._update_fun(self)
    249 def _update_fun(self):
    250     if not self.f_updated:
--> 251         self._update_fun_impl()
    252         self.f_updated = True

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:155, in ScalarFunction.__init__.<locals>.update_fun()
    154 def update_fun():
--> 155     self.f = fun_wrapped(self.x)

File /usr/share/miniconda3/envs/pydda-docs/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:137, in ScalarFunction.__init__.<locals>.fun_wrapped(x)
    133 self.nfev += 1
    134 # Send a copy because the user may overwrite it.
    135 # Overwriting results in undefined behaviour because
    136 # fun(self.x) will change self.x, with the two no longer linked.
--> 137 fx = fun(np.copy(x), *args)
    138 # Make sure the function returns a true scalar
    139 if not np.isscalar(fx):

File ~/work/PyDDA/PyDDA/pydda/cost_functions/cost_functions.py:178, in J_function(winds, parameters)
    168 winds = np.reshape(
    169     winds,
    170     (
   (...)
    175     ),
    176 )
    177 # Had to change to float because Jax returns device array (use np.float_())
--> 178 Jvel = _cost_functions_numpy.calculate_radial_vel_cost_function(
    179     parameters.vrs,
    180     parameters.azs,
    181     parameters.els,
    182     winds[0],
    183     winds[1],
    184     winds[2],
    185     parameters.wts,
    186     rmsVr=parameters.rmsVr,
    187     weights=parameters.weights,
    188     coeff=parameters.Co,
    189 )
    190 # print("apples Jvel", Jvel)
    192 if parameters.Cm > 0:
    193     # Had to change to float because Jax returns device array (use np.float_())

File ~/work/PyDDA/PyDDA/pydda/cost_functions/_cost_functions_numpy.py:62, in calculate_radial_vel_cost_function(vrs, azs, els, u, v, w, wts, rmsVr, weights, coeff)
     58 lambda_o = coeff / (rmsVr * rmsVr)
     59 for i in range(len(vrs)):
     60     v_ar = (
     61         np.cos(els[i]) * np.sin(azs[i]) * u
---> 62         + np.cos(els[i]) * np.cos(azs[i]) * v
     63         + np.sin(els[i]) * (w - np.abs(wts[i]))
     64     )
     65     the_weight = weights[i]
     66     J_o += lambda_o * np.sum(np.square(vrs[i] - v_ar) * the_weight)

KeyboardInterrupt: