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')