QVP plot

QPE & QVPs¶

Overview¶

In this notebook, we demonstrate how utilizing radar data in the Analysis-Ready Cloud-Optimized (ARCO) format enables efficient computation of Quantitative Precipitation Estimates (QPE) and Quasi-Vertical Profiles (QVP). The ARCO format ensures that radar data is pre-processed, clean, and well-organized, significantly reducing the time spent on data preparation and cleaning. By leveraging ARCO radar data, we can focus more on scientific analysis.

QPE Demo
QVP Demo

Prerequisites¶

Concepts	Importance	Notes
Intro to Xarray	Necessary	Basic features
Radar Cookbook	Necessary	Radar basics
Intro to Zarr	Necessary	Zarr basics
Intro to Hvplot	Necessary	Interactive visualization basics

Time to learn: 45 minutes

Imports¶

import s3fs
import xarray as xr
import fsspec
import numpy as np
import hvplot.xarray
import matplotlib.pyplot as plt
import holoviews as hv

ARCO radar dataset¶

We build a dataset (object storage) containing both X and C-band radar data. This dataset is located at the Pythia S3 Bucket.

## S3 bucket connection
URL = 'https://js2.jetstream-cloud.org:8001/'
path = f'pythia/radar/erad2024'
fs = s3fs.S3FileSystem(anon=True, client_kwargs=dict(endpoint_url=URL))
file =  s3fs.S3Map(f"{path}/zarr_radar/erad_2024.zarr", s3=fs)

Let’s read our ARCO dataset using xr.backends.api.open_datatree

%%time
dtree = xr.backends.api.open_datatree(
    file, 
    engine='zarr', 
    consolidated=True, 
    chunks={}
)

CPU times: user 2.99 s, sys: 303 ms, total: 3.29 s
Wall time: 9.58 s

dtree

list(dtree.children)

['cband', 'xband']

C-band radar data¶

We can access a nested dataset by using the path as follows:

dtree["cband"]

# all sweeps/nodes within the cband datatree
list(dtree["cband"].children)

['sweep_0',
 'sweep_1',
 'sweep_10',
 'sweep_11',
 'sweep_12',
 'sweep_13',
 'sweep_14',
 'sweep_15',
 'sweep_16',
 'sweep_17',
 'sweep_18',
 'sweep_19',
 'sweep_2',
 'sweep_3',
 'sweep_4',
 'sweep_5',
 'sweep_6',
 'sweep_7',
 'sweep_8',
 'sweep_9']

Accessing one of the sweeps (“sweep_2”):

dsc = dtree["cband/sweep_2"].ds

display(dsc)

Let’s use hvplot.quadmesh to create interactive plots

# set constant colorbar limits
ref_C = dsc.reflectivity.compute().hvplot.quadmesh(
    x="x",
    y="y",
    groupby="volume_time",
    clim=(-10, 60),
    cmap="ChaseSpectral",
    width=600,
    height=500,
    widget_type="scrubber",
    widget_location="bottom",
    rasterize=True,
)

ref_C

X-band radar data¶

Similarly, we can access X-band radar data as follows:

list(dtree["xband"].children)

['sweep_0',
 'sweep_1',
 'sweep_2',
 'sweep_3',
 'sweep_4',
 'sweep_5',
 'sweep_6',
 'sweep_7']

dsx = dtree["xband/sweep_0"].ds

display(dsx)

Interactive plot using hvplot.quadmesh

# set constant colorbar limits
ref_X = dsx.DBZ.compute().hvplot.quadmesh(
    x="x",
    y="y",
    groupby="volume_time",
    clim=(-10, 60),
    cmap="ChaseSpectral",
    width=600,
    height=500,
    widget_type="scrubber",
    widget_location="bottom",
    rasterize=True,
)
ref_X

Quantitave Precipitation Estimation (QPE)¶

QPE is a critical process in meteorology, providing measurements of rainfall intensity and accumulation from radar data. One of the foundational methods for QPE is based on the Z-R relationship, which establishes a statistical relationship between radar reflectivity (Z) and rainfall rate (R). This empirical relationship, derived from observations, is commonly expressed as

