Hide keyboard shortcuts

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 

3 

4.. versionadded:: 0.9 

5 

6""" 

7import numpy as np 

8from .interpnd import LinearNDInterpolator, NDInterpolatorBase, \ 

9 CloughTocher2DInterpolator, _ndim_coords_from_arrays 

10from scipy.spatial import cKDTree 

11 

12__all__ = ['griddata', 'NearestNDInterpolator', 'LinearNDInterpolator', 

13 'CloughTocher2DInterpolator'] 

14 

15#------------------------------------------------------------------------------ 

16# Nearest-neighbor interpolation 

17#------------------------------------------------------------------------------ 

18 

19 

20class NearestNDInterpolator(NDInterpolatorBase): 

21 """ 

22 NearestNDInterpolator(x, y) 

23 

24 Nearest-neighbor interpolation in N dimensions. 

25 

26 .. versionadded:: 0.9 

27 

28 Methods 

29 ------- 

30 __call__ 

31 

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. 

42 

43 .. versionadded:: 0.14.0 

44 tree_options : dict, optional 

45 Options passed to the underlying ``cKDTree``. 

46 

47 .. versionadded:: 0.17.0 

48 

49 

50 Notes 

51 ----- 

52 Uses ``scipy.spatial.cKDTree`` 

53 

54 """ 

55 

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) 

64 

65 def __call__(self, *args): 

66 """ 

67 Evaluate interpolator at given points. 

68 

69 Parameters 

70 ---------- 

71 xi : ndarray of float, shape (..., ndim) 

72 Points where to interpolate data at. 

73 

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] 

80 

81 

82#------------------------------------------------------------------------------ 

83# Convenience interface function 

84#------------------------------------------------------------------------------ 

85 

86def griddata(points, values, xi, method='linear', fill_value=np.nan, 

87 rescale=False): 

88 """ 

89 Interpolate unstructured D-D data. 

90 

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 

101 

102 ``nearest`` 

103 return the value at the data point closest to 

104 the point of interpolation. See `NearestNDInterpolator` for 

105 more details. 

106 

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. 

111 

112 ``cubic`` (1-D) 

113 return the value determined from a cubic 

114 spline. 

115 

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. 

130 

131 .. versionadded:: 0.14.0 

132 

133 Returns 

134 ------- 

135 ndarray 

136 Array of interpolated values. 

137 

138 Notes 

139 ----- 

140 

141 .. versionadded:: 0.9 

142 

143 Examples 

144 -------- 

145 

146 Suppose we want to interpolate the 2-D function 

147 

148 >>> def func(x, y): 

149 ... return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2 

150 

151 on a grid in [0, 1]x[0, 1] 

152 

153 >>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j] 

154 

155 but we only know its values at 1000 data points: 

156 

157 >>> points = np.random.rand(1000, 2) 

158 >>> values = func(points[:,0], points[:,1]) 

159 

160 This can be done with `griddata` -- below we try out all of the 

161 interpolation methods: 

162 

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') 

167 

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: 

171 

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() 

188 

189 """ 

190 

191 points = _ndim_coords_from_arrays(points) 

192 

193 if points.ndim < 2: 

194 ndim = points.ndim 

195 else: 

196 ndim = points.shape[-1] 

197 

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))