Move functions from dp to util and add tests

wradlib/dp.py

@@ -54,7 +54,6 @@
 
 import numpy as np
 from scipy.interpolate import interp1d
-from scipy.signal import medfilt
 from scipy.stats import linregress
 from scipy.ndimage.filters import convolve1d
 from . import util as util
@@ -376,7 +375,7 @@ def unfold_phi(phidp, rho, width=5, copy=False):
     if copy:
         phidp = phidp.copy()
     rho = rho.reshape((-1, shape[-1]))
-    gradphi = gradient_from_smoothed(phidp)
+    gradphi = util.gradient_from_smoothed(phidp)
 
     beams, rs = phidp.shape
 
@@ -425,7 +424,7 @@ def unfold_phi_naive(phidp, rho, width=5, copy=False):
     if copy:
         phidp = phidp.copy()
     rho = rho.reshape((-1, shape[-1]))
-    gradphi = gradient_from_smoothed(phidp)
+    gradphi = util.gradient_from_smoothed(phidp)
 
     beams, rs = phidp.shape
 
@@ -558,36 +557,5 @@ def texture(data):
     return np.sqrt(num / xa_valid_count)
 
 
-# TO UTILS
-def medfilt_along_axis(x, N, axis=-1):
-    """Applies median filter smoothing on one axis of an N-dimensional array.
-    """
-    kernel_size = np.array(x.shape)
-    kernel_size[:] = 1
-    kernel_size[axis] = N
-    return medfilt(x, kernel_size)
-
-
-# TO UTILS
-def gradient_along_axis(x):
-    """Computes gradient along last axis of an N-dimensional array
-    """
-    axis = -1
-    newshape = np.array(x.shape)
-    newshape[axis] = 1
-    diff_begin = (x[..., 1] - x[..., 0]).reshape(newshape)
-    diff_end = (x[..., -1] - x[..., -2]).reshape(newshape)
-    diffs = ((x - np.roll(x, 2, axis)) / 2.)
-    diffs = np.append(diffs[..., 2:], diff_end, axis=axis)
-    return np.insert(diffs, 0, diff_begin, axis=axis)
-
-
-# TO UTILS
-def gradient_from_smoothed(x, N=5):
-    """Computes gradient of smoothed data along final axis of an array
-    """
-    return gradient_along_axis(medfilt_along_axis(x, N)).astype("f4")
-
-
 if __name__ == '__main__':
     print('wradlib: Calling module <dp> as main...')

wradlib/tests/test_util.py

@@ -108,6 +108,27 @@ def test_roll2d_polar(self):
         np.testing.assert_equal(result4[:, :-1],
                                 np.roll(data, -1, axis=1)[:, :-1])
 
+    def test_medfilt_along_axis(self):
+        x = np.arange(10).reshape((2, 5)).astype("f4")
+        shouldbe = np.array([[0., 1., 2., 3., 3.],
+                             [5., 6., 7., 8., 8.]])
+        result = util.medfilt_along_axis(x, 3)
+        np.testing.assert_allclose(result, shouldbe)
+
+    def test_gradient_along_axis(self):
+        x = np.arange(10).reshape((2, 5)).astype("f4") ** 2
+        result = util.gradient_along_axis(x)
+        shouldbe = np.array([[1., 11., 2., 4., 6., 7.],
+                             [1., 11., 12., 14., 16., 17.]])
+        np.testing.assert_allclose(result, shouldbe)
+
+    def test_gradient_from_smoothed(self):
+        x = np.arange(10).reshape((2, 5)).astype("f4") ** 2
+        result = util.gradient_from_smoothed(x)
+        shouldbe = np.array([[1., 11., 2., 1.5, 0., 0.],
+                             [1., 11., 12., 6.5, 0., 0.]])
+        np.testing.assert_allclose(result, shouldbe)
+
 
 class TestUtil(unittest.TestCase):
     def setUp(self):

wradlib/util.py

@@ -29,6 +29,7 @@
 import numpy as np
 from scipy.ndimage import filters
 from osgeo import ogr
+from scipy.signal import medfilt
 
 
 class OptionalModuleStub(object):
@@ -741,5 +742,33 @@ def calculate_polynomial(data, w):
     return poly
 
 
+def medfilt_along_axis(x, n, axis=-1):
+    """Applies median filter smoothing on one axis of an N-dimensional array.
+    """
+    kernel_size = np.array(x.shape)
+    kernel_size[:] = 1
+    kernel_size[axis] = n
+    return medfilt(x, kernel_size)
+
+
+def gradient_along_axis(x):
+    """Computes gradient along last axis of an N-dimensional array
+    """
+    axis = -1
+    newshape = np.array(x.shape)
+    newshape[axis] = 1
+    diff_begin = (x[..., 1] - x[..., 0]).reshape(newshape)
+    diff_end = (x[..., -1] - x[..., -2]).reshape(newshape)
+    diffs = ((x - np.roll(x, 2, axis)) / 2.)
+    diffs = np.append(diffs[..., 2:], diff_end, axis=axis)
+    return np.insert(diffs, 0, diff_begin, axis=axis)
+
+
+def gradient_from_smoothed(x, n=5):
+    """Computes gradient of smoothed data along final axis of an array
+    """
+    return gradient_along_axis(medfilt_along_axis(x, n)).astype("f4")
+
+
 if __name__ == '__main__':
     print('wradlib: Calling module <util> as main...')