Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/interpolate/ndgriddata.py : 17%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Convenience interface to N-D interpolation
4.. versionadded:: 0.9
6"""
7import numpy as np
8from .interpnd import LinearNDInterpolator, NDInterpolatorBase, \
9 CloughTocher2DInterpolator, _ndim_coords_from_arrays
10from scipy.spatial import cKDTree
12__all__ = ['griddata', 'NearestNDInterpolator', 'LinearNDInterpolator',
13 'CloughTocher2DInterpolator']
15#------------------------------------------------------------------------------
16# Nearest-neighbor interpolation
17#------------------------------------------------------------------------------
20class NearestNDInterpolator(NDInterpolatorBase):
21 """
22 NearestNDInterpolator(x, y)
24 Nearest-neighbor interpolation in N dimensions.
26 .. versionadded:: 0.9
28 Methods
29 -------
30 __call__
32 Parameters
33 ----------
34 x : (Npoints, Ndims) ndarray of floats
35 Data point coordinates.
36 y : (Npoints,) ndarray of float or complex
37 Data values.
38 rescale : boolean, optional
39 Rescale points to unit cube before performing interpolation.
40 This is useful if some of the input dimensions have
41 incommensurable units and differ by many orders of magnitude.
43 .. versionadded:: 0.14.0
44 tree_options : dict, optional
45 Options passed to the underlying ``cKDTree``.
47 .. versionadded:: 0.17.0
50 Notes
51 -----
52 Uses ``scipy.spatial.cKDTree``
54 """
56 def __init__(self, x, y, rescale=False, tree_options=None):
57 NDInterpolatorBase.__init__(self, x, y, rescale=rescale,
58 need_contiguous=False,
59 need_values=False)
60 if tree_options is None:
61 tree_options = dict()
62 self.tree = cKDTree(self.points, **tree_options)
63 self.values = np.asarray(y)
65 def __call__(self, *args):
66 """
67 Evaluate interpolator at given points.
69 Parameters
70 ----------
71 xi : ndarray of float, shape (..., ndim)
72 Points where to interpolate data at.
74 """
75 xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1])
76 xi = self._check_call_shape(xi)
77 xi = self._scale_x(xi)
78 dist, i = self.tree.query(xi)
79 return self.values[i]
82#------------------------------------------------------------------------------
83# Convenience interface function
84#------------------------------------------------------------------------------
86def griddata(points, values, xi, method='linear', fill_value=np.nan,
87 rescale=False):
88 """
89 Interpolate unstructured D-D data.
91 Parameters
92 ----------
93 points : 2-D ndarray of floats with shape (n, D), or length D tuple of 1-D ndarrays with shape (n,).
94 Data point coordinates.
95 values : ndarray of float or complex, shape (n,)
96 Data values.
97 xi : 2-D ndarray of floats with shape (m, D), or length D tuple of ndarrays broadcastable to the same shape.
98 Points at which to interpolate data.
99 method : {'linear', 'nearest', 'cubic'}, optional
100 Method of interpolation. One of
102 ``nearest``
103 return the value at the data point closest to
104 the point of interpolation. See `NearestNDInterpolator` for
105 more details.
107 ``linear``
108 tessellate the input point set to N-D
109 simplices, and interpolate linearly on each simplex. See
110 `LinearNDInterpolator` for more details.
112 ``cubic`` (1-D)
113 return the value determined from a cubic
114 spline.
116 ``cubic`` (2-D)
117 return the value determined from a
118 piecewise cubic, continuously differentiable (C1), and
119 approximately curvature-minimizing polynomial surface. See
120 `CloughTocher2DInterpolator` for more details.
121 fill_value : float, optional
122 Value used to fill in for requested points outside of the
123 convex hull of the input points. If not provided, then the
124 default is ``nan``. This option has no effect for the
125 'nearest' method.
126 rescale : bool, optional
127 Rescale points to unit cube before performing interpolation.
128 This is useful if some of the input dimensions have
129 incommensurable units and differ by many orders of magnitude.
131 .. versionadded:: 0.14.0
133 Returns
134 -------
135 ndarray
136 Array of interpolated values.
138 Notes
139 -----
141 .. versionadded:: 0.9
143 Examples
144 --------
146 Suppose we want to interpolate the 2-D function
148 >>> def func(x, y):
149 ... return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
151 on a grid in [0, 1]x[0, 1]
153 >>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]
155 but we only know its values at 1000 data points:
157 >>> points = np.random.rand(1000, 2)
158 >>> values = func(points[:,0], points[:,1])
160 This can be done with `griddata` -- below we try out all of the
161 interpolation methods:
163 >>> from scipy.interpolate import griddata
164 >>> grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
165 >>> grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
166 >>> grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')
168 One can see that the exact result is reproduced by all of the
169 methods to some degree, but for this smooth function the piecewise
170 cubic interpolant gives the best results:
172 >>> import matplotlib.pyplot as plt
173 >>> plt.subplot(221)
174 >>> plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
175 >>> plt.plot(points[:,0], points[:,1], 'k.', ms=1)
176 >>> plt.title('Original')
177 >>> plt.subplot(222)
178 >>> plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
179 >>> plt.title('Nearest')
180 >>> plt.subplot(223)
181 >>> plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
182 >>> plt.title('Linear')
183 >>> plt.subplot(224)
184 >>> plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
185 >>> plt.title('Cubic')
186 >>> plt.gcf().set_size_inches(6, 6)
187 >>> plt.show()
189 """
191 points = _ndim_coords_from_arrays(points)
193 if points.ndim < 2:
194 ndim = points.ndim
195 else:
196 ndim = points.shape[-1]
198 if ndim == 1 and method in ('nearest', 'linear', 'cubic'):
199 from .interpolate import interp1d
200 points = points.ravel()
201 if isinstance(xi, tuple):
202 if len(xi) != 1:
203 raise ValueError("invalid number of dimensions in xi")
204 xi, = xi
205 # Sort points/values together, necessary as input for interp1d
206 idx = np.argsort(points)
207 points = points[idx]
208 values = values[idx]
209 if method == 'nearest':
210 fill_value = 'extrapolate'
211 ip = interp1d(points, values, kind=method, axis=0, bounds_error=False,
212 fill_value=fill_value)
213 return ip(xi)
214 elif method == 'nearest':
215 ip = NearestNDInterpolator(points, values, rescale=rescale)
216 return ip(xi)
217 elif method == 'linear':
218 ip = LinearNDInterpolator(points, values, fill_value=fill_value,
219 rescale=rescale)
220 return ip(xi)
221 elif method == 'cubic' and ndim == 2:
222 ip = CloughTocher2DInterpolator(points, values, fill_value=fill_value,
223 rescale=rescale)
224 return ip(xi)
225 else:
226 raise ValueError("Unknown interpolation method %r for "
227 "%d dimensional data" % (method, ndim))