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 <[email protected]>
Co-authored-by: Max Grover <[email protected]>

continuous_integration/environment-ci.yml

@@ -28,3 +28,4 @@ dependencies:
     - cibuildwheel
     - pooch
     - versioneer
+    - git+https://github.com/openradar/open-radar-data

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

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 ([email protected])
+# 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

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 ([email protected])
+# 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')