ADD: Add in tests for KDP/PHI DP functions (#1184)

* rename compare kdp code

* remove old matplotlib artificat

* don't use parallel processing

* reduce size of test array

* ensure even windsize

* only test first kdp function

* add in second kdp method

* add in last kdp retrieval

* add image test for kdp processing

* remove plt.show() and return fig

* add test at beginning of function

continuous_integration/environment-ci.yml

@@ -14,6 +14,7 @@ dependencies:
   - cvxopt
   - xarray>=0.21.1
   - pytest-cov
+  - pytest-mpl
   - coveralls
   - flake8
   - pytest

pyart/retrieve/tests/test_compare_kdp_proc.py

@@ -0,0 +1,85 @@
+"""
+A script for comparing the e Vulpiani, Maesaka and Schneedbeli methods
+for KDP estimation with sample data.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import pytest
+
+from pyart.retrieve.kdp_proc import kdp_maesaka, kdp_vulpiani, kdp_schneebeli
+from pyart.testing import sample_objects
+from pyart.config import get_field_name
+
+
+def _make_real_psidp_radar():
+    """
+    Create single-ray radar with linear differential phase profile with
+    specified slope.
+    ---
+    Returns
+    -------
+    radar : Radar
+            PyART radar instance with differential phase profile in deg.
+    """
+    psidp = np.array([[-2.33313751e+00, 1.80617523e+00, 7.17742920e-01,
+                     1.82811661e+01, 1.89352417e+01, 1.67904205e+01]])
+    psidp = np.ma.array(psidp)
+    radar = sample_objects.make_empty_ppi_radar(len(psidp[0]), 1, 1)
+    psidp_dict = {
+                  'data': psidp,
+                  }
+    radar.add_field(get_field_name('differential_phase'), psidp_dict)
+    # Define real ranges
+    radar.range['data'] = 75*np.arange(0, len(psidp[0]))
+    return radar
+
+
+@pytest.mark.mpl_image_compare(tolerance=30)
+def test_compare_kdp_estimation_methods():
+    # Get profile of noisy psidp
+    prof_psidp = _make_real_psidp_radar()
+    # Maesaka method
+    kdp_mae, phidpf_mae, phidp_mae = kdp_maesaka(prof_psidp,
+                                                 maxiter=1000,
+                                                 check_outliers=False,
+                                                 parallel=False)
+
+    # Vulpiani method (note windsize is just a guess here..)
+    kdp_vulp, phidp_vulp = kdp_vulpiani(prof_psidp,
+                                        windsize=2,
+                                        n_iter=20,
+                                        parallel=False,
+                                        band='X')
+    # Kalman filter method
+    kdp_schnee, kdp_std_schnee, phidp_schnee = kdp_schneebeli(prof_psidp,
+                                                              parallel=False,
+                                                              band='X')
+    # Create figure
+    fig = plt.figure(figsize=(10, 10))
+    plt.subplot(2, 1, 1)
+    plt.grid(True)
+    plt.title('Kdp estimation')
+    ranges = prof_psidp.range['data']
+    plt.plot(ranges, kdp_mae['data'][0])
+    plt.plot(ranges, kdp_vulp['data'][0])
+    plt.plot(ranges, kdp_schnee['data'][0])
+    plt.xlabel('Range [m]')
+    plt.ylabel('Kdp [deg/km]')
+    plt.legend(['Maesaka', 'Vulpiani', 'Schneebeli'], loc=0)
+    plt.subplot(2, 1, 2)
+    plt.grid(True)
+    plt.title('Reconstructed Phidp')
+    ranges = prof_psidp.range['data']
+    phidp_mae = 0.5 * (phidp_mae['data'][0] + phidp_mae['data'][0])
+    plt.plot(ranges, phidp_mae)
+    plt.plot(ranges, phidp_vulp['data'][0])
+    plt.plot(ranges, phidp_schnee['data'][0])
+    plt.plot(ranges, prof_psidp.fields['differential_phase']['data'][0])
+    plt.xlabel('Range [m]')
+    plt.ylabel('Diff. phase [deg]')
+    plt.legend(['Maesaka', 'Vulpiani', 'Schneebeli', 'Real Psidp'], loc=0)
+    try:
+        return fig
+    finally:
+        plt.close(fig)