ARM Logo

Quantitative Precipitation Estimate (QPE) Analysis with Py-ART¶

Overview¶

Within this notebook, we will cover:

Calculation of QPE from various radar fields
Genearting a gridded QPE product
Comparison against operational models (MRMS)

Prerequisites¶

Concepts	Importance	Notes
Py-ART Basics	Helpful	Basic features
Weather Radar Basics	Helpful	Background Information
Matplotlib Basics	Helpful	Basic plotting
NumPy Basics	Helpful	Basic arrays
Xarray Basics	Helpful	Multi-dimensional arrays

Time to learn: 15 minutes

Imports¶

import os
import warnings
import glob

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import pyart
from pyart.testing import get_test_data
import xradar as xd

warnings.filterwarnings('ignore')


## 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

Investigate the Available Data¶

file_list = sorted(glob.glob("/data/project/ARM_Summer_School_2025/radar/csapr2/ppi/*20250519*"))
file_list

[]

dt = xd.io.open_cfradial1_datatree(file_list[0])
dt

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[3], line 1
----> 1 dt = xd.io.open_cfradial1_datatree(file_list[0])
      2 dt

IndexError: list index out of range

radar = pyart.xradar.Xradar(dt)
display = pyart.graph.RadarDisplay(radar)

radar

# Generate matplotlib figure and axe array objects
fig, axarr = plt.subplots(2, 2, figsize=[20, 12])
plt.subplots_adjust(hspace=0.35)

# reflectivity
display.plot('reflectivity', sweep=0, ax=axarr[0, 0], cmap='ChaseSpectral')

# differential reflectivity
display.plot("differential_reflectivity", sweep=0, ax=axarr[0, 1], cmap="RdBu_r", vmin=-2, vmax=2)

# doppler velocity
display.plot("mean_doppler_velocity", sweep=0, ax=axarr[1, 0], cmap="RdBu_r", vmin=-16, vmax=16)

# differential phase
display.plot("specific_differential_phase", sweep=0, ax=axarr[1, 1], cmap="ChaseSpectral")

QPE Functions¶

def reflectivity_rain(radar, refl="reflectivity", alpha=0.0376, beta=0.6112):
    """
    Function to calculate rainfall rates from radar reflectivity factor

    Inputs
    ------
    radar : Py-ART Radar Object
        Py-ART radar object to extract reflectivity field from
    refl : str
        Specific name of reflectivity field within radar object
    alpha : float
        fit parameter
    beta : float
        fit parameter

    Outputs
    -------
    radar : Py-ART Radar Object
        Py-ART radar object with rainfall estimate from reflectivity included
    """
    # define a gatefilter to apply the relationship to
    gatefilter_z = pyart.correct.GateFilter(radar)
    gatefilter_z.exclude_above(refl, 35)
    # Apply the gatefilter to the rain rate
    masked_z = np.ma.masked_array(radar.fields[refl]['data'], mask=gatefilter_z.gate_excluded) 
    # Apply the R(Z) relationship
    rr_data = alpha * np.ma.power(np.ma.power(10.0, 0.1 * masked_z), beta)
    # define the dictionary structure for the rain rate data
    rain = pyart.config.get_metadata("radar_estimated_rain_rate")
    rain["long_name"] = "R(Z) Radar Estimated Rain Rate"
    rain["standard_name"] = "R(Z) Radar Estimated Rain Rate"
    rain["data"] = rr_data
    # add the field back into the radar object
    radar.add_field("rain_z", rain)

    return radar

def kdp_rain(radar, phase="specific_differential_phase", alpha=25.1, beta=0.777):
    """
    Function to calculate rainfall rates from specific differential phase
    Inputs
    ------
    radar : Py-ART Radar Object
        Py-ART radar object to extract reflectivity field from
    refl : str
        Specific name of reflectivity field within radar object
    alpha : float
        fit parameter
    beta : float
        fit parameter

    Outputs
    -------
    radar : Py-ART Radar Object
        Py-ART radar object with rainfall estimate from reflectivity included
    """
    # define a gatefilter to apply the relationship to
    gatefilter_kdp = pyart.correct.GateFilter(radar)
    gatefilter_kdp.exclude_below('reflectivity', 35)
    # Apply the gatefilter to the rain rate
    masked_z = np.ma.masked_array(radar.fields[phase]['data'], mask=gatefilter_kdp.gate_excluded) 
    # define the reflectivity data
    reflect = radar.fields[phase]["data"]
    rr_data = alpha * np.ma.power(np.ma.power(10.0, 0.1 * masked_z), beta)
    # define the dictionary structure for the rain rate data
    rain = pyart.config.get_metadata("radar_estimated_rain_rate")
    rain["long_name"] = "R(KDP) Radar Estimated Rain Rate"
    rain["standard_name"] = "R(KDP) Radar Estimated Rain Rate"
    rain["data"] = rr_data
    # add the field back into the radar object
    radar.add_field("rain_kdp", rain)

    return radar

