Coverage for /usr/lib/python3/dist-packages/scipy/interpolate/_ndgriddata.py: 17%

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 """NearestNDInterpolator(x, y). 

22 

23 Nearest-neighbor interpolation in N > 1 dimensions. 

24 

25 .. versionadded:: 0.9 

26 

27 Methods 

28 ------- 

29 __call__ 

30 

31 Parameters 

32 ---------- 

33 x : (npoints, ndims) 2-D ndarray of floats 

34 Data point coordinates. 

35 y : (npoints, ) 1-D ndarray of float or complex 

36 Data values. 

37 rescale : boolean, optional 

38 Rescale points to unit cube before performing interpolation. 

39 This is useful if some of the input dimensions have 

40 incommensurable units and differ by many orders of magnitude. 

41 

42 .. versionadded:: 0.14.0 

43 tree_options : dict, optional 

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

45 

46 .. versionadded:: 0.17.0 

47 

48 See Also 

49 -------- 

50 griddata : 

51 Interpolate unstructured D-D data. 

52 LinearNDInterpolator : 

53 Piecewise linear interpolant in N dimensions. 

54 CloughTocher2DInterpolator : 

55 Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D. 

56 interpn : Interpolation on a regular grid or rectilinear grid. 

57 RegularGridInterpolator : Interpolation on a regular or rectilinear grid 

58 in arbitrary dimensions (`interpn` wraps this 

59 class). 

60 

61 Notes 

62 ----- 

63 Uses ``scipy.spatial.cKDTree`` 

64 

65 .. note:: For data on a regular grid use `interpn` instead. 

66 

67 Examples 

68 -------- 

69 We can interpolate values on a 2D plane: 

70 

71 >>> from scipy.interpolate import NearestNDInterpolator 

72 >>> import numpy as np 

73 >>> import matplotlib.pyplot as plt 

74 >>> rng = np.random.default_rng() 

75 >>> x = rng.random(10) - 0.5 

76 >>> y = rng.random(10) - 0.5 

77 >>> z = np.hypot(x, y) 

78 >>> X = np.linspace(min(x), max(x)) 

79 >>> Y = np.linspace(min(y), max(y)) 

80 >>> X, Y = np.meshgrid(X, Y) # 2D grid for interpolation 

81 >>> interp = NearestNDInterpolator(list(zip(x, y)), z) 

82 >>> Z = interp(X, Y) 

83 >>> plt.pcolormesh(X, Y, Z, shading='auto') 

84 >>> plt.plot(x, y, "ok", label="input point") 

85 >>> plt.legend() 

86 >>> plt.colorbar() 

87 >>> plt.axis("equal") 

88 >>> plt.show() 

89 

90 """ 

91 

92 def __init__(self, x, y, rescale=False, tree_options=None): 

93 NDInterpolatorBase.__init__(self, x, y, rescale=rescale, 

94 need_contiguous=False, 

95 need_values=False) 

96 if tree_options is None: 

97 tree_options = dict() 

98 self.tree = cKDTree(self.points, **tree_options) 

99 self.values = np.asarray(y) 

100 

101 def __call__(self, *args): 

102 """ 

103 Evaluate interpolator at given points. 

104 

105 Parameters 

106 ---------- 

107 x1, x2, ... xn : array-like of float 

108 Points where to interpolate data at. 

109 x1, x2, ... xn can be array-like of float with broadcastable shape. 

110 or x1 can be array-like of float with shape ``(..., ndim)`` 

111 

112 """ 

113 xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1]) 

114 xi = self._check_call_shape(xi) 

115 xi = self._scale_x(xi) 

116 dist, i = self.tree.query(xi) 

117 return self.values[i] 

118 

119 

120#------------------------------------------------------------------------------ 

121# Convenience interface function 

122#------------------------------------------------------------------------------ 

123 

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

125 rescale=False): 