𝑍=\alpha R^{\beta}

(1)

where α and β are constants.

def rain_depth(z: xr.DataArray, a: float=200.0, b: float=1.6, t:int=5) -> xr.DataArray:
    """
    Estimates rainfall depth using radar reflectivity and a Z-R relationship.
    
    This function computes Quantitative Precipitation Estimation (QPE) by converting
    radar reflectivity (Z) into rainfall rate (R) using the Z-R relationship and
    then integrating over time to estimate the total rainfall depth.
    
    Parameters:
    -----------
    z : xr.DataArray
        Radar reflectivity in dBZ. This should be a multi-dimensional Xarray DataArray.
    a : float, optional
        The alpha (a) parameter in the Z-R relationship. Default is 200.0, corresponding
        to the Marshall and Palmer (1948) relationship.
    b : float, optional
        The beta (b) parameter in the Z-R relationship. Default is 1.6, also from the
        Marshall and Palmer (1948) relationship.
    t : int, optional
        Time integration period in minutes, used to convert rainfall rates into
        accumulated depth. Default is 5 minutes.
    
    Returns:
    --------
    xr.DataArray
        A DataArray representing the estimated rainfall depth in the same dimensions
        as the input radar reflectivity. The units of the returned DataArray will be
        consistent with the time integration provided (e.g., mm for 5-minute accumulation).
    
    Notes:
    ------
    - The Z-R relationship used is of the form Z = a * R^b, where Z is in linear units.
    - The function first converts the radar reflectivity from dBZ to linear units (Z),
      then computes the rainfall rate (R), and finally multiplies by the time integration
      period to obtain the rainfall depth.
    
    Example:
    --------
    To compute the rainfall depth over a 5-minute period using reflectivity data:
    
    >>> rainfall_depth = rain_depth(z, a=200.0, b=1.6, t=5)
    
    This will return the estimated rainfall depth in millimeters, assuming the default
    parameters for the Marshall and Palmer (1955) Z-R relationship.
    
    """
    # Convert reflectivity from dBZ to linear units
    z_lin = 10 ** (z / 10) 
    # Compute rainfall depth using the Z-R relationship and time integration
    
    return ((1 / a) ** (1 / b) * z_lin ** (1 / b)) * (t / 60) # rainfall depth

Let’s apply this fucntion to our radar dataset

r_depth = rain_depth(dsc.reflectivity)

r_depth

As the coordinates depend on the volume_time dimension, they will disappear when aggregating along this dimension. Therefore, we need to store these coordinates so that we can add them back to our xr.DataArray

x = dsc.x.isel(volume_time=0)
y = dsc.y.isel(volume_time=0)
z = dsc.z.isel(volume_time=0)

Let’s compute the total rainfall depth across the entire dataset.

r_total = r_depth.sum("volume_time")

Adding coordinates back

r_total = r_total.assign_coords({"x":x, "y":y, "z":z})

Plotting the total rainfall depth

%%time
fig, ax = plt.subplots(figsize=(6, 5))
c = r_total.plot.pcolormesh(
    x='x', 
    y='y', 
    cmap='ChaseSpectral', 
    vmin=0,
    vmax=50,
    add_colorbar=False
)
plt.colorbar(c, ax=ax, label="Rainfall depth [mm]")
ax.set_title("")

CPU times: user 842 ms, sys: 24.4 ms, total: 867 ms
Wall time: 1.23 s

Now it’s your turn. Try computing the QPE for the X-band radar dataset

# Computing rainfall depth

# saving coordinates

# Acummulating rainfall depths

# adding back x,y, and z coordinates

# Plotting rainfall accumulation

Quasi-Vertical Profile (QVP)¶

Quasi-Vertical Profiles (QVP) are a radar analysis technique that provides vertical profiles of atmospheric phenomena by averaging radar data over a specific azimuthal range. This method simplifies the study of storm structures, revealing vertical distributions of reflectivity, velocity, and other key parameters. QVPs are valuable for understanding storm dynamics and the development of severe weather.

The following function will help us computing QVPs

