From befcd363e3536eadbbf06d76542fb08711cf8799 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Tue, 8 Nov 2022 15:49:42 +0100 Subject: [PATCH] DOC: Added two pyart examples (VAD and hydrometeors) (#1317) * Added two pyart examples (VAD and hydrometeors) * Added two pyart examples (VAD and hydrometeors) * Added two pyart examples (VAD and hydrometeors) * ENH: Add link to open radar data github repo. * STY: PEP8 fixes for both examples. * Update examples/retrieve/plot_hydrometeor.py * Apply suggestions from mgrover1 * add open-radar-data to docs environment Co-authored-by: zssherman Co-authored-by: Max Grover --- continuous_integration/environment-ci.yml | 1 + doc/environment.yml | 1 + examples/retrieve/plot_hydrometeor.py | 84 +++++++++++++++++++++++ examples/retrieve/plot_vad.py | 49 +++++++++++++ 4 files changed, 135 insertions(+) create mode 100644 examples/retrieve/plot_hydrometeor.py create mode 100644 examples/retrieve/plot_vad.py diff --git a/continuous_integration/environment-ci.yml b/continuous_integration/environment-ci.yml index 0b61cb07ab..e0737dd0ab 100644 --- a/continuous_integration/environment-ci.yml +++ b/continuous_integration/environment-ci.yml @@ -28,3 +28,4 @@ dependencies: - cibuildwheel - pooch - versioneer + - git+https://github.com/openradar/open-radar-data diff --git a/doc/environment.yml b/doc/environment.yml index ab191aeeeb..11318c57e2 100644 --- a/doc/environment.yml +++ b/doc/environment.yml @@ -27,6 +27,7 @@ dependencies: - shapely<1.8.3 - pip - pip: + - git+https://github.com/openradar/open-radar-data - pydata-sphinx-theme<0.9.0 - sphinx_gallery - sphinx-copybutton diff --git a/examples/retrieve/plot_hydrometeor.py b/examples/retrieve/plot_hydrometeor.py new file mode 100644 index 0000000000..425204dcbd --- /dev/null +++ b/examples/retrieve/plot_hydrometeor.py @@ -0,0 +1,84 @@ + +""" +========================================= +Calculate and Plot hydrometeor classification +========================================= + +Calculates a hydrometeor classification and displays the results +""" + +# Author: Daniel Wolfensberger (daniel.wolfensberger@meteoswiss.ch) +# License: BSD 3 clause + +import numpy as np +import matplotlib as mpl +import matplotlib.pyplot as plt + +import pyart +from open_radar_data import DATASETS + +# Read in a sample file +filename = DATASETS.fetch('MLL2217907250U.003.nc') +radar = pyart.io.read_cfradial(filename) + +# Read temperature preinterpolated from NWP model +filename = DATASETS.fetch('20220628072500_savevol_COSMO_LOOKUP_TEMP.nc') +nwp_temp = pyart.io.read_cfradial(filename) + +# Add temperature to radar object as new field +radar.add_field('temperature', nwp_temp.fields['temperature']) + +# Compute attenuation +out = pyart.correct.calculate_attenuation_zphi( + radar, + phidp_field='uncorrected_differential_phase', + temp_field='temperature', + temp_ref='temperature') +spec_at, pia, cor_z, spec_diff_at, pida, cor_zdr = out +radar.add_field('corrected_reflectivity', cor_z) +radar.add_field('corrected_differential_reflectivity', cor_zdr) +radar.add_field('specific_attenuation', spec_at) + +# Compute KDP +kdp, _, _ = pyart.retrieve.kdp_maesaka( + radar, psidp_field='uncorrected_differential_phase') +radar.add_field('specific_differential_phase', kdp) + +# Compute hydrometeor classification +hydro = pyart.retrieve.hydroclass_semisupervised( + radar, + refl_field='corrected_reflectivity', + zdr_field='corrected_differential_reflectivity', + kdp_field='specific_differential_phase', + rhv_field='uncorrected_cross_correlation_ratio', + temp_field='temperature') + +radar.add_field('radar_echo_classification', hydro) + +# Display hydrometeor classification with categorical colormap +fig, ax = plt.subplots(1,1, figsize=(6, 6)) +display = pyart.graph.RadarDisplay(radar) + +labels = ['NC','AG', 'CR', 'LR', 'RP', 'RN', 'VI', 'WS', 'MH', 'IH/HDG'] +ticks = np.arange(len(labels)) +boundaries = np.arange(-0.5, len(labels)) +norm = mpl.colors.BoundaryNorm(boundaries, 256) + +cax = display.plot_ppi('radar_echo_classification', 0, ax=ax, + norm=norm, ticks=ticks, ticklabs=labels) + +ax.set_xlim([-50,50]) +ax.set_ylim([-50,50]) +ax.set_aspect('equal', 'box') + +# For info +# NC = not classified +# AG = aggregates +# CR = ice crystals +# LR = light rain +# RP = rimed particles +# RN = rain +# VI = vertically oriented ice +# WS = wet snow +# MH = melting hail +# IH/HDG = dry hail / high density graupel diff --git a/examples/retrieve/plot_vad.py b/examples/retrieve/plot_vad.py new file mode 100644 index 0000000000..8806094cc5 --- /dev/null +++ b/examples/retrieve/plot_vad.py @@ -0,0 +1,49 @@ +""" +========================================= +Calculate and Plot VAD profile +========================================= + +Calculates a VAD and plots a vertical profile of wind +""" + +# Author: Daniel Wolfensberger (daniel.wolfensberger@meteoswiss.ch) +# License: BSD 3 clause + +import numpy as np +import matplotlib.pyplot as plt + +import pyart +from open_radar_data import DATASETS + +# Read in a sample file +filename = DATASETS.fetch('MLA2119412050U.nc') +radar = pyart.io.read_cfradial(filename) + +# Loop on all sweeps and compute VAD +zlevels = np.arange(100, 5000, 100) # height above radar +u_allsweeps = [] +v_allsweeps = [] + +for idx in range(radar.nsweeps): + radar_1sweep = radar.extract_sweeps([idx]) + vad = pyart.retrieve.vad_browning( + radar_1sweep, + 'corrected_velocity', + z_want=zlevels) + u_allsweeps.append(vad.u_wind) + v_allsweeps.append(vad.v_wind) + +# Average U and V over all sweeps and compute magnitude and angle +u_avg = np.nanmean(np.array(u_allsweeps), axis=0) +v_avg = np.nanmean(np.array(v_allsweeps), axis=0) +orientation = np.rad2deg(np.arctan2(-u_avg, -v_avg))%360 +speed = np.sqrt(u_avg**2 + v_avg**2) + +# Display vertical profile of wind +fig, ax = plt.subplots(1, 2, sharey=True) +ax[0].plot(speed*2, zlevels+radar.altitude['data']) +ax[1].plot(orientation, zlevels+radar.altitude['data']) +ax[0].set_xlabel('Wind speed [m/s]') +ax[1].set_xlabel('Wind direction [deg]') +ax[0].set_ylabel('Altitude [m]') +fig.suptitle('Wind profile obtained from VAD')