Example on integrating radar and HRRR data#
This is an example of how to retrieve winds in Hurricane Florence. In this example, we use data from 2 NEXRAD radars as well as from the HRRR to retrieve the winds
This example has been updated to use _Herbie to retrieve the HRRR data. In addition, _pooch is used to retrieve the gridded data for the example. Herbie is not required to run PyDDA, but must be installed to run this example.
Author: Robert C. Jackson
import pydda
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
from herbie import Herbie
H = Herbie("2018-09-14 06:00", model="hrrr", product="prs", fxx=0)
H.download()
grid_mhx_path = pydda.tests.get_sample_file("grid_mhx.nc")
grid_ltx_path = pydda.tests.get_sample_file("grid_ltx.nc")
grid_mhx = pydda.io.read_grid(grid_mhx_path)
grid_ltx = pydda.io.read_grid(grid_ltx_path)
grid_mhx = pydda.constraints.add_hrrr_constraint_to_grid(grid_mhx, H.grib)
grid_mhx = pydda.initialization.make_constant_wind_field(grid_mhx, (0.0, 0.0, 0.0))
out_grids, _ = pydda.retrieval.get_dd_wind_field(
[grid_mhx, grid_ltx],
Co=1e-2,
Cm=128.0,
Cmod=1e-4,
Cx=1e-4,
Cy=1e-4,
Cz=1e-4,
max_iterations=100,
mask_outside_opt=True,
vel_name="corrected_velocity",
engine="scipy",
model_fields=["hrrr"],
)
fig = plt.figure(figsize=(25, 15))
ax = plt.axes(projection=ccrs.PlateCarree())
ax = pydda.vis.plot_horiz_xsection_barbs_map(
out_grids,
ax=ax,
bg_grid_no=-1,
level=3,
barb_spacing_x_km=20.0,
barb_spacing_y_km=20.0,
cmap="ChaseSpectral",
)
ax.set_xticks(np.arange(-80, -75, 0.5))
ax.set_yticks(np.arange(33.0, 35.5, 0.5))
plt.title(out_grids[0].time.attrs["units"][13:] + " winds at 0.5 km")
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.2
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'
✅ Found ┊ model=hrrr ┊ product=prs ┊ 2018-Sep-14 06:00 UTC F00 ┊ GRIB2 @ local ┊ IDX @ aws
False
Calculating weights for radars 0 and 1
Calculating weights for radars 1 and 0
Calculating weights for models...
Starting solver
rmsVR = 21.58818045293143
Total points: 1284618
The max of w_init is 0.0
Total number of model points: 4363281
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[1], line 19
17 grid_mhx = pydda.constraints.add_hrrr_constraint_to_grid(grid_mhx, H.grib)
18 grid_mhx = pydda.initialization.make_constant_wind_field(grid_mhx, (0.0, 0.0, 0.0))
---> 19 out_grids, _ = pydda.retrieval.get_dd_wind_field(
20 [grid_mhx, grid_ltx],
21 Co=1e-2,
22 Cm=128.0,
23 Cmod=1e-4,
24 Cx=1e-4,
25 Cy=1e-4,
26 Cz=1e-4,
27 max_iterations=100,
28 mask_outside_opt=True,
29 vel_name="corrected_velocity",
30 engine="scipy",
31 model_fields=["hrrr"],
32 )
34 fig = plt.figure(figsize=(25, 15))
35 ax = plt.axes(projection=ccrs.PlateCarree())
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:400, 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)
398 raise ValueError('length of x0 != length of bounds')
399 else:
--> 400 bounds = np.array(old_bound_to_new(bounds))
402 # check bounds
403 if (bounds[0] > bounds[1]).any():
File /usr/share/miniconda/envs/pydda-docs/lib/python3.12/site-packages/scipy/optimize/_constraints.py:441, in old_bound_to_new(bounds)
437 lb, ub = zip(*bounds)
439 # Convert occurrences of None to -inf or inf, and replace occurrences of
440 # any numpy array x with x.item(). Then wrap the results in numpy arrays.
--> 441 lb = np.array([float(_arr_to_scalar(x)) if x is not None else -np.inf
442 for x in lb])
443 ub = np.array([float(_arr_to_scalar(x)) if x is not None else np.inf
444 for x in ub])
446 return lb, ub
KeyboardInterrupt: