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 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 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.359750110231316
Total points: 399749
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|4080.4252| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000
The gradient of the cost functions is 0.12208368625646192
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
10| 91.8283| 22.9782| 0.0001| 0.0000| 0.0000| 0.0000| 0.0000| 80.8833
Max change in w: nan
The gradient of the cost functions is 0.11469278158160577
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
20| 85.8892| 20.6679| 0.0001| 0.0000| 0.0000| 0.0000| 0.0000| 94.1776
The gradient of the cost functions is 0.07331884235807853
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
30| 84.9256| 21.3641| 0.0001| 0.0000| 0.0000| 0.0000| 0.0000| 96.2424
The gradient of the cost functions is 0.13249227771053818
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
40| 84.3181| 22.0388| 0.0001| 0.0000| 0.0000| 0.0000| 0.0000| 97.0717
The gradient of the cost functions is 0.0673313941252296
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
50| 84.1381| 21.9409| 0.0001| 0.0000| 0.0000| 0.0000| 0.0000| 97.1398
Max change in w: nan
The gradient of the cost functions is 0.0656208074782532
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
60| 83.8890| 22.0606| 0.0001| 0.0000| 0.0000| 0.0000| 0.0000| 97.3290
The gradient of the cost functions is 0.06240248465031262
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
70| 83.4992| 22.4180| 0.0001| 0.0000| 0.0000| 0.0000| 0.0000| 97.6397
The gradient of the cost functions is 0.06205928785181723
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[1], line 17
13 grid4 = pydda.io.read_grid(grid4_path)
14
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,
File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1745, in get_dd_wind_field(Grids, u_init, v_init, w_init, engine, **kwargs)
1738 w_init = new_grids[0]["w"].values.squeeze()
1740 if (
1741 engine.lower() == "scipy"
1742 or engine.lower() == "jax"
1743 or engine.lower() == "auglag"
1744 ):
-> 1745 return _get_dd_wind_field_scipy(
1746 new_grids, u_init, v_init, w_init, engine, **kwargs
1747 )
1748 elif engine.lower() == "tensorflow":
1749 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 ~/miniconda3/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 ~/miniconda3/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_lbfgsb_py.py:420, in _minimize_lbfgsb(fun, x0, args, jac, bounds, maxcor, ftol, gtol, eps, maxfun, maxiter, callback, maxls, finite_diff_rel_step, workers, **unknown_options)
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)
421 elif task[0] == 1:
422 # new iteration
423 n_iterations += 1
File ~/miniconda3/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:412, in ScalarFunction.fun_and_grad(self, x)
410 if not np.array_equal(x, self.x):
411 self._update_x(x)
--> 412 self._update_fun()
413 self._update_grad()
414 return self.f, self.g
File ~/miniconda3/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:362, in ScalarFunction._update_fun(self)
360 def _update_fun(self):
361 if not self.f_updated:
--> 362 fx = self._wrapped_fun(self.x)
363 self._nfev += 1
364 if fx < self._lowest_f:
File ~/miniconda3/envs/pydda-docs/lib/python3.14/site-packages/scipy/_lib/_util.py:544, in _ScalarFunctionWrapper.__call__(self, x)
541 def __call__(self, x):
542 # Send a copy because the user may overwrite it.
543 # The user of this class might want `x` to remain unchanged.
--> 544 fx = self.f(np.copy(x), *self.args)
545 self.nfev += 1
547 # Make sure the function returns a true scalar
File ~/miniconda3/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_optimize.py:80, in MemoizeJac.__call__(self, x, *args)
78 def __call__(self, x, *args):
79 """ returns the function value """
---> 80 self._compute_if_needed(x, *args)
81 return self._value
File ~/miniconda3/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_optimize.py:74, in MemoizeJac._compute_if_needed(self, x, *args)
72 if not np.all(x == self.x) or self._value is None or self.jac is None:
73 self.x = np.asarray(x).copy()
---> 74 fg = self.fun(x, *args)
75 self.jac = fg[1]
76 self._value = fg[0]
File ~/work/PyDDA/PyDDA/pydda/cost_functions/cost_functions.py:1009, in J_and_grad(winds, parameters)
997 return J_function(winds, parameters), grad_J(winds, parameters)
999 shaped_winds = np.reshape(
1000 winds,
1001 (
(...) 1006 ),
1007 )
1008 radial_cost, radial_gradient = (
-> 1009 _cost_functions_numpy.calculate_radial_vel_cost_and_gradient(
1010 parameters.vrs,
1011 parameters.els,
1012 parameters.azs,
1013 shaped_winds[0],
1014 shaped_winds[1],
1015 shaped_winds[2],
1016 parameters.wts,
1017 parameters.weights,
1018 parameters.rmsVr,
1019 coeff=parameters.Co,
1020 upper_bc=parameters.upper_bc,
1021 upper_bc_mask=parameters.upper_bc_mask,
1022 parallel=False,
1023 radar_geometry=parameters.radar_geometry,
1024 )
1025 )
1026 previous_cache = getattr(parameters, "_radial_eval_cache", None)
1027 parameters._radial_eval_cache = {
1028 "source_winds": winds,
1029 "cost": radial_cost,
1030 "gradient": radial_gradient,
1031 }
File ~/work/PyDDA/PyDDA/pydda/cost_functions/_cost_functions_numpy.py:312, in calculate_radial_vel_cost_and_gradient(vrs, els, azs, u, v, w, wts, weights, rmsVr, coeff, upper_bc, upper_bc_mask, parallel, radar_geometry)
310 p_x1 += 2 * residual * x_projection * weight * lambda_o
311 p_y1 += 2 * residual * y_projection * weight * lambda_o
--> 312 p_z1 += 2 * residual * sin_el * weight * lambda_o
314 p_z1[0, :, :] = 0
315 p_z1 = _apply_upper_bc(p_z1, upper_bc, upper_bc_mask)
KeyboardInterrupt: