Adds tests showing out-of-range data returning from scipy interpolate…

improver/utilities/interpolation.py

@@ -53,7 +53,7 @@ def interpolate_missing_data(
         was possible to fill these in.
     """
     if valid_points is None:
-        valid_points = np.full_like(data, True, dtype=np.bool)
+        valid_points = np.full_like(data, True, dtype=bool)
 
     # Interpolate linearly across the remaining points
     index = ~np.isnan(data)
@@ -187,6 +187,8 @@ def process(
                 out=np.full(cslice.shape, np.nan),
                 where=valid_points,
             )
+            min_difference = np.nanmin(difference_field)
+            max_difference = np.nanmax(difference_field)
             interpolated_difference = interpolate_missing_data(
                 difference_field, valid_points=valid_points
             )
@@ -198,6 +200,11 @@ def process(
                 interpolated_difference = interpolate_missing_data(
                     difference_field, valid_points=~remain_invalid, method="nearest"
                 )
+            # It is possible for the interpolated differences to be outside of the range
+            # of the source data (machine-precision). Enforce the original data range.
+            interpolated_difference = np.clip(
+                interpolated_difference, min_difference, max_difference
+            )
 
             result = cslice.copy()
             result.data[invalid_points] = (

improver_tests/utilities/test_InterpolateUsingDifference.py

@@ -8,7 +8,7 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_array_almost_equal, assert_array_equal
 
 from improver.synthetic_data.set_up_test_cubes import (
     add_coordinate,
@@ -286,6 +286,41 @@ def test_convert_units(self):
         self.assertEqual(result.coords(), self.sleet_rain.coords())
         self.assertEqual(result.metadata, self.sleet_rain.metadata)
 
+    def test_range_enforcement(self):
+        """Test interpolation on a case where the result is known to be outside of the
+        input data range."""
+        data = np.zeros(
+            (18, 18), dtype=np.float32
+        )  # The smallest array where this behaviour has been found
+        data[1:-1, 1:-1] = np.nan
+        data[0, 4] = 100
+        sleet_rain = np.ma.masked_invalid(data)
+
+        self.snow_sleet = set_up_variable_cube(
+            np.zeros_like(data),
+            name="altitude_of_snow_falling_level",
+            units="m",
+            spatial_grid="equalarea",
+        )
+        self.sleet_rain = set_up_variable_cube(
+            sleet_rain,
+            name="altitude_of_rain_falling_level",
+            units="m",
+            spatial_grid="equalarea",
+        )
+
+        expected = np.zeros_like(data)
+        expected[0, 4] = 100
+        expected[1, 3] = 75
+        expected[2, 2] = 50
+        expected[3, 1] = 25
+
+        result = InterpolateUsingDifference().process(self.sleet_rain, self.snow_sleet)
+
+        assert_array_almost_equal(result.data, expected, decimal=5)
+        assert (result.data >= np.nanmin(self.sleet_rain.data)).all()
+        assert (result.data <= np.nanmax(self.sleet_rain.data)).all()
+
 
 if __name__ == "__main__":
     unittest.main()

improver_tests/utilities/test_interpolate.py

@@ -43,6 +43,30 @@ def setUp(self):
         self.valid_data_for_limit_test = np.full((5, 5), True)
         self.valid_data_for_limit_test[:, 1:4] = False
 
+    def test_mostly_zeros(self):
+        """Test when all-but-one of the points around the missing data are the same.
+        The point of this test is to highlight a case where values outside of the max:min
+        range of the input can be found, if the test tolerance is sufficiently tight.
+        If this test fails with a newer version of Scipy, then the enforcement of this range
+        in improver.utilitiess.interpolation.InterpolateUsingDifference needs revisiting."""
+        data = np.zeros(
+            (18, 18)
+        )  # The smallest array where this behaviour has been found
+        data[1:-1, 1:-1] = np.nan
+        data[0, 4] = 100
+        expected = np.zeros_like(data)
+        expected[0, 4] = 100
+        expected[1, 3] = 75
+        expected[2, 2] = 50
+        expected[3, 1] = 25
+        expected[2, 3] = -1.11022302e-14
+        expected[3, 2] = -4.44089210e-14
+        expected[4, 1] = -2.22044605e-14
+
+        data_updated = interpolate_missing_data(data)
+
+        self.assertArrayAlmostEqual(data_updated, expected, decimal=21)
+
     def test_basic_linear(self):
         """Test when all the points around the missing data are the same."""
         data = np.ones((3, 3))