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.1.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...
TensorFlow detected. Checking for tensorflow-probability...
TensorFlow-probability detected. TensorFlow engine enabled!
✅ 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.58818045293141
Total points: 1284618
The max of w_init is 0.0

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[1], line 19
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,
   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())

File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:1474, 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"
):
-> 1474     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:527, 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)
print("The max of w_init is", wcurrmax)
iterations = 0
--> 527 bounds = [(-x, x) for x in max_wind_mag * np.ones(winds.shape)]
if model_fields is not None:
   for i, the_field in enumerate(model_fields):

File ~/work/PyDDA/PyDDA/pydda/retrieval/wind_retrieve.py:527, in <listcomp>(.0)
print("The max of w_init is", wcurrmax)
iterations = 0
--> 527 bounds = [(-x, x) for x in max_wind_mag * np.ones(winds.shape)]
if model_fields is not None:
   for i, the_field in enumerate(model_fields):

KeyboardInterrupt: 