diff --git a/examples/retrieve/hydrometeor-class.ipynb b/examples/retrieve/hydrometeor-class.ipynb deleted file mode 100644 index 5421256d03..0000000000 --- a/examples/retrieve/hydrometeor-class.ipynb +++ /dev/null @@ -1,167 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "be2d2ecb-8e76-432c-b7e1-3b92eb61e425", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import cartopy.crs as ccrs\n", - "import cartopy\n", - "import matplotlib.pyplot as plt\n", - "\n", - "import pyart" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "40d2e0f8-4695-4f3c-9bad-c866536f260c", - "metadata": {}, - "outputs": [], - "source": [ - "# Read in some test data\n", - "filename = pyart.testing.get_test_data(\"swx_20120520_0641.nc\")\n", - "radar = pyart.io.read(filename)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3ecbd17b-8f94-43b7-84dd-a69f7e071646", - "metadata": {}, - "outputs": [], - "source": [ - "list(radar.fields)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "584a1954-c1f0-4690-8423-c0053d54b712", - "metadata": {}, - "outputs": [], - "source": [ - "radar.fields[\"diff_phase\"]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6cb555a2-549e-4e10-adb8-938be66aeb4e", - "metadata": {}, - "outputs": [], - "source": [ - "gatefilter = pyart.filters.moment_and_texture_based_gate_filter(\n", - " radar,\n", - " zdr_field=\"diff_reflectivity\",\n", - " rhv_field=\"copol_coeff\",\n", - " phi_field=\"dp_phase_shift\",\n", - " refl_field=\"corrected_reflectivity_horizontal\",\n", - ")\n", - "# gatefilter.exclude_below('signal_to_noise_ratio', 10)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3178c6c3-4b84-4712-af6c-27d04297599c", - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(1, 2, figsize=(12, 6), sharex=True, sharey=True)\n", - "display = pyart.graph.RadarDisplay(radar)\n", - "display.plot_ppi(\n", - " \"corrected_reflectivity_horizontal\",\n", - " 0,\n", - " vmin=0,\n", - " vmax=60.0,\n", - " ax=ax[0],\n", - " colorbar_label=\"Raw Ref\",\n", - " cmap=\"pyart_HomeyerRainbow\",\n", - ")\n", - "display.plot_ppi(\n", - " \"corrected_reflectivity_horizontal\",\n", - " 0,\n", - " vmin=0,\n", - " vmax=60.0,\n", - " gatefilter=gtfilter,\n", - " cmap=\"pyart_HomeyerRainbow\",\n", - " ax=ax[1],\n", - " colorbar_label=\"Filtered Ref\",\n", - ")\n", - "ax[0].set_xlim([-50, 50])\n", - "ax[0].set_ylim([-50, 50])\n", - "ax[0].set_aspect(\"equal\", \"box\")\n", - "ax[1].set_aspect(\"equal\", \"box\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "74c8b2e9-dd8c-49de-a0fe-b49c9413f3a8", - "metadata": {}, - "outputs": [], - "source": [ - "list(radar.fields)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a4be687b-cf52-4e85-bd22-91b4054b7b20", - "metadata": {}, - "outputs": [], - "source": [ - "hydro = pyart.retrieve.hydroclass_semisupervised(\n", - " radar,\n", - " refl_field=\"corrected_reflectivity_horizontal\",\n", - " zdr_field=\"diff_reflectivity\",\n", - " kdp_field=\"diff_phase\",\n", - " rhv_field=\"copol_coeff\",\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "95d0f228-d999-4038-87e5-6cf06e622c9d", - "metadata": {}, - "outputs": [], - "source": [ - "radar.in" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6f232563-c9b1-4c75-875a-34a7bb51e9f5", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.5" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/examples/retrieve/plot_hydrometeor_class_x_band.py b/examples/retrieve/plot_hydrometeor_class_x_band.py new file mode 100644 index 0000000000..1841562fd8 --- /dev/null +++ b/examples/retrieve/plot_hydrometeor_class_x_band.py @@ -0,0 +1,227 @@ +""" +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="pyart_Carbone42") # ZDR +plt.xlim(-20, 20) +plt.ylim(-20, 20) + +ax3 = plt.subplot(1, 3, 3) +ax3 = display.plot("corrected_specific_diff_phase", cmap="pyart_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", +) + +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, +) + +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, +) + +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", + "Fuchsia", +] +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, 10, 0.9)) + cb.ax.set_yticklabels( + [ + "Drizzle", + "Rain", + "Ice Crystals", + "Aggregates", + "Wet Snow", + "Vertical Ice", + "LD Graupel", + "HD Graupel", + "Hail", + "Big Drops", + ] + ) + 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]) diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index a084597171..695605da32 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -609,6 +609,7 @@ def hydroclass_semisupervised( kdp_field=None, temp_field=None, hydro_field=None, + radar_freq=None, ): """ Classifies precipitation echoes following the approach by Besic et al @@ -634,6 +635,9 @@ def hydroclass_semisupervised( Output. Field name which represents the hydrometeor class field. A value of None will use the default field name as defined in the Py-ART configuration file. + radar_freq : str, optional + Radar frequency in Hertz (Hz) used for classification. + This parameter will be ignored, if the radar object has frequency information. Returns ------- @@ -650,8 +654,9 @@ def hydroclass_semisupervised( Notes ----- The default hydrometeor classification is valid for C-band radars. For X-band radars, - if frequency information is not present in the `radar.instrument_parameters`, a warning that the - algorithm is defaulting to the C band is printed. + if frequency information is not present in the `radar.instrument_parameters`, the user-supplied + `radar_freq` will be used with a warning. If both `radar.instrument_parameters` and + `radar_freq` parameter are missing, the algorithm defaults to the C band. If the radar frequency information is missing from the radar object, you can add it in `radar.instrument_parameters`, as follows: @@ -668,22 +673,20 @@ def hydroclass_semisupervised( # select the centroids as a function of frequency band if mass_centers is None: # assign coefficients according to radar frequency - if radar.instrument_parameters is not None: - if "frequency" in radar.instrument_parameters: - mass_centers = _get_mass_centers( - radar.instrument_parameters["frequency"]["data"][0] - ) - else: - mass_centers = _mass_centers_table()["C"] - warn( - "Radar frequency unknown. " - "Default coefficients for C band will be applied." - ) + if radar.instrument_parameters and "frequency" in radar.instrument_parameters: + frequency = radar.instrument_parameters["frequency"]["data"][0] + mass_centers = _get_mass_centers(frequency) + warn(f"Using radar frequency from instrument parameters: {frequency}") + elif radar_freq is not None: + mass_centers = _get_mass_centers(radar_freq) + warn( + f"Radar instrument parameters are empty. Using user-supplied radar frequency: {radar_freq}" + ) else: mass_centers = _mass_centers_table()["C"] warn( - "Radar instrument parameters is empty. So frequency is " - "unknown. Default coefficients for C band will be applied." + "Radar instrument parameters and radar_freq param are empty." + "So frequency is unknown. Default coefficients for C band will be applied." ) # parse the field parameters