def compute_qvp(ds: xr.Dataset, var="reflectivity")-> xr.DataArray:
    """
    Computes a Quasi-Vertical Profile (QVP) from a radar time-series dataset.
    
    This function averages the specified variable over the azimuthal dimension
    to produce a QVP. If the variable is in dBZ (a logarithmic scale), it converts
    the values to linear units before averaging and then converts the result 
    back to dBZ.
    
    Parameters:
    -----------
    ds : xr.Dataset
        The Xarray Dataset containing the radar data. This dataset should include
        multiple sweeps, azimuth angles, and range gates.
    var : str, optional
        The variable to be averaged to create the QVP. 
        Default is "reflectivity".
    
    Returns:
    --------
    xr.DataArray
        A DataArray representing the QVP for the specified variable. The result
        is averaged over azimuth and adjusted for height using the mean sweep 
        elevation angle.
    
    Notes:
    ------
    - If the variable is in dBZ units, the function converts it to linear units 
      before averaging to ensure accurate results, then converts it back to dBZ.
    - The QVP is calculated by adjusting the range gates to height using the 
      sine of the mean sweep elevation angle.
    
    Example:
    --------
    To compute a QVP for reflectivity:
    
    >>> qvp_reflectivity = compute_qvp(ds, var="reflectivity")
    
    The resulting QVP will be in dBZ and aligned along the height dimension.
    
    """
    
    units: str = ds[var].attrs['units']
    if units.startswith("dB"):
        qvp = 10 ** (ds[var] / 10)
        qvp = qvp.mean("azimuth")
        qvp = 10 * np.log10(qvp)
    else:
        qvp = ds[var]
        qvp = qvp.mean("azimuth")
        
    qvp = qvp.assign_coords({"range":(qvp.range.values * 
                                 np.sin(ds.sweep_fixed_angle.values.mean() *  
                                        np.pi / 180.))})
    qvp = qvp.rename(f"qvp_{var}")
    qvp = qvp.rename({"range": "height"})
    return qvp

Let’s compute the QVP for the X-band radar dataset at the highest elevation angle (“sweep_7”).

ds_x = dtree["xband/sweep_7"].ds

display(ds_x)

We can use our compute_qvp function to create the reflectivity QVP.

qvp_ref = compute_qvp(ds_x, var="DBZ")

display(qvp_ref)

We can visualize the QVP results using the .plot functionality within xarray library.

qvp_ref.sel(height=slice(0, 1.2e4)).plot(
    x="volume_time",
    y="height",
    cmap="ChaseSpectral",
    vmin=-10,
    vmax=50
)

Let’s try to create a figure similar to the one in Ryzhkov et al. (2016) by estimating the QVP for differential reflectivity (ZDR), the cross-correlation coefficient (RHOHV), and the differential phase (PHIDP).

qvp_zdr = compute_qvp(ds_x, var="ZDR")
qvp_rhohv = compute_qvp(ds_x, var="RHOHV")
qvp_phidp = compute_qvp(ds_x, var="PHIDP")

Let’s create the figure

%%time
fig, axs = plt.subplots(2, 2, figsize=(12, 5), sharey=True, sharex=True)

cf = qvp_ref.sel(height=slice(0, 1.2e4)).plot.contourf(
    x="volume_time",
    y="height",
    cmap="ChaseSpectral",
    vmin=-10,
    vmax=50, 
    ax=axs[0][0],
    levels=np.linspace(-10, 50, 61),
    add_colorbar=False,
)
axs[0][0].set_title(r"$Z$")
axs[0][0].set_xlabel("")
axs[0][0].set_ylabel(r"$Height \ [m]$")

cbar = plt.colorbar(cf, ax=axs[0][0], 
                    label=r"$Reflectivity \ [dBZ]$", 
                   )


cf1 = qvp_zdr.sel(height=slice(0, 1.2e4)).plot.contourf(
    x="volume_time",
    y="height",
    cmap="ChaseSpectral",
    vmin=-1,
    vmax=5, 
    ax=axs[0][1],
    levels=np.linspace(-1, 5, 11),
    add_colorbar=False,
)
axs[0][1].set_title(r"$Z_{DR}$")
axs[0][1].set_xlabel("")
axs[0][1].set_ylabel(r"")

