PyDDA is an open source Python package for retrieving winds using Multiple Doppler analyses. It uses the 3D variational technique for retrieving winds from multiple Doppler radars. It uses Py-ART Grid Objects as inputs. Therefore, preprocessed and gridded data are needed for PyDDA to run.
3D variational technique¶
PyDDA uses a 3D variational technique to retrieve the 3D wind field. We will leave students that are interested in more information about these technqiues two references to read at the end of this notebook. A basic introduction to 3D variational analysis is given here.
PyDDA minimizes a cost function that corresponds to various penalties including:
which corresponds to the mass continuity equation.
RMSE between radar winds and analysis winds.
RMSE between sounding winds and analysis winds.
which corresponds to the smoothness of wind field to eliminate high-frequency noise that can result from numerical instability.
The cost function to be minimized is a weighted sum of the various cost functions in PyDDA and are represented in Equation (1):
$J = c_{m}J_{m} + c_{o}J_{o} + c_{b}J_{b} + c_{s}J_{s} + ...$ (1)
References¶
Imports¶
Let’s import the necessary libraries. For now, we’ll need PyART, glob, matplotlib, and PyDDA.
import glob
import pyart
import matplotlib.pyplot as plt
import pydda
import warnings
import numpy as np
import pandas as pd
import fsspec
import cartopy.crs as ccrs
warnings.filterwarnings("ignore")Case study¶
We will examine the same case study over Italy that has been used for this entire short course. For this case in May 2024, we had weak convective cells in the region of
For this case, we had coverage of the storm from two radars, a MeteoSwiss C-band radar near Monte Lema, CH and an X-band radar from ARPA Lombardia located about 50 km away in Milano, IT.
C-band Radar¶
We need to load and preprocess the C-band radar data first. In order to do so, we first need to load the MeteoSwiss C-band radar file.
# Set the URL and path for the cloud
URL = "https://js2.jetstream-cloud.org:8001/"
path = f"pythia/radar/erad2024"
fs = fsspec.filesystem("s3", anon=True, client_kwargs=dict(endpoint_url=URL))
fs.glob(f"{path}/*")
files = fs.glob("pythia/radar/erad2024/dda_data/*")
local_files = [
fsspec.open_local(
f"simplecache::{URL}{i}", s3={"anon": True}, filecache={"cache_storage": "."}
)
for i in files
]
radar = pyart.io.read(local_files[1])radar.info()Dealiasing¶
The next step is to dealias the radar data. We will follow the steps that were shown eariler in this course to dealias the radar data. Dealiasing is covered in detail in the Radar Cookbook.
nyquist = radar.instrument_parameters["nyquist_velocity"]["data"][0]
vel_dealias = pyart.correct.dealias_region_based(
radar,
vel_field="velocity",
nyquist_vel=nyquist,
centered=True,
)
radar.add_field("corrected_velocity", vel_dealias, replace_existing=True)Plot the data to make sure dealiasing succeeded.
display = pyart.graph.RadarDisplay(radar)
display.plot("reflectivity", vmin=-20, vmax=70, cmap="pyart_ChaseSpectral", sweep=0)display = pyart.graph.RadarDisplay(radar)
display.plot("corrected_velocity", vmin=-30, vmax=30, cmap="twilight_shifted", sweep=0)Some regions did not dealias correctly, like in the region around (-100 km, -100 km) where gate-to-gate shear exceeds the Nyquist interval in the absence of . The region-based algorithm is very sensitive to speckles that are present in the above image. Let’s despeckle the above radar data to remove some artifacts.
gatefilter = pyart.filters.GateFilter(radar)
gatefilter.exclude_above("spectrum_width", 7)
# gatefilter.exclude_below('cross_correlation_ratio', 0.8)
gatefilter = pyart.correct.despeckle_field(
radar, "reflectivity", gatefilter=gatefilter, size=36
)
display = pyart.graph.RadarDisplay(radar)
display.plot(
"reflectivity",
vmin=-20,
vmax=70,
gatefilter=gatefilter,
cmap="pyart_ChaseSpectral",
sweep=0,
)Let’s also tweak the settings of the region based dealiasing. Currently, the dealiasing will attempt to dealias regions that are separated by 100 masked gates in the radial. We will disable this feature by setting skip_between_rays to 0 so that all dealiasing across filtered gates along the azimuth is disabled.
nyquist = radar.instrument_parameters["nyquist_velocity"]["data"][0]
vel_dealias = pyart.correct.dealias_region_based(
radar,
vel_field="velocity",
nyquist_vel=nyquist,
gatefilter=gatefilter,
centered=True,
skip_between_rays=0,
)
radar.add_field("corrected_velocity", vel_dealias, replace_existing=True)Let’s plot all of the sweeps to make sure dealiasing worked correctly. A good indicator of checking this is to check if there are discontinuities in the radial velocities between sweeps on the order of the Nyquist velocity. If there are such discontinuties, further tweaking of the parameters or grounding the retrieval to a sounding may be needed. Using these parameters, the above example radar data is properly dealiased.
display = pyart.graph.RadarDisplay(radar)
fig, ax = plt.subplots(4, 5, figsize=(15, 10))
for i in range(radar.nsweeps):
display.plot(
"corrected_velocity",
vmin=-30,
vmax=30,
cmap="twilight_shifted",
sweep=i,
colorbar_label="m/s",
ax=ax[int(i / 5), i % 5],
)
ax[int(i / 5), i % 5].set_title("%2.1f deg" % radar.fixed_angle["data"][i])
ax[int(i / 5), i % 5].set_xlabel("X [km]")
ax[int(i / 5), i % 5].set_ylabel("Y [km]")
fig.tight_layout()pyart.io.write_cfradial("MonteLema.20240522.151606_dealiased.nc", radar)Gridding¶
PyDDA requires data to be gridded to Cartesian coordinates in order to retrieve the 3D wind fields. Therefore, we will use Py-ART’s grid_from_radars function in order to do the gridding. You usually want to have a grid resolution such that your features of interest are covered by four grid points. In this case, we’re at 1 km horizontal and 0.5 km vertical resolution. For more information on gridding with Py-ART, see the Radar Cookbook.
grid_limits = ((0.0, 15000.0), (-100_000.0, 100_000.0), (-100_000.0, 100_000.0))
grid_shape = (31, 201, 201)
radar.fields["corrected_velocity"]["data"] = np.ma.masked_where(
gatefilter.gate_excluded, radar.fields["corrected_velocity"]["data"]
)
radar.fields["reflectivity"]["data"] = np.ma.masked_where(
gatefilter.gate_excluded, radar.fields["reflectivity"]["data"]
)
cband_grid = pyart.map.grid_from_radars(
[radar], grid_limits=grid_limits, grid_shape=grid_shape
)
cband_ds = cband_grid.to_xarray()
cband_dsLet’s make sure the grid looks good!
cband_ds.isel(z=1).reflectivity.plot(
x="lon", y="lat", vmin=-20, vmax=70, cmap="pyart_ChaseSpectral"
)cband_ds.isel(y=35).reflectivity.plot(
x="x", y="z", vmin=-20, vmax=70, cmap="pyart_ChaseSpectral"
)X-band Data¶
Next, we need to load the X-band radar data from the archive.
radar_xband = pyart.io.read(local_files[0])Read a single file, the one closes to the UAH volume scan used before
radar.info()Visualize the data to make sure we have the correct scan.
display = pyart.graph.RadarDisplay(radar_xband)
display.plot("DBZ", vmin=0, vmax=70, cmap="pyart_ChaseSpectral", sweep=0)
plt.ylim(-100, 100)
plt.xlim(-100, 100)The X-band radar data also need to be dealiased. We will still need to filter out the noise in the above radar data.
# Use the C-band radar lat/lon as the center for the grid
gatefilter = pyart.filters.GateFilter(radar_xband)
gatefilter.exclude_above("WIDTH", 7)
gatefilter.exclude_below("RHOHV", 0.8)
gatefilter = pyart.correct.despeckle_field(
radar_xband, "DBZ", gatefilter=gatefilter, size=100
)
vel_dealias = pyart.correct.dealias_region_based(
radar_xband,
vel_field="VEL",
nyquist_vel=nyquist,
centered=True,
skip_between_rays=0,
skip_along_ray=0,
interval_splits=3,
)
radar_xband.add_field("corrected_velocity", vel_dealias, replace_existing=True)Let’s view the sweeps to make sure that dealiasing is working correctly.
display = pyart.graph.RadarDisplay(radar_xband)
fig, ax = plt.subplots(2, 4, figsize=(15, 5))
for i in range(radar_xband.nsweeps):
display.plot(
"corrected_velocity",
vmin=-30,
vmax=30,
cmap="twilight_shifted",
sweep=i,
colorbar_label="m/s",
ax=ax[int(i / 4), i % 4],
)
ax[int(i / 4), i % 4].set_title("%2.1f deg" % radar.fixed_angle["data"][i])
ax[int(i / 4), i % 4].set_xlabel("X [km]")
ax[int(i / 4), i % 4].set_ylabel("Y [km]")
fig.tight_layout()pyart.io.write_cfradial("ARPA_Lombardia.20240522.151546_dealiased.nc", radar_xband)grid_lat = radar.latitude["data"][0]
grid_lon = radar.longitude["data"][0]
xband_grid = pyart.map.grid_from_radars(
[radar_xband],
grid_limits=grid_limits,
grid_shape=grid_shape,
grid_origin=(grid_lat, grid_lon),
)
# Convert to xarray and remove the time dimension
xband_ds = xband_grid.to_xarray().squeeze()Visualize the grids¶
Let’s see what our output data looks like!
xband_ds.DBZ.isel(z=1).plot(x="lon", y="lat", cmap="Spectral_r", vmin=-20, vmax=50)Velocities look good!
xband_ds.corrected_velocity.isel(z=4).plot(
x="lon", y="lat", cmap="twilight_shifted", vmin=-30, vmax=30
)xband_ds.corrected_velocity.isel(y=55).plot(
x="x", y="z", cmap="twilight_shifted", vmin=-20, vmax=20
)PyDDA initialization¶
The 3DVAR wind retrieval first requires an initial guess at the wind field in order to start the cost function minimization process. PyDDA has support for using WRF and sounding data as an initial guess of the wind field as well as constant wind fields.
| Initalization functions in pydda.initialization module: | Functionality |
|---|---|
| make_constant_wind_field(Grid[, wind, vel_field]) | This function makes a constant wind field given a wind vector |
| make_wind_field_from_profile(Grid, profile) | This function makes a 3D wind field from a sounding. |
| make_background_from_wrf(Grid, file_path, ...) | This function makes an initalization field based off of the u and w from a WRF run. |
| make_initialization_from_era_interim(Grid[, ...]) | This function will read ERA Interim in NetCDF format and add it to the Py-ART grid specified by Grid. |
For this example, using a zero initialization field is sufficient to create an accurate representation of the wind field. Depending on your radar set up, you may need to adjust the input initial wind field to avoid artifacts at the edges of the multi-Doppler lobes.
cband_ds = pydda.io.read_from_pyart_grid(cband_grid)
xband_ds = pydda.io.read_from_pyart_grid(xband_grid)
cband_ds = pydda.initialization.make_constant_wind_field(
cband_ds, (0.0, 0.0, 0.0), "velocity"
)
xband_ds["reflectivity"] = xband_ds["DBZ"]PyDDA wind retrieval¶
The core wind retrieval function in PyDDA is done using retrieval.get_dd_wind_field. It has many potential keyword inputs that the user can enter. In this example, we are specifying:
| Input to pydda.initialization module: | Meaning | Value |
|---|---|---|
| Grids | The input grids to analyze. | [uah_grid, nexrad_grid] |
| u_init | Initial guess of u field. | u_init |
| v_init | Initial guess of u field. | v_init |
| w_init | Initial guess of u field. | w_init |
| Co | Weight for cost function related to radar observations | 1.0 |
| Cm | Weight of cost function related to mass continuity equation | 256.0 |
| Cx | Weight of cost function for smoothess in the x-direction | 1e-3 |
| Cy | Weight of cost function for smoothess in the y-direction | 1e-3 |
| Cz | Weight of cost function for smoothess in the z-direction | 1e-3 |
| Cb | Weight of cost function for sounding (background) constraint | 0 |
| frz | The freezing level in meters. This is to tell PyDDA where to use ice particle fall speeds in the wind retrieval verus liquid. | 5000. |
| filter_window | The window to apply the low pass filter on | 5 |
| mask_outside_opt | Mask all winds outside the Dual Doppler lobes | True |
| vel_name | The name of the velocity field in the radar data | ‘corrected_velocity’ |
| wind_tol | Stop optimization when the change in wind speeds between iterations is less than this value | |
| engine | PyDDA supports three backends for optimization: SciPy, JAX, and TensorFlow. | “scipy” |
grids, params = pydda.retrieval.get_dd_wind_field(
[cband_ds, xband_ds],
Co=1,
Cm=1024.0,
Cx=1e-2,
Cy=1e-2,
Cz=1e-2,
frz=5000.0,
mask_outside_opt=True,
upper_bc=1,
vel_name="corrected_velocity",
wind_tol=0.5,
engine="scipy",
)Visualize the results¶
Let’s visualize the results. There are two ways in which this data can be visualized. One way is by using PyDDA’s visualization routines. You can also use xarray to visualize the output grids.
ds = grids[0]
dsIf you are not able to get the above example to run on your laptop, pre-generated grid files are available here.
ds.reflectivity.sel(z=2000, method="nearest").plot(
cmap="pyart_ChaseSpectral", vmin=-10, vmax=60
)
ds.isel(time=0).sel(z=2000, method="nearest").w.plot.contour(
x="x", y="y", levels=np.arange(1, 5, 1)
)
plt.xlim([-50000, 25000])
plt.ylim([-100000, -25000])Saving the grids¶
We can save the grids using either xarray’s or Py-ART’s functionality. To save the output grids, use pyart.io.write_grid (PyART) functions. For example:
# Keep these lines if we can't update PyDDA in time for ERAD
grids[0].attrs["radar_name"] = "C-band"
del grids[0]["time"].attrs["units"]
grids[1].attrs["radar_name"] = "X-band"
del grids[1]["time"].attrs["units"]
grids[0].to_netcdf("output_grid_Cband.nc")
grids[1].to_netcdf("output_grid_Xband.nc")In order to load the grids again, we can just use PyART’s pyart.io.read_grid procedure.
# For live demos, load these since retrieval takes about 7 minutes
grids = [
pydda.io.read_grid("output_grid_Cband.nc"),
pydda.io.read_grid("output_grid_Xband.nc"),
]PyDDA visualization routines¶
PyDDA’s visualization routines support the native PyART grids that are output by PyDDA. These routines have an advantage over xarray’s plotting routines for adjusting your barb and quiver size by specifying their using parameters that are in scales of kilometers. This makes it easier to plot barb and quiver plots compared to using xarray’s functionality.
For example, the documentation for pydda.vis.plot_horiz_xsection_quiver is given below.
?pydda.vis.plot_horiz_xsection_quiverPyDDA has the following visualization routines for your sets of grids:
| Procedure | Description |
|---|---|
| plot_horiz_xsection_barbs(Grids[, ax, ...]) | Horizontal cross section of winds from wind fields generated by PyDDA using barbs. |
| plot_xz_xsection_barbs(Grids[, ax, ...]) | Cross section of winds from wind fields generated by PyDDA in the X-Z plane using barbs. |
| plot_yz_xsection_barbs(Grids[, ax, ...]) | Cross section of winds from wind fields generated by PyDDA in the Y-Z plane using barbs. |
| plot_horiz_xsection_barbs_map(Grids[, ax, ...]) | Horizontal cross section of winds from wind fields generated by PyDDA onto a geographical map using barbs. |
| plot_horiz_xsection_streamlines(Grids[, ax, ...]) | Horizontal cross section of winds from wind fields generated by PyDDA using streamlines. |
| plot_xz_xsection_streamlines(Grids[, ax, ...]) | Cross section of winds from wind fields generated by PyDDA in the X-Z plane using streamlines. |
| plot_yz_xsection_streamlines(Grids[, ax, ...]) | Cross section of winds from wind fields generated by PyDDA in the Y-Z plane using streamlines. |
| plot_horiz_xsection_streamlines_map(Grids[, ...]) | Horizontal cross section of winds from wind fields generated by PyDDA using streamlines. |
| plot_horiz_xsection_quiver(Grids[, ax, ...]) | Horizontal cross section of winds from wind fields generated by PyDDA using quivers. |
| plot_xz_xsection_quiver(Grids[, ax, ...]) | Cross section of winds from wind fields generated by PyDDA in the X-Z plane using quivers. |
| plot_yz_xsection_quiver(Grids[, ax, ...]) | Cross section of winds from wind fields generated by PyDDA in the Y-Z plane using quivers. |
| plot_horiz_xsection_quiver_map(Grids[, ax, ...]) | Horizontal cross section of winds from wind fields generated by PyDDA using quivers onto a geographical map. |
Let’s show a quiver plot of this storm!
We have specified the quivers to be 4 km apart and moved the key to the bottom right with the specific length indicating 20 m/s winds. Let’s look at the 3 km level.
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
pydda.vis.plot_horiz_xsection_quiver(
grids,
quiver_spacing_x_km=2.5,
quiver_spacing_y_km=2.5,
quiver_width=0.005,
vmax=60,
quiverkey_len=10.0,
w_vel_contours=np.arange(1, 5, 1),
level=4,
cmap="pyart_ChaseSpectral",
ax=ax,
quiverkey_loc="bottom_right",
)
plt.xlim([-75, 25])
plt.ylim([-100, -25])We can zoom in and modify the plot using standard matplotlib functions on the axis handle.
It is much easier to see updrafts being placed just to the outside of the strongest precipitation, with potential new growth in the north of the domain with updraft velocities > 7 m/s. The precipitation is downwind of the updraft as we would expect.
Updrafts are right tilted due to the horizontal wind shear. The horizontal wind shear also causes the most intense precipitation to be downshear of the updraft. This therefore shows us that we have a good quality wind retrieval below about 5 km in altitude.
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
pydda.vis.plot_xz_xsection_quiver(
grids,
quiver_spacing_x_km=2.0,
quiver_spacing_z_km=1.0,
quiver_width=0.005,
quiverkey_len=10.0,
w_vel_contours=np.arange(1, 5, 1),
level=37,
cmap="pyart_ChaseSpectral",
ax=ax,
quiverkey_loc="top_right",
)
ax.set_xlim([-75, -25])
ax.set_ylim([0, 13])Let’s view a horizontal cross section with barbs!
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
pydda.vis.plot_horiz_xsection_barbs(
grids,
barb_spacing_x_km=6.0,
barb_spacing_y_km=6.0,
w_vel_contours=np.arange(1, 10, 1),
level=6,
cmap="pyart_ChaseSpectral",
ax=ax,
)
# ax.set_xlim([-60, 0])
# ax.set_ylim([0, 70])cband_ds["time"][0].dt.year.values