Added a test for 0 gradient on the linear extrapolation.

Author: pelson

stratify/_vinterp.pyx

@@ -368,15 +368,21 @@ cdef class LinearExtrapolator(Extrapolator):
         cdef unsigned int m = fz_src.shape[0]
         cdef unsigned int n_src_pts = fz_src.shape[1]
         cdef unsigned int p0, p1, i
-        cdef double frac
+        cdef double frac, z_step
 
         if direction < 0:
             p0, p1 = 0, 1
         else:
             p0, p1 = n_src_pts - 2, n_src_pts - 1
 
-        frac = ((level - z_src[p0]) /
-                (z_src[p1] - z_src[p0]))
+        # Compute the normalised distance of the target point between p0 and p1
+        z_step = z_src[p1] - z_src[p0]
+        if z_step == 0:
+            # If there is nothing between the last two points then we
+            # extrapolate using a 0 gradient.
+            frac = 0
+        else:
+            frac = ((level - z_src[p0]) / z_step)
 
         for i in range(m):
            fz_target[i] = fz_src[i, p0] + frac * (fz_src[i, p1] - fz_src[i, p0])

stratify/tests/test_vinterp.py

@@ -148,6 +148,12 @@ def test_on_the_mark(self):
         assert_array_equal(self.interpolate([0, 1, 2, 3, 4]),
                            [0, 10, 20, 30, 40])
 
+    def test_zero_gradient(self):
+        assert_array_equal(
+            stratify.interpolate([1], [0, 1, 1, 2], [10, 20, 30, 40],
+                                 interpolation='linear'),
+            [20])
+
     def test_inbetween(self):
         assert_array_equal(self.interpolate([0.5, 1.25, 2.5, 3.75]),
                            [5, 12.5, 25, 37.5])
@@ -258,6 +264,12 @@ def test_below(self):
     def test_above(self):
         assert_array_almost_equal(self.interpolate([15.123]), [151.23])
 
+    def test_zero_gradient(self):
+        assert_array_almost_equal(
+            stratify.interpolate([2], [0, 0], [1, 1],
+                                 extrapolation='linear'),
+            [1])
+
     def test_npts(self):
         interpolation = IndexInterpolator()
         extrapolation = stratify.EXTRAPOLATE_LINEAR
@@ -369,7 +381,7 @@ class Test_interpolate(unittest.TestCase):
     def test_target_z_3d_axis_0(self):
         z_target = z_source = f_source = np.arange(3) * np.ones([4, 2, 3])
         result= vinterp.interpolate(z_target, z_source, f_source,
-                                    axis=0, extrapolation='linear')
+                                    extrapolation='linear')
         assert_array_equal(result, f_source)