1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

""" 

Convenience interface to N-D interpolation 

 

.. versionadded:: 0.9 

 

""" 

from __future__ import division, print_function, absolute_import 

 

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

 

Nearest-neighbour interpolation in N dimensions. 

 

.. versionadded:: 0.9 

 

Methods 

------- 

__call__ 

 

Parameters 

---------- 

x : (Npoints, Ndims) ndarray of floats 

Data point coordinates. 

y : (Npoints,) ndarray of float or complex 

Data values. 

rescale : boolean, optional 

Rescale points to unit cube before performing interpolation. 

This is useful if some of the input dimensions have 

incommensurable units and differ by many orders of magnitude. 

 

.. versionadded:: 0.14.0 

tree_options : dict, optional 

Options passed to the underlying ``cKDTree``. 

 

.. versionadded:: 0.17.0 

 

 

Notes 

----- 

Uses ``scipy.spatial.cKDTree`` 

 

""" 

 

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

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

need_contiguous=False, 

need_values=False) 

if tree_options is None: 

tree_options = dict() 

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

self.values = y 

 

def __call__(self, *args): 

""" 

Evaluate interpolator at given points. 

 

Parameters 

---------- 

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

Points where to interpolate data at. 

 

""" 

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

xi = self._check_call_shape(xi) 

xi = self._scale_x(xi) 

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

return self.values[i] 

 

 

#------------------------------------------------------------------------------ 

# Convenience interface function 

#------------------------------------------------------------------------------ 

 

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

rescale=False): 

""" 

Interpolate unstructured D-dimensional data. 

 

Parameters 

---------- 

points : ndarray of floats, shape (n, D) 

Data point coordinates. Can either be an array of 

shape (n, D), or a tuple of `ndim` arrays. 

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

Data values. 

xi : 2-D ndarray of float or tuple of 1-D array, shape (M, D) 

Points at which to interpolate data. 

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

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

tessellate 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 determined 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. 

fill_value : float, optional 

Value used to fill in for requested points outside of the 

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

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

'nearest' method. 

rescale : bool, optional 

Rescale points to unit cube before performing interpolation. 

This is useful if some of the input dimensions have 

incommensurable units and differ by many orders of magnitude. 

 

.. versionadded:: 0.14.0 

 

Returns 

------- 

ndarray 

Array of interpolated values. 

 

Notes 

----- 

 

.. versionadded:: 0.9 

 

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) 

 

if points.ndim < 2: 

ndim = points.ndim 

else: 

ndim = points.shape[-1] 

 

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

from .interpolate import interp1d 

points = points.ravel() 

if isinstance(xi, tuple): 

if len(xi) != 1: 

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

xi, = xi 

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

idx = np.argsort(points) 

points = points[idx] 

values = values[idx] 

if method == 'nearest': 

fill_value = 'extrapolate' 

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

fill_value=fill_value) 

return ip(xi) 

elif method == 'nearest': 

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

return ip(xi) 

elif method == 'linear': 

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

rescale=rescale) 

return ip(xi) 

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

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

rescale=rescale) 

return ip(xi) 

else: 

raise ValueError("Unknown interpolation method %r for " 

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