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

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'

✅ 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
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))
---> 19 out_grids, _ = pydda.retrieval.get_dd_wind_field(
   [grid_mhx, grid_ltx],
   Co=1e-2,
   Cm=128.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)
   w_init = new_grids[0]["w"].values.squeeze()
if (
   engine.lower() == "scipy"
   or engine.lower() == "jax"
   or engine.lower() == "auglag"
):
-> 1502     return _get_dd_wind_field_scipy(
       new_grids, u_init, v_init, w_init, engine, **kwargs
   )
elif engine.lower() == "tensorflow":
   return _get_dd_wind_field_tensorflow(
       new_grids, u_init, v_init, w_init, **kwargs
   )

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)
parameters.print_out = False
if engine.lower() == "scipy":
--> 614     winds = fmin_l_bfgs_b(
       J_function,
       winds,
       args=(parameters,),
       maxiter=max_iterations,
       pgtol=tolerance,
       bounds=bounds,
       fprime=grad_J,
       callback=_vert_velocity_callback,
   )
else:
   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)
callback = _wrap_callback(callback)
opts = {'disp': disp,
       'iprint': iprint,
       'maxcor': m,
   (...)    278         'callback': callback,
       'maxls': maxls}
--> 281 res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
                      **opts)
d = {'grad': res['jac'],
    'task': res['message'],
    'funcalls': res['nfev'],
    'nit': res['nit'],
    'warnflag': res['status']}
f = res['fun']

File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/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)
   raise ValueError('length of x0 != length of bounds')
else:
--> 400     bounds = np.array(old_bound_to_new(bounds))
   # check bounds
   if (bounds[0] > bounds[1]).any():

File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_constraints.py:450, in old_bound_to_new(bounds)
lb, ub = zip(*bounds)
# Convert occurrences of None to -inf or inf, and replace occurrences of
# any numpy array x with x.item(). Then wrap the results in numpy arrays.
--> 450 lb = np.array([float(_arr_to_scalar(x)) if x is not None else -np.inf
              for x in lb])
ub = np.array([float(_arr_to_scalar(x)) if x is not None else np.inf
              for x in ub])
return lb, ub

File /usr/share/miniconda/envs/pydda-docs/lib/python3.14/site-packages/scipy/optimize/_constraints.py:16, in _arr_to_scalar(x)
from ._optimize import OptimizeWarning
from scipy._lib._sparse import issparse
---> 16 def _arr_to_scalar(x):
   # If x is a numpy array, return x.item().  This will
   # fail if the array has more than one element.
   return x.item() if isinstance(x, np.ndarray) else x
class NonlinearConstraint:

KeyboardInterrupt: 