# Apply the Radar estimated rain rates 
radar = reflectivity_rain(radar)
radar = kdp_rain(radar)
radar["sweep_1"]

display = pyart.graph.RadarDisplay(radar)
# Generate matplotlib figure and axe array objects
fig, axarr = plt.subplots(1, 2, figsize=[14, 5])
plt.subplots_adjust(wspace=0.2, hspace=0.35)

# reflectivity
display.plot('rain_z', sweep=0, ax=axarr[0], cmap='ChaseSpectral')

# differential reflectivity
display.plot("rain_kdp", sweep=0, ax=axarr[1], cmap="RdBu_r")

Combined Radar Estimated Rainfall Product¶

radar.fields["rain_z"]["data"].mask

combined_data = np.where(~radar.fields["rain_z"]["data"].mask, radar.fields["rain_z"]["data"], radar.fields["rain_kdp"]["data"])
# Combine the masks using logical OR (mask where either is masked)
combined_mask = np.ma.mask_or(radar.fields["rain_z"]["data"].mask, radar.fields["rain_kdp"]["data"].mask)

# Merge data and apply the combined mask
merged = np.ma.array(combined_data, mask=combined_mask)

combined_data

# define the dictionary structure for the rain rate data
rain = pyart.config.get_metadata("radar_estimated_rain_rate")
rain["long_name"] = "R(Z+KDP) Radar Estimated Rain Rate"
rain["standard_name"] = "R(Z+KDP) Radar Estimated Rain Rate"
rain["data"] = combined_data
# add the field back into the radar object
radar.add_field("rain_combined", rain)

display = pyart.graph.RadarDisplay(radar)
# Generate matplotlib figure and axe array objects
fig, axarr = plt.subplots(1, 2, figsize=[14, 5])
plt.subplots_adjust(wspace=0.2, hspace=0.35)

# reflectivity
display.plot('rain_combined', sweep=0, ax=axarr[0], cmap='ChaseSpectral', vmax=10)
# reflectivity
display.plot('reflectivity', sweep=0, ax=axarr[1], cmap='ChaseSpectral')

Create a Gridded QPE Product for Future Comparison with Model/MRMS¶

def bnf_grid(radar, 
             z_limits=(250., 15_000.), 
             y_limits=(-30_000., 30_000), 
             x_limits=(-30_000., 30_000),
             resolution=250
):
    """
    Function to create a Py-ART grid object from a given radar file

    Inputs
    ------
    radar : Py-ART radar object
        Py-ART radar object to create a grid object from
    z_limits : tuple
        vertical dimension grid limits
    y_limits : tuple
        longitude dimension grid limits
    x_limits : tuple
        latitude dimension grid limits
    resolution : float
        desired resolution of our grid object

    Calls
    -----
    compute_number_of_points
        compute number of gates in each direction

    Outputs
    -------
    grid : Py-ART Grid object
    """
    def compute_number_of_points(extent, resolution):
        return int((extent[1] - extent[0])/resolution)

    z_grid_points = compute_number_of_points(z_limits, resolution)
    x_grid_points = compute_number_of_points(x_limits, resolution)
    y_grid_points = compute_number_of_points(y_limits, resolution)

    grid = pyart.map.grid_from_radars([radar],
                                      grid_shape=(z_grid_points,
                                                  y_grid_points,
                                                  x_grid_points),
                                      grid_limits=(z_grid_limits,
                                                   y_grid_limits,
                                                   x_grid_limits),
    )

    return grid

grid = bnf_grid(radar)

display = pyart.graph.GridMapDisplay(grid)
display.plot_grid('rain_combined',
                  level=1,
                  vmin=0.1,
                  vmax=15,
                  cmap='HomeyerRainbow')

Radar Data with Py-ART

Py-ART Corrections

Timeseries Data with ACT

ACT Basics with BNF