126 """ 

127 Interpolate unstructured D-D data. 

128 

129 Parameters 

130 ---------- 

131 points : 2-D ndarray of floats with shape (n, D), or length D tuple of 1-D ndarrays with shape (n,). 

132 Data point coordinates. 

133 values : ndarray of float or complex, shape (n,) 

134 Data values. 

135 xi : 2-D ndarray of floats with shape (m, D), or length D tuple of ndarrays broadcastable to the same shape. 

136 Points at which to interpolate data. 

137 method : {'linear', 'nearest', 'cubic'}, optional 

138 Method of interpolation. One of 

139 

140 ``nearest`` 

141 return the value at the data point closest to 

142 the point of interpolation. See `NearestNDInterpolator` for 

143 more details. 

144 

145 ``linear`` 

146 tessellate the input point set to N-D 

147 simplices, and interpolate linearly on each simplex. See 

148 `LinearNDInterpolator` for more details. 

149 

150 ``cubic`` (1-D) 

151 return the value determined from a cubic 

152 spline. 

153 

154 ``cubic`` (2-D) 

155 return the value determined from a 

156 piecewise cubic, continuously differentiable (C1), and 

157 approximately curvature-minimizing polynomial surface. See 

158 `CloughTocher2DInterpolator` for more details. 

159 fill_value : float, optional 

160 Value used to fill in for requested points outside of the 

161 convex hull of the input points. If not provided, then the 

162 default is ``nan``. This option has no effect for the 

163 'nearest' method. 

164 rescale : bool, optional 

165 Rescale points to unit cube before performing interpolation. 

166 This is useful if some of the input dimensions have 

167 incommensurable units and differ by many orders of magnitude. 

168 

169 .. versionadded:: 0.14.0 

170 

171 Returns 

172 ------- 

173 ndarray 

174 Array of interpolated values. 

175 

176 See Also 

177 -------- 

178 LinearNDInterpolator : 

179 Piecewise linear interpolant in N dimensions. 

180 NearestNDInterpolator : 

181 Nearest-neighbor interpolation in N dimensions. 

182 CloughTocher2DInterpolator : 

183 Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D. 

184 interpn : Interpolation on a regular grid or rectilinear grid. 

185 RegularGridInterpolator : Interpolation on a regular or rectilinear grid 

186 in arbitrary dimensions (`interpn` wraps this 

187 class). 

188 

189 Notes 

190 ----- 

191 

192 .. versionadded:: 0.9 

193 

194 .. note:: For data on a regular grid use `interpn` instead. 

195 

196 Examples 

197 -------- 

198 

199 Suppose we want to interpolate the 2-D function 

200 

201 >>> import numpy as np 

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

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

204 

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

206 

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

208 

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

210 

211 >>> rng = np.random.default_rng() 

212 >>> points = rng.random((1000, 2)) 

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

214 

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

216 interpolation methods: 

217 

218 >>> from scipy.interpolate import griddata 

219 >>> grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest') 

220 >>> grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear') 

221 >>> grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic') 

222 

223 One can see that the exact result is reproduced by all of the 

224 methods to some degree, but for this smooth function the piecewise 

225 cubic interpolant gives the best results: 

226 

227 >>> import matplotlib.pyplot as plt 

228 >>> plt.subplot(221) 

229 >>> plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower') 

230 >>> plt.plot(points[:,0], points[:,1], 'k.', ms=1) 

231 >>> plt.title('Original') 

232 >>> plt.subplot(222) 

233 >>> plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower') 

234 >>> plt.title('Nearest') 

235 >>> plt.subplot(223) 

236 >>> plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower') 

237 >>> plt.title('Linear') 

238 >>> plt.subplot(224) 

239 >>> plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower') 

240 >>> plt.title('Cubic') 

241 >>> plt.gcf().set_size_inches(6, 6) 

242 >>> plt.show() 

243 

244 """ 

245 

246 points = _ndim_coords_from_arrays(points) 

247 

248 if points.ndim < 2: 

249 ndim = points.ndim 

250 else: 

251 ndim = points.shape[-1] 

252 

253 if ndim == 1 and method in ('nearest', 'linear', 'cubic'): 

254 from ._interpolate import interp1d 

255 points = points.ravel() 

256 if isinstance(xi, tuple): 

257 if len(xi) != 1: 

258 raise ValueError("invalid number of dimensions in xi") 

259 xi, = xi 

260 # Sort points/values together, necessary as input for interp1d 

261 idx = np.argsort(points) 

262 points = points[idx] 

263 values = values[idx] 

264 if method == 'nearest': 

265 fill_value = 'extrapolate' 

266 ip = interp1d(points, values, kind=method, axis=0, bounds_error=False, 

267 fill_value=fill_value) 

268 return ip(xi) 

269 elif method == 'nearest': 

270 ip = NearestNDInterpolator(points, values, rescale=rescale) 

271 return ip(xi) 

272 elif method == 'linear': 

273 ip = LinearNDInterpolator(points, values, fill_value=fill_value, 

274 rescale=rescale) 

275 return ip(xi) 

276 elif method == 'cubic' and ndim == 2: 

277 ip = CloughTocher2DInterpolator(points, values, fill_value=fill_value, 

278 rescale=rescale) 

279 return ip(xi) 

280 else: 

281 raise ValueError("Unknown interpolation method %r for " 

282 "%d dimensional data" % (method, ndim))