Interpolation switched to RegularGridInterpolator; works for ND datasets

scipy.interpolate.interp2d is deprecated

VERSION

@@ -1 +1 @@
-0.9.0.dev23
+0.9.0.dev24

aspecd/processing.py

@@ -2638,16 +2638,9 @@ class Interpolation(SingleProcessingStep):
     * `<https://stackoverflow.com/a/32763635>`_
 
 
-    .. important::
-        Currently, interpolation works *only* for **1D and 2D** datasets,
-        not for higher-dimensional datasets. This may, however, change in
-        the future.
-
     .. todo::
         * Make type of interpolation controllable
 
-        * Check for ways to make it work with ND, N>2
-
 
     Attributes
     ----------
@@ -2727,6 +2720,14 @@ class Interpolation(SingleProcessingStep):
 
     .. versionadded:: 0.2
 
+    .. versionchanged:: 0.8.3
+        Interpolation for *N*\ D datasets with arbitrary dimension *N*
+
+    .. versionchanged:: 0.8.3
+        Change interpolation method for 2D data from deprecated
+        :class:`scipy.interpolate.interp2d` to
+        :class:`scipy.interpolate.RegularGridInterpolator`
+
     """
 
     def __init__(self):
@@ -2738,27 +2739,6 @@ def __init__(self):
         self.parameters["unit"] = "index"
         self._axis_values = []
 
-    @staticmethod
-    def applicable(dataset):
-        """
-        Check whether processing step is applicable to the given dataset.
-
-        Interpolation is currently only applicable to datasets with one- and
-        two-dimensional data.
-
-        Parameters
-        ----------
-        dataset : :class:`aspecd.dataset.Dataset`
-            dataset to check
-
-        Returns
-        -------
-        applicable : :class:`bool`
-            `True` if successful, `False` otherwise.
-
-        """
-        return len(dataset.data.axes) <= 3
-
     def _sanitise_parameters(self):
         if not self.parameters["range"]:
             raise ValueError('No range provided to interpolate for')
@@ -2777,18 +2757,14 @@ def _sanitise_parameters(self):
 
     def _perform_task(self):
         self._get_axis_values()
-        if self.dataset.data.data.ndim == 1:
-            interp = interpolate.interp1d(self.dataset.data.axes[0].values,
-                                          self.dataset.data.data)
-            self.dataset.data.data = interp(self._axis_values[0])
-        elif self.dataset.data.data.ndim == 2:
-            # Note: interp2d uses Cartesian indexing (x,y => col, row),
-            #       not matrix indexing (row, col)
-            interp = interpolate.interp2d(self.dataset.data.axes[1].values,
-                                          self.dataset.data.axes[0].values,
-                                          self.dataset.data.data)
-            self.dataset.data.data = interp(self._axis_values[1],
-                                            self._axis_values[0])
+        interp = interpolate.RegularGridInterpolator(
+            [x.values for x in self.dataset.data.axes[:-1]],
+            self.dataset.data.data
+        )
+        grid = np.meshgrid(*self._axis_values, indexing='ij')
+        test_points = np.array([x.ravel() for x in grid]).T
+        shape = [len(x) for x in self._axis_values]
+        self.dataset.data.data = interp(test_points).reshape(shape)
         for dim in range(self.dataset.data.data.ndim):
             self.dataset.data.axes[dim].values = self._axis_values[dim]
 

docs/changelog.rst

@@ -11,6 +11,31 @@ Version 0.9.0
 *Not yet released*
 
 
+Version 0.8.3
+=============
+
+*Not yet released*
+
+Fixes
+-----
+
+* Exporter tasks (:class:`aspecd.tasks.ExportTask`) automatically save datasets with default name if no target is provided.
+* Correct setting of contour plot properties with newer versions of Matplotlib
+
+
+Changes
+-------
+
+* :class:`aspecd.processing.Interpolation` changed interpolation method for 2D data from deprecated :class:`scipy.interpolate.interp2d` to :class:`scipy.interpolate.RegularGridInterpolator`
+
+
+New features
+------------
+
+* :class:`aspecd.processing.Interpolation` works for *N*\ D datasets with arbitrary dimension *N*
+* :class:`aspecd.tasks.Recipe` with new setting ``autosave_datasets`` (default: ``True``)
+
+
 Version 0.8.2
 =============
 

requirements.txt

@@ -1,5 +1,5 @@
 numpy
-scipy
+scipy>=0.14
 matplotlib
 jinja2>=3.0
 oyaml

setup.py

@@ -49,7 +49,7 @@ def read(fname):
         'jinja2>=3.0',
         'matplotlib',
         'numpy',
-        'scipy',
+        'scipy>=0.14',
         'oyaml',
         'asdf',
         'bibrecord',

tests/test_processing.py

@@ -1796,6 +1796,13 @@ def setUp(self):
             np.reshape(np.tile(np.linspace(0, 5, 11), 21), (21, 11))
         self.dataset2d.data.axes[0].values = np.linspace(30, 40, 21)
         self.dataset2d.data.axes[1].values = np.linspace(5, 15, 11)
+        self.dataset3d = aspecd.dataset.Dataset()
+        self.dataset3d.data.data = \
+            np.reshape(np.tile(np.linspace(0, 5, 11), [21, 16]).ravel(),
+                       (21, 11, 16))
+        self.dataset3d.data.axes[0].values = np.linspace(30, 40, 21)
+        self.dataset3d.data.axes[1].values = np.linspace(5, 15, 11)
+        self.dataset3d.data.axes[2].values = np.linspace(2.5, 10, 16)
 
     def test_instantiate_class(self):
         pass
@@ -1806,12 +1813,6 @@ def test_has_appropriate_description(self):
     def test_is_undoable(self):
         self.assertTrue(self.processing.undoable)
 
-    def test_with_dataset_with_more_than_2d_raises(self):
-        dataset = aspecd.dataset.Dataset()
-        dataset.data.data = np.random.random([11, 11, 11])
-        with self.assertRaises(aspecd.exceptions.NotApplicableToDatasetError):
-            dataset.process(self.processing)
-
     def test_interpolate_without_range_raises(self):
         self.processing.parameters["npoints"] = 21
         with self.assertRaisesRegex(ValueError, 'No range provided'):
@@ -1912,7 +1913,7 @@ def test_interpolate_2d_data_interpolates_data(self):
         self.processing.parameters["npoints"] = [41, 21]
         self.dataset2d.process(self.processing)
         self.assertListEqual(list(np.linspace(0, 5, 21)),
-                             list(self.dataset2d.data.data[1, :]))
+                             list(self.dataset2d.data.data[0, :]))
 
     def test_interpolate_2d_data_interpolates_axis(self):
         self.processing.parameters["range"] = [[0, 20], [0, 10]]
@@ -1921,6 +1922,28 @@ def test_interpolate_2d_data_interpolates_axis(self):
         self.assertListEqual(list(np.linspace(30, 40, 41)),
                              list(self.dataset2d.data.axes[0].values))
 
+    def test_interpolate_3d_data_with_missing_range_raises(self):
+        self.processing.parameters["range"] = [[0, 20], [0, 10]]
+        self.processing.parameters["npoints"] = [41, 21, 21]
+        with self.assertRaisesRegex(IndexError, 'List of ranges does not fit '
+                                                'data dimensions'):
+            self.dataset3d.process(self.processing)
+
+    def test_interpolate_3d_data_with_missing_npoints_raises(self):
+        self.processing.parameters["range"] = [[0, 20], [0, 10], [5, 15]]
+        self.processing.parameters["npoints"] = [41, 21]
+        with self.assertRaisesRegex(IndexError, 'List of npoints does not fit '
+                                                'data dimensions'):
+            self.dataset3d.process(self.processing)
+
+    def test_interpolate_3d_data_interpolates_data(self):
+        self.processing.parameters["range"] = [[0, 20], [0, 10], [0, 15]]
+        self.processing.parameters["npoints"] = [41, 21, 16]
+        old_data = self.dataset3d.data.data[0, :, 0]
+        self.dataset3d.process(self.processing)
+        self.assertListEqual(list(old_data),
+                             list(self.dataset3d.data.data[0, ::2, 0]))
+
 
 class TestFiltering(unittest.TestCase):