cbar = plt.colorbar(cf1, ax=axs[0][1], 
                    label=r"$Diff. \ Reflectivity \ [dB]$", 
                   )

cf2 = qvp_rhohv.sel(height=slice(0, 1.2e4)).plot.contourf(
    x="volume_time",
    y="height",
    cmap="ChaseSpectral",
    vmin=0,
    vmax=1, 
    ax=axs[1][0],
    levels=np.linspace(0, 1, 101),
    add_colorbar=False,
)
axs[1][0].set_title(r"$\rho _{HV}$")
axs[1][0].set_ylabel(r"$Height \ [m]$")
axs[1][0].set_xlabel(r"$Time \ [UTC]$")

cbar = plt.colorbar(cf2, ax=axs[1][0], 
                    label=r"$Cross-Correlation \ Coef.$", 
                   )
cf3 = qvp_phidp.sel(height=slice(0, 1.2e4)).plot.contourf(
    x="volume_time",
    y="height",
    cmap="jet",
    vmin=0,
    vmax=100, 
    ax=axs[1][1],
    levels=np.linspace(0, 100, 101),
    add_colorbar=False,

)
axs[1][1].set_title(r"$\theta _{DP}$")
axs[1][1].set_xlabel(r"$Time \ [UTC]$")
axs[1][1].set_ylabel(r"")

cbar = plt.colorbar(cf3, ax=axs[1][1], 
                    label=r"$Differential \ Phase \ [deg]$", 
                   )

fig.tight_layout()
plt.savefig('../images/QVP.svg',  bbox_inches='tight')

CPU times: user 2.57 s, sys: 136 ms, total: 2.71 s
Wall time: 5.4 s

Now it’s your turn. Try computing the QPE for the C-band radar dataset using the 20 deg elevation angle (“sweep_15”)

## Select the sweep_15 from the dtree

## compute the QVP for "reflectivity", "differential_reflectivity", 
## "uncorrected_cross_correlation_ratio", and "uncorrected_differential_phase"

# create the figure

Summary¶

In this notebook, we successfully computed Quantitative Precipitation Estimation (QPE) and Quasi-Vertical Profiles (QVP) for both X-band and C-band radar data using the Analysis-Ready Cloud-Optimized (ARCO) dataset. By leveraging the ARCO format, we were able to streamline the data processing, allowing us to efficiently apply our custom functions for QPE and QVP computation. This approach demonstrated the effectiveness of ARCO datasets in facilitating advanced radar data analysis with minimal preprocessing effort.

Resources and references¶

Ryzhkov, A., P. Zhang, H. Reeves, M. Kumjian, T. Tschallener, S. Trömel, and C. Simmer, 2016: Quasi-Vertical Profiles—A New Way to Look at Polarimetric Radar Data. J. Atmos. Oceanic Technol., 33, 551–562, Ryzhkov et al. (2016)
Marshall, J. S.; Palmer, W. M. (1948). “The distribution of raindrops with size”. Journal of Meteorology. 5 (4): 165–166. Marshall & Palmer (1948)
Xradar
Radar cookbook
Py-Art landing page
Wradlib landing page

References¶

Ryzhkov, A., Zhang, P., Reeves, H., Kumjian, M., Tschallener, T., Trömel, S., & Simmer, C. (2016). Quasi-Vertical Profiles—A New Way to Look at Polarimetric Radar Data. Journal of Atmospheric and Oceanic Technology, 33(3), 551–562. 10.1175/jtech-d-15-0020.1
Marshall, J. S., & Palmer, W. M. K. (1948). THE DISTRIBUTION OF RAINDROPS WITH SIZE. Journal of Meteorology, 5(4), 165–166. https://doi.org/10.1175/1520-0469(1948)005<;0165:tdorws>2.0.co;2

ARCO Radar Datasets

Analysis-Ready Cloud-Optimized Datasets

Spaceborne - Ground Radar Matching and Calibration

Spaceborne-Ground Radar Calibration (Application)