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
Welcome to PyDDA 2.4.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...
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.1119091690376132
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
10| 92.4335| 17.0200| 0.0001| 0.0000| 0.0000| 0.0000| 0.0000| 48.1870
Max change in w: nan
---------------------------------------------------------------------------
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:1502, in get_dd_wind_field(Grids, u_init, v_init, w_init, engine, **kwargs)
1495 w_init = new_grids[0]["w"].values.squeeze()
1497 if (
1498 engine.lower() == "scipy"
1499 or engine.lower() == "jax"
1500 or engine.lower() == "auglag"
1501 ):
-> 1502 return _get_dd_wind_field_scipy(
1503 new_grids, u_init, v_init, w_init, engine, **kwargs
1504 )
1505 elif engine.lower() == "tensorflow":
1506 return _get_dd_wind_field_tensorflow(
1507 new_grids, u_init, v_init, w_init, **kwargs
1508 )
File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:614, 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, parallel)
612 parameters.print_out = False
613 if engine.lower() == "scipy":
--> 614 winds = fmin_l_bfgs_b(
615 J_function,
616 winds,
617 args=(parameters,),
618 maxiter=max_iterations,
619 pgtol=tolerance,
620 bounds=bounds,
621 fprime=grad_J,
622 callback=_vert_velocity_callback,
623 )
624 else:
626 def loss_and_gradient(x):
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_lbfgsb_py.py:281, in fmin_l_bfgs_b(func, x0, fprime, args, approx_grad, bounds, m, factr, pgtol, epsilon, iprint, maxfun, maxiter, disp, callback, maxls)
269 callback = _wrap_callback(callback)
270 opts = {'disp': disp,
271 'iprint': iprint,
272 'maxcor': m,
(...) 278 'callback': callback,
279 'maxls': maxls}
--> 281 res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
282 **opts)
283 d = {'grad': res['jac'],
284 'task': res['message'],
285 'funcalls': res['nfev'],
286 'nit': res['nit'],
287 'warnflag': res['status']}
288 f = res['fun']
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_lbfgsb_py.py:469, in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, workers, **unknown_options)
461 _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr, pgtol, wa,
462 iwa, task, lsave, isave, dsave, maxls, ln_task)
464 if task[0] == 3:
465 # The minimization routine wants f and g at the current x.
466 # Note that interruptions due to maxfun are postponed
467 # until the completion of the current minimization iteration.
468 # Overwrite f and g:
--> 469 f, g = func_and_grad(x)
470 elif task[0] == 1:
471 # new iteration
472 n_iterations += 1
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:413, in ScalarFunction.fun_and_grad(self, x)
411 self._update_x(x)
412 self._update_fun()
--> 413 self._update_grad()
414 return self.f, self.g
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:375, in ScalarFunction._update_grad(self)
373 if self._orig_grad in FD_METHODS:
374 self._update_fun()
--> 375 self.g = self._wrapped_grad(self.x, f0=self.f)
376 self.g_updated = True
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:39, in _ScalarGradWrapper.__call__(self, x, f0, **kwds)
35 def __call__(self, x, f0=None, **kwds):
36 # Send a copy because the user may overwrite it.
37 # The user of this class might want `x` to remain unchanged.
38 if callable(self.grad):
---> 39 g = np.atleast_1d(self.grad(np.copy(x), *self.args))
40 elif self.grad in FD_METHODS:
41 g, dct = approx_derivative(
42 self.fun,
43 x,
44 f0=f0,
45 **self.finite_diff_options,
46 )
File ~/work/PyDDA/PyDDA/pydda/cost_functions/cost_functions.py:517, in grad_J(winds, parameters)
515 if parameters.parallel:
516 futures = []
--> 517 with ThreadPoolExecutor() as pool:
518 futures.append(
519 pool.submit(
520 _cost_functions_numpy.calculate_grad_radial_vel,
(...) 533 )
534 )
535 if parameters.Cm > 0:
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/concurrent/futures/_base.py:667, in Executor.__exit__(self, exc_type, exc_val, exc_tb)
666 def __exit__(self, exc_type, exc_val, exc_tb):
--> 667 self.shutdown(wait=True)
668 return False
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/concurrent/futures/thread.py:273, in ThreadPoolExecutor.shutdown(self, wait, cancel_futures)
271 if wait:
272 for t in self._threads:
--> 273 t.join()
File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/threading.py:1133, in Thread.join(self, timeout)
1130 if timeout is not None:
1131 timeout = max(timeout, 0)
-> 1133 self._os_thread_handle.join(timeout)
KeyboardInterrupt: