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
Failed to import TF-Keras. Please note that TF-Keras is not installed by default when you install TensorFlow Probability. This is so that JAX-only users do not have to install TensorFlow or TF-Keras. To use TensorFlow Probability with TensorFlow, please install the tf-keras or tf-keras-nightly package.
This can be be done through installing the tensorflow-probability[tf] extra.
Welcome to PyDDA 2.1.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...
TensorFlow detected. Checking for tensorflow-probability...
Failed to import TF-Keras. Please note that TF-Keras is not installed by default when you install TensorFlow Probability. This is so that JAX-only users do not have to install TensorFlow or TF-Keras. To use TensorFlow Probability with TensorFlow, please install the tf-keras or tf-keras-nightly package.
This can be be done through installing the tensorflow-probability[tf] extra.
Unable to load both TensorFlow and tensorflow-probability. TensorFlow engine disabled.
No module named 'tf_keras'
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.242250167336594
Total points: 281276
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|2864.7900| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000
The gradient of the cost functions is 0.12289888187724508
Nfeval | Jvel | Jmass | Jsmooth | Jbg | Jvort | Jmodel | Jpoint | Max w
10| 117.2266| 23.2756| 0.0000| 0.0000| 0.0000| 0.0000| 0.0000| 18.4482
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[1], line 17
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,
21 Cx=10,
22 Cy=10,
23 Cz=10,
24 vel_name="VRADH_corr",
25 refl_field="DBZH",
26 mask_outside_opt=True,
27 wind_tol=0.5,
28 max_iterations=200,
29 engine="scipy",
30 )
31 # Make a neat plot
32 fig = plt.figure(figsize=(10, 7))
File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1472, in get_dd_wind_field(Grids, u_init, v_init, w_init, engine, **kwargs)
1465 w_init = new_grids[0]["w"].values.squeeze()
1467 if (
1468 engine.lower() == "scipy"
1469 or engine.lower() == "jax"
1470 or engine.lower() == "auglag"
1471 ):
-> 1472 return _get_dd_wind_field_scipy(
1473 new_grids, u_init, v_init, w_init, engine, **kwargs
1474 )
1475 elif engine.lower() == "tensorflow":
1476 return _get_dd_wind_field_tensorflow(
1477 new_grids, u_init, v_init, w_init, **kwargs
1478 )
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 callback=_vert_velocity_callback,
611 )
612 else:
614 def loss_and_gradient(x):
File /usr/share/miniconda/envs/pydda-docs/lib/python3.12/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.12/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.12/site-packages/scipy/optimize/_differentiable_functions.py:404, in ScalarFunction.fun_and_grad(self, x)
402 self._update_x(x)
403 self._update_fun()
--> 404 self._update_grad()
405 return self.f, self.g
File /usr/share/miniconda/envs/pydda-docs/lib/python3.12/site-packages/scipy/optimize/_differentiable_functions.py:366, in ScalarFunction._update_grad(self)
364 if self._orig_grad in FD_METHODS:
365 self._update_fun()
--> 366 self.g = self._wrapped_grad(self.x, f0=self.f)
367 self.g_updated = True
File /usr/share/miniconda/envs/pydda-docs/lib/python3.12/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: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:128, in calculate_grad_radial_vel(vrs, els, azs, u, v, w, wts, weights, rmsVr, coeff, upper_bc)
120 for i in range(len(vrs)):
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
KeyboardInterrupt: