Hydrometeor Classification with Custom Frequency Settings#

This script shows how to use hydrometeor classification for X-band radar data. We are reading radar data, plotting some variables of interest and applying the classification to identify types of precipitation.

.. note:: The script initially attempts hydrometeor classification without specific radar frequency information for band selection.

import matplotlib.pyplot as plt
import numpy as np
from open_radar_data import DATASETS

import pyart

filename = DATASETS.fetch("gucxprecipradarcmacppiS2.c1.20220314.025840.nc")
radar = pyart.io.read_cfradial(filename)

figure = plt.figure(figsize=(15, 4))

ax1 = plt.subplot(1, 3, 1)
display = pyart.graph.RadarDisplay(radar)
ax1 = display.plot("DBZ", vmin=0, vmax=50)  # DBZ corrected_reflectivity
plt.xlim(-20, 20)
plt.ylim(-20, 20)

ax2 = plt.subplot(1, 3, 2)
ax2 = display.plot("corrected_differential_reflectivity", cmap="Carbone42")  # ZDR
plt.xlim(-20, 20)
plt.ylim(-20, 20)

ax3 = plt.subplot(1, 3, 3)
ax3 = display.plot("corrected_specific_diff_phase", cmap="Carbone42")  # KDP
plt.xlim(-20, 20)

# ### When instrument parameters does not have radar frequency info.

print(radar.instrument_parameters)


# This shows an issue where radar frequency information is missing. Without this hydrometeor classification will default to C-band.

# Get classification
hydromet_class = pyart.retrieve.hydroclass_semisupervised(
    radar,
    refl_field="corrected_reflectivity",
    zdr_field="corrected_differential_reflectivity",
    kdp_field="filtered_corrected_specific_diff_phase",
    rhv_field="RHOHV",
    temp_field="sounding_temperature",
)["hydro"]

radar.add_field("hydro_classification", hydromet_class, replace_existing=True)


# Use `radar_freq` parameters
# To address this issue, radar frequency information can be supplied to the function with `radar_freq` parameter.


# Get classification
hydromet_class = pyart.retrieve.hydroclass_semisupervised(
    radar,
    refl_field="corrected_reflectivity",
    zdr_field="corrected_differential_reflectivity",
    kdp_field="filtered_corrected_specific_diff_phase",
    rhv_field="RHOHV",
    temp_field="sounding_temperature",
    radar_freq=9.2e9,
)["hydro"]

radar.add_field("hydro_classification", hydromet_class, replace_existing=True)


# Add radar frequency to the radar object
# Incorporating radar frequency into the radar object enhances processing pipeline.

# Add X-band frequency information to radar.instrument_parameters
radar.instrument_parameters["frequency"] = {
    "long_name": "Radar frequency",
    "units": "Hz",
    "data": [9.2e9],
}

radar.instrument_parameters


# Let's run the classification again and the warning should change telling the radar frequency from instrument parameters is used.


hydromet_class = pyart.retrieve.hydroclass_semisupervised(
    radar,
    refl_field="corrected_reflectivity",
    zdr_field="corrected_differential_reflectivity",
    kdp_field="filtered_corrected_specific_diff_phase",
    rhv_field="RHOHV",
    temp_field="sounding_temperature",
    radar_freq=9.2e9,
)["hydro"]

radar.add_field("hydro_classification", hydromet_class, replace_existing=True)


# Note that the frequency used here is from the radar object, not the user supplied.


# plotting

import matplotlib.colors as colors

hid_colors = [
    "White",
    "LightBlue",
    "MediumBlue",
    "DarkOrange",
    "LightPink",
    "Cyan",
    "DarkGray",
    "Lime",
    "Yellow",
    "Red",
]
cmaphid = colors.ListedColormap(hid_colors)
cmapmeth = colors.ListedColormap(hid_colors[0:6])
cmapmeth_trop = colors.ListedColormap(hid_colors[0:7])


def adjust_fhc_colorbar_for_pyart(cb):
    cb.set_ticks(np.arange(1.4, 9, 0.9))
    cb.ax.set_yticklabels(
        [
            "Aggregates",
            "Ice Crystals",
            "Light rain",
            "Rain",
            "Vertical Ice",
            "Wet snow",
            "HD Graupel",
            "Melting hail",
            "Dry hail or high-density graupel",
        ]
    )
    cb.ax.set_ylabel("")
    cb.ax.tick_params(length=0)
    return cb


def adjust_meth_colorbar_for_pyart(cb, tropical=False):
    if not tropical:
        cb.set_ticks(np.arange(1.25, 5, 0.833))
        cb.ax.set_yticklabels(
            ["R(Kdp, Zdr)", "R(Kdp)", "R(Z, Zdr)", "R(Z)", "R(Zrain)"]
        )
    else:
        cb.set_ticks(np.arange(1.3, 6, 0.85))
        cb.ax.set_yticklabels(
            ["R(Kdp, Zdr)", "R(Kdp)", "R(Z, Zdr)", "R(Z_all)", "R(Z_c)", "R(Z_s)"]
        )
    cb.ax.set_ylabel("")
    cb.ax.tick_params(length=0)
    return cb


def two_panel_plot(
    radar,
    sweep=0,
    var1="corrected_reflectivity",
    vmin1=0,
    vmax1=65,
    cmap1="RdYlBu_r",
    units1="dBZ",
    var2="corrected_differential_reflectivity",
    vmin2=-5,
    vmax2=5,
    cmap2="RdYlBu_r",
    units2="dB",
    return_flag=False,
    xlim=[-150, 150],
    ylim=[-150, 150],
):
    display = pyart.graph.RadarDisplay(radar)
    fig = plt.figure(figsize=(13, 5))
    ax1 = fig.add_subplot(121)
    display.plot_ppi(
        var1,
        sweep=sweep,
        vmin=vmin1,
        vmax=vmax1,
        cmap=cmap1,
        colorbar_label=units1,
        mask_outside=True,
    )
    display.set_limits(xlim=xlim, ylim=ylim)
    ax2 = fig.add_subplot(122)
    display.plot_ppi(
        var2,
        sweep=sweep,
        vmin=vmin2,
        vmax=vmax2,
        cmap=cmap2,
        colorbar_label=units2,
        mask_outside=True,
    )
    display.set_limits(xlim=xlim, ylim=ylim)
    if return_flag:
        return fig, ax1, ax2, display


lim = [-20, 20]
fig, ax1, ax2, display = two_panel_plot(
    radar,
    sweep=0,
    var1="corrected_reflectivity",
    var2="hydro_classification",
    vmin2=0,
    vmax2=10,
    cmap2=cmaphid,
    units2="",
    return_flag=True,
    xlim=lim,
    ylim=lim,
)
display.cbs[1] = adjust_fhc_colorbar_for_pyart(display.cbs[1])