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.3
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 1 and 0
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
---------------------------------------------------------------------------
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: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:513, in grad_J(winds, parameters)
503 elif parameters.engine == "scipy":
504 winds = np.reshape(
505 winds,
506 (
(...)
511 ),
512 )
--> 513 grad = _cost_functions_numpy.calculate_grad_radial_vel(
514 parameters.vrs,
515 parameters.els,
516 parameters.azs,
517 winds[0],
518 winds[1],
519 winds[2],
520 parameters.wts,
521 parameters.weights,
522 parameters.rmsVr,
523 coeff=parameters.Co,
524 upper_bc=parameters.upper_bc,
525 )
527 if parameters.Cm > 0:
528 grad += _cost_functions_numpy.calculate_mass_continuity_gradient(
529 winds[0],
530 winds[1],
(...)
537 upper_bc=parameters.upper_bc,
538 )
File ~/work/PyDDA/PyDDA/pydda/cost_functions/_cost_functions_numpy.py:131, in calculate_grad_radial_vel(vrs, els, azs, u, v, w, wts, weights, rmsVr, coeff, upper_bc)
121 v_ar = (
122 np.cos(els[i]) * np.sin(azs[i]) * u
123 + np.cos(els[i]) * np.cos(azs[i]) * v
124 + np.sin(els[i]) * (w - np.abs(wts[i]))
125 )
127 x_grad = (
128 2 * (v_ar - vrs[i]) * np.cos(els[i]) * np.sin(azs[i]) * weights[i]
129 ) * lambda_o
130 y_grad = (
--> 131 2 * (v_ar - vrs[i]) * np.cos(els[i]) * np.cos(azs[i]) * weights[i]
132 ) * lambda_o
133 z_grad = (2 * (v_ar - vrs[i]) * np.sin(els[i]) * weights[i]) * lambda_o
135 p_x1 += x_grad
KeyboardInterrupt: