# [Scipy-svn] r5156 - trunk/scipy/interpolate

scipy-svn@scip... scipy-svn@scip...
Fri Nov 21 15:08:09 CST 2008

```Author: ptvirtan
Date: 2008-11-21 15:07:57 -0600 (Fri, 21 Nov 2008)
New Revision: 5156

Modified:
trunk/scipy/interpolate/interpolate.py
Log:
Implement kind='nearest' interpolation, code adapted from #305

Modified: trunk/scipy/interpolate/interpolate.py
===================================================================
--- trunk/scipy/interpolate/interpolate.py	2008-11-21 10:24:29 UTC (rev 5155)
+++ trunk/scipy/interpolate/interpolate.py	2008-11-21 21:07:57 UTC (rev 5156)
@@ -8,7 +8,7 @@

from numpy import shape, sometrue, rank, array, transpose, searchsorted, \
ones, logical_or, atleast_1d, atleast_2d, meshgrid, ravel, \
-                  dot, poly1d, asarray
+                  dot, poly1d, asarray, intp
import numpy as np
import scipy.special as spec
import math
@@ -198,14 +198,14 @@
self.bounds_error = bounds_error
self.fill_value = fill_value

-        if kind in ['zero', 'slinear', 'quadratic', 'cubic', 'nearest']:
+        if kind in ['zero', 'slinear', 'quadratic', 'cubic']:
order = {'nearest':0, 'zero':0,'slinear':1,
kind = 'spline'
elif isinstance(kind, int):
order = kind
kind = 'spline'
-        elif kind != 'linear':
+        elif kind not in ('linear', 'nearest'):
raise NotImplementedError("%s is unsupported: Use fitpack "\
"routines for other types." % kind)
x = array(x, copy=self.copy)
@@ -220,7 +220,7 @@
self.axis = axis % len(y.shape)
self._kind = kind

-        if kind == 'linear':
+        if kind in ('linear', 'nearest'):
# Make a "view" of the y array that is rotated to the interpolation
# axis.
axes = range(y.ndim)
@@ -229,7 +229,11 @@
oriented_y = y.transpose(axes)
minval = 2
len_y = oriented_y.shape[-1]
-            self._call = self._call_linear
+            if kind == 'linear':
+                self._call = self._call_linear
+            elif kind == 'nearest':
+                self.x_bds = (x[1:] + x[:-1]) / 2.0
+                self._call = self._call_nearest
else:
axes = range(y.ndim)
del axes[self.axis]
@@ -242,10 +246,10 @@

len_x = len(x)
if len_x != len_y:
-            raise ValueError("x and y arrays must be equal in length along"
-            "interpolation axis.")
+            raise ValueError("x and y arrays must be equal in length along "
+                             "interpolation axis.")
if len_x < minval:
-            raise ValueError("x and y arrays must have at " \
+            raise ValueError("x and y arrays must have at "
"least %d entries" % minval)
self.x = x
self.y = oriented_y
@@ -280,6 +284,23 @@

return y_new

+    def _call_nearest(self, x_new):
+        """ Find nearest neighbour interpolated y_new = f(x_new)."""
+
+        # 2. Find where in the averaged data the values to interpolate
+        #    would be inserted.
+        #    Note: use side='left' (right) to searchsorted() to define the
+        #    halfway point to be nearest to the left (right) neighbour
+        x_new_indices = searchsorted(self.x_bds, x_new, side='left')
+
+        # 3. Clip x_new_indices so that they are within the range of x indices.
+        x_new_indices = x_new_indices.clip(0,  len(self.x)-1).astype(intp)
+
+        # 4. Calculate the actual value for each entry in x_new.
+        y_new = self.y[..., x_new_indices]
+
+        return y_new
+
def _call_spline(self, x_new):
x_new =np.asarray(x_new)
result = spleval(self._spline,x_new.ravel())
@@ -327,7 +348,7 @@
else:
y_new[...] = self.fill_value
return asarray(y_new)
-        elif self._kind == 'linear':
+        elif self._kind in ('linear', 'nearest'):
y_new[..., out_of_bounds] = self.fill_value
axes = range(ny - nx)
axes[self.axis:self.axis] = range(ny - nx, ny)

```