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=2,
wind_tol=0.5,
max_iterations=150,
engine="scipy",
parallel=False,
)
# 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 Office of Science as part of
## the Atmospheric Radiation Measurement (ARM) 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.5.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...
Unable to load both TensorFlow and tensorflow-probability. TensorFlow engine disabled.
No module named 'tensorflow'
False
Calculating weights for radars 0 and 1
Calculating weights for radars 1 and 0
Calculating weights for models...
Starting solver
rmsVR = 6.827303971100176
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.660301869612016
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
10| 1.8730| 41.7898| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 11.3127
Max change in w: 10.571
The gradient of the cost functions is 0.1395435664190129
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
20| 0.3720| 19.6893| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 11.9101
Max change in w: 6.110
The gradient of the cost functions is 0.10614552703397478
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
30| 0.1939| 11.4887| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 11.8181
Max change in w: 3.992
The gradient of the cost functions is 0.08095070591641679
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
40| 0.2357| 7.6193| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 14.0516
Max change in w: 4.072
The gradient of the cost functions is 0.06602647187439353
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
50| 0.0689| 5.3911| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 18.7420
Max change in w: 5.336
The gradient of the cost functions is 0.03798101974532811
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
60| 0.0451| 4.0065| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 22.0508
Max change in w: 3.816
The gradient of the cost functions is 0.03507602692863761
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
70| 0.0505| 3.1445| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 24.9823
Max change in w: 2.847
The gradient of the cost functions is 0.03643970677079748
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
80| 0.0820| 2.4891| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 27.8134
Max change in w: 2.942
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[1], line 14
10 )
11
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,
File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1746, in get_dd_wind_field(Grids, u_init, v_init, w_init, engine, **kwargs)
1739 w_init = new_grids[0]["w"].values.squeeze()
1741 if (
1742 engine.lower() == "scipy"
1743 or engine.lower() == "jax"
1744 or engine.lower() == "auglag"
1745 ):
-> 1746 return _get_dd_wind_field_scipy(
1747 new_grids, u_init, v_init, w_init, engine, **kwargs
1748 )
1749 elif engine.lower() == "tensorflow":
1750 if kwargs.get("Cterrain", 0.0) > 0:
File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:788, 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, Cterrain, terrain, Cvad, cvtol, gtol, Jveltol, Ut, Vt, low_pass_filter, mask_outside_opt, weights_obs, weights_model, weights_bg, max_iterations, mask_w_outside_opt, filter_type, filter_window, filter_order, leise_nstep, min_bca, max_bca, upper_bc, lower_bc, above, model_fields, output_cost_functions, roi, wind_tol, tolerance, const_boundary_cond, max_wind_mag, parallel)
786 parameters.print_out = False
787 if engine.lower() == "scipy":
--> 788 winds = fmin_l_bfgs_b(
789 J_and_grad,
790 winds,
791 args=(parameters,),
792 maxiter=max_iterations,
793 pgtol=tolerance,
794 bounds=bounds,
795 callback=_vert_velocity_callback,
796 )
797 else:
799 def loss_and_gradient(x):
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_lbfgsb_py.py:259, in fmin_l_bfgs_b(func, x0, fprime, args, approx_grad, bounds, m, factr, pgtol, epsilon, maxfun, maxiter, callback, maxls)
249 callback = _wrap_callback(callback)
250 opts = {'maxcor': m,
251 'ftol': factr * np.finfo(float).eps,
252 'gtol': pgtol,
(...) 256 'callback': callback,
257 'maxls': maxls}
--> 259 res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
260 **opts)
261 d = {'grad': res['jac'],
262 'task': res['message'],
263 'funcalls': res['nfev'],
264 'nit': res['nit'],
265 'warnflag': res['status']}
266 f = res['fun']
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_lbfgsb_py.py:412, in _minimize_lbfgsb(fun, x0, args, jac, bounds, maxcor, ftol, gtol, eps, maxfun, maxiter, callback, maxls, finite_diff_rel_step, workers, **unknown_options)
410 g = g.astype(np.float64)
411 # x, f, g, wa, iwa, task, csave, lsave, isave, dsave = \
--> 412 _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr, pgtol, wa,
413 iwa, task, lsave, isave, dsave, maxls, ln_task)
415 if task[0] == 3:
416 # The minimization routine wants f and g at the current x.
417 # Note that interruptions due to maxfun are postponed
418 # until the completion of the current minimization iteration.
419 # Overwrite f and g:
420 f, g = func_and_grad(x)
KeyboardInterrupt: