[Scipy-svn] r6657 - in trunk/scipy/interpolate: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Aug 31 17:58:55 EDT 2010
Author: ptvirtan
Date: 2010-08-31 16:58:55 -0500 (Tue, 31 Aug 2010)
New Revision: 6657
Added:
trunk/scipy/interpolate/griddatand.py
trunk/scipy/interpolate/tests/test_griddatand.py
Modified:
trunk/scipy/interpolate/__init__.py
Log:
interpolate: add griddata, convenience function for N-d interpolation
Add a convenience function for easy N-d interpolation.
Modified: trunk/scipy/interpolate/__init__.py
===================================================================
--- trunk/scipy/interpolate/__init__.py 2010-08-31 21:58:25 UTC (rev 6656)
+++ trunk/scipy/interpolate/__init__.py 2010-08-31 21:58:55 UTC (rev 6657)
@@ -14,6 +14,8 @@
from polyint import *
+from griddatand import *
+
__all__ = filter(lambda s:not s.startswith('_'),dir())
from numpy.testing import Tester
test = Tester().test
Added: trunk/scipy/interpolate/griddatand.py
===================================================================
--- trunk/scipy/interpolate/griddatand.py (rev 0)
+++ trunk/scipy/interpolate/griddatand.py 2010-08-31 21:58:55 UTC (rev 6657)
@@ -0,0 +1,172 @@
+"""
+Convenience interface to N-D interpolation
+
+.. versionadded:: 0.9
+
+"""
+
+import numpy as np
+from interpnd import LinearNDInterpolator, NDInterpolatorBase, \
+ CloughTocher2DInterpolator, _ndim_coords_from_arrays
+from scipy.spatial import cKDTree
+
+__all__ = ['griddata', 'NearestNDInterpolator', 'LinearNDInterpolator',
+ 'CloughTocher2DInterpolator']
+
+#------------------------------------------------------------------------------
+# Nearest-neighbour interpolation
+#------------------------------------------------------------------------------
+
+class NearestNDInterpolator(NDInterpolatorBase):
+ """
+ NearestNDInterpolator(points, values)
+
+ Nearest-neighbour interpolation in N dimensions.
+
+ .. versionadded:: 0.9
+
+ Parameters
+ ----------
+ points : ndarray of floats, shape (npoints, ndims)
+ Data point coordinates.
+ values : ndarray of float or complex, shape (npoints, ...)
+ Data values.
+
+ Notes
+ -----
+ Uses ``scipy.spatial.cKDTree``
+
+ """
+
+ def __init__(self, x, y):
+ x = _ndim_coords_from_arrays(x)
+ self._check_init_shape(x, y)
+ self.tree = cKDTree(x)
+ self.points = x
+ self.values = y
+
+ def __call__(self, xi):
+ """
+ Evaluate interpolator at given points.
+
+ Parameters
+ ----------
+ xi : ndarray of float, shape (..., ndim)
+ Points where to interpolate data at.
+
+ """
+ xi = self._check_call_shape(xi)
+ dist, i = self.tree.query(xi)
+ return self.values[i]
+
+
+#------------------------------------------------------------------------------
+# Convenience interface function
+#------------------------------------------------------------------------------
+
+def griddata(points, values, xi, method='linear'):
+ """
+ griddata(points, values, xi, method='linear')
+
+ Interpolate unstructured N-dimensional data.
+
+ .. versionadded:: 0.9
+
+ Parameters
+ ----------
+ points : ndarray of floats, shape (npoints, ndims)
+ Data point coordinates. Can either be a ndarray of
+ size (npoints, ndim), or a tuple of `ndim` arrays.
+ values : ndarray of float or complex, shape (npoints, ...)
+ Data values.
+ xi : ndarray of float, shape (..., ndim)
+ Points where to interpolate data at.
+
+ method : {'linear', 'nearest', 'cubic'}
+ Method of interpolation. One of
+
+ - ``nearest``: return the value at the data point closest to
+ the point of interpolation. See `NearestNDInterpolator` for
+ more details.
+
+ - ``linear``: tesselate the input point set to n-dimensional
+ simplices, and interpolate linearly on each simplex. See
+ `LinearNDInterpolator` for more details.
+
+ - ``cubic`` (1-D): return the value detemined from a cubic
+ spline.
+
+ - ``cubic`` (2-D): return the value determined from a
+ piecewise cubic, continuously differentiable (C1), and
+ approximately curvature-minimizing polynomial surface. See
+ `CloughTocher2DInterpolator` for more details.
+
+
+ Examples
+ --------
+
+ Suppose we want to interpolate the 2-D function
+
+ >>> def func(x, y):
+ >>> return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
+
+ on a grid in [0, 1]x[0, 1]
+
+ >>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]
+
+ but we only know its values at 1000 data points:
+
+ >>> points = np.random.rand(1000, 2)
+ >>> values = func(points[:,0], points[:,1])
+
+ This can be done with `griddata` -- below we try out all of the
+ interpolation methods:
+
+ >>> from scipy.interpolate import griddata
+ >>> grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
+ >>> grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
+ >>> grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')
+
+ One can see that the exact result is reproduced by all of the
+ methods to some degree, but for this smooth function the piecewise
+ cubic interpolant gives the best results:
+
+ >>> import matplotlib.pyplot as plt
+ >>> plt.subplot(221)
+ >>> plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
+ >>> plt.plot(points[:,0], points[:,1], 'k.', ms=1)
+ >>> plt.title('Original')
+ >>> plt.subplot(222)
+ >>> plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
+ >>> plt.title('Nearest')
+ >>> plt.subplot(223)
+ >>> plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
+ >>> plt.title('Linear')
+ >>> plt.subplot(224)
+ >>> plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
+ >>> plt.title('Cubic')
+ >>> plt.gcf().set_size_inches(6, 6)
+ >>> plt.show()
+
+ """
+
+ points = _ndim_coords_from_arrays(points)
+ xi = _ndim_coords_from_arrays(xi)
+
+ ndim = points.shape[-1]
+
+ if ndim == 1 and method in ('nearest', 'linear', 'cubic'):
+ ip = interp1d(points, values, kind=method, axis=0, bounds_error=False)
+ return ip(xi)
+ elif method == 'nearest':
+ ip = NearestNDInterpolator(points, values)
+ return ip(xi)
+ elif method == 'linear':
+ ip = LinearNDInterpolator(points, values)
+ return ip(xi)
+ elif method == 'cubic' and ndim == 2:
+ ip = CloughTocher2DInterpolator(points, values)
+ return ip(xi)
+ else:
+ raise ValueError("Unknown interpolation method %r for "
+ "%d dimensional data" % (method, ndim))
Added: trunk/scipy/interpolate/tests/test_griddatand.py
===================================================================
--- trunk/scipy/interpolate/tests/test_griddatand.py (rev 0)
+++ trunk/scipy/interpolate/tests/test_griddatand.py 2010-08-31 21:58:55 UTC (rev 6657)
@@ -0,0 +1,57 @@
+import numpy as np
+from numpy.testing import *
+
+from scipy.interpolate.griddatand import griddata
+
+
+class TestGriddata(object):
+
+ def test_alternative_call(self):
+ x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
+ dtype=np.double)
+ y = (np.arange(x.shape[0], dtype=np.double)[:,None]
+ + np.array([0,1])[None,:])
+
+ for method in ('nearest', 'linear', 'cubic'):
+ yi = griddata((x[:,0], x[:,1]), y, (x[:,0], x[:,1]), method=method)
+ assert_almost_equal(y, yi, err_msg=method)
+
+ def test_multivalue_2d(self):
+ x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
+ dtype=np.double)
+ y = (np.arange(x.shape[0], dtype=np.double)[:,None]
+ + np.array([0,1])[None,:])
+
+ for method in ('nearest', 'linear', 'cubic'):
+ yi = griddata(x, y, x, method=method)
+ assert_almost_equal(y, yi, err_msg=method)
+
+ def test_multipoint_2d(self):
+ x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
+ dtype=np.double)
+ y = np.arange(x.shape[0], dtype=np.double)
+
+ xi = x[:,None,:] + np.array([0,0,0])[None,:,None]
+
+ for method in ('nearest', 'linear', 'cubic'):
+ yi = griddata(x, y, xi, method=method)
+
+ assert_equal(yi.shape, (5, 3), err_msg=method)
+ assert_almost_equal(yi, np.tile(y[:,None], (1, 3)), err_msg=method)
+
+ def test_complex_2d(self):
+ x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
+ dtype=np.double)
+ y = np.arange(x.shape[0], dtype=np.double)
+ y = y - 2j*y[::-1]
+
+ xi = x[:,None,:] + np.array([0,0,0])[None,:,None]
+
+ for method in ('nearest', 'linear', 'cubic'):
+ yi = griddata(x, y, xi, method=method)
+
+ assert_equal(yi.shape, (5, 3), err_msg=method)
+ assert_almost_equal(yi, np.tile(y[:,None], (1, 3)), err_msg=method)
+
+if __name__ == "__main__":
+ run_module_suite()
More information about the Scipy-svn
mailing list