Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-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) ndarray of floats 

34 Data point coordinates. 

35 y : (Npoints,) 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 

49 Notes 

50 ----- 

51 Uses ``scipy.spatial.cKDTree`` 

52 

53 Examples 

54 -------- 

55 We can interpolate values on a 2D plane: 

56 

57 >>> from scipy.interpolate import NearestNDInterpolator 

58 >>> import numpy as np 

59 >>> import matplotlib.pyplot as plt 

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

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

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

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

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

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

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

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

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

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

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

71 >>> plt.legend() 

72 >>> plt.colorbar() 

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

74 >>> plt.show() 

75 

76 See also 

77 -------- 

78 griddata : 

79 Interpolate unstructured D-D data. 

80 LinearNDInterpolator : 

81 Piecewise linear interpolant in N dimensions. 

82 CloughTocher2DInterpolator : 

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

84 

85 """ 

86 

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

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

89 need_contiguous=False, 

90 need_values=False) 

91 if tree_options is None: 

92 tree_options = dict() 

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

94 self.values = np.asarray(y) 

95 

96 def __call__(self, *args): 

97 """ 

98 Evaluate interpolator at given points. 

99 

100 Parameters 

101 ---------- 

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

103 Points where to interpolate data at. 

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

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

106 

107 """ 

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

109 xi = self._check_call_shape(xi) 

110 xi = self._scale_x(xi) 

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

112 return self.values[i] 

113 

114 

115#------------------------------------------------------------------------------ 

116# Convenience interface function 

117#------------------------------------------------------------------------------ 

118 

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

120 rescale=False): 

121 """ 

122 Interpolate unstructured D-D data. 

123 

124 Parameters 

125 ---------- 

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

127 Data point coordinates. 

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

129 Data values. 

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

131 Points at which to interpolate data. 

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

133 Method of interpolation. One of 

134 

135 ``nearest`` 

136 return the value at the data point closest to 

137 the point of interpolation. See `NearestNDInterpolator` for 

138 more details. 

139 

140 ``linear`` 

141 tessellate the input point set to N-D 

142 simplices, and interpolate linearly on each simplex. See 

143 `LinearNDInterpolator` for more details. 

144 

145 ``cubic`` (1-D) 

146 return the value determined from a cubic 

147 spline. 

148 

149 ``cubic`` (2-D) 

150 return the value determined from a 

151 piecewise cubic, continuously differentiable (C1), and 

152 approximately curvature-minimizing polynomial surface. See 

153 `CloughTocher2DInterpolator` for more details. 

154 fill_value : float, optional 

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

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

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

158 'nearest' method. 

159 rescale : bool, optional 

160 Rescale points to unit cube before performing interpolation. 

161 This is useful if some of the input dimensions have 

162 incommensurable units and differ by many orders of magnitude. 

163 

164 .. versionadded:: 0.14.0 

165 

166 Returns 

167 ------- 

168 ndarray 

169 Array of interpolated values. 

170 

171 Notes 

172 ----- 

173 

174 .. versionadded:: 0.9 

175 

176 For data on a regular grid use `interpn` instead. 

177 

178 Examples 

179 -------- 

180 

181 Suppose we want to interpolate the 2-D function 

182 

183 >>> import numpy as np 

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

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

186 

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

188 

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

190 

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

192 

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

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

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

196 

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

198 interpolation methods: 

199 

200 >>> from scipy.interpolate import griddata 

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

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

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

204 

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

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

207 cubic interpolant gives the best results: 

208 

209 >>> import matplotlib.pyplot as plt 

210 >>> plt.subplot(221) 

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

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

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

214 >>> plt.subplot(222) 

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

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

217 >>> plt.subplot(223) 

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

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

220 >>> plt.subplot(224) 

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

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

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

224 >>> plt.show() 

225 

226 See Also 

227 -------- 

228 LinearNDInterpolator : 

229 Piecewise linear interpolant in N dimensions. 

230 NearestNDInterpolator : 

231 Nearest-neighbor interpolation in N dimensions. 

232 CloughTocher2DInterpolator : 

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

234 

235 """ 

236 

237 points = _ndim_coords_from_arrays(points) 

238 

239 if points.ndim < 2: 

240 ndim = points.ndim 

241 else: 

242 ndim = points.shape[-1] 

243 

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

245 from ._interpolate import interp1d 

246 points = points.ravel() 

247 if isinstance(xi, tuple): 

248 if len(xi) != 1: 

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

250 xi, = xi 

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

252 idx = np.argsort(points) 

253 points = points[idx] 

254 values = values[idx] 

255 if method == 'nearest': 

256 fill_value = 'extrapolate' 

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

258 fill_value=fill_value) 

259 return ip(xi) 

260 elif method == 'nearest': 

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

262 return ip(xi) 

263 elif method == 'linear': 

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

265 rescale=rescale) 

266 return ip(xi) 

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

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

269 rescale=rescale) 

270 return ip(xi) 

271 else: 

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

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