Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_rgi.py: 16%
202 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1__all__ = ['RegularGridInterpolator', 'interpn']
3import itertools
5import numpy as np
7from .interpnd import _ndim_coords_from_arrays
8from ._cubic import PchipInterpolator
9from ._rgi_cython import evaluate_linear_2d, find_indices
10from ._bsplines import make_interp_spline
11from ._fitpack2 import RectBivariateSpline
14def _check_points(points):
15 descending_dimensions = []
16 grid = []
17 for i, p in enumerate(points):
18 # early make points float
19 # see https://github.com/scipy/scipy/pull/17230
20 p = np.asarray(p, dtype=float)
21 if not np.all(p[1:] > p[:-1]):
22 if np.all(p[1:] < p[:-1]):
23 # input is descending, so make it ascending
24 descending_dimensions.append(i)
25 p = np.flip(p)
26 else:
27 raise ValueError(
28 "The points in dimension %d must be strictly "
29 "ascending or descending" % i)
30 # see https://github.com/scipy/scipy/issues/17716
31 p = np.ascontiguousarray(p)
32 grid.append(p)
33 return tuple(grid), tuple(descending_dimensions)
36def _check_dimensionality(points, values):
37 if len(points) > values.ndim:
38 raise ValueError("There are %d point arrays, but values has %d "
39 "dimensions" % (len(points), values.ndim))
40 for i, p in enumerate(points):
41 if not np.asarray(p).ndim == 1:
42 raise ValueError("The points in dimension %d must be "
43 "1-dimensional" % i)
44 if not values.shape[i] == len(p):
45 raise ValueError("There are %d points and %d values in "
46 "dimension %d" % (len(p), values.shape[i], i))
49class RegularGridInterpolator:
50 """
51 Interpolation on a regular or rectilinear grid in arbitrary dimensions.
53 The data must be defined on a rectilinear grid; that is, a rectangular
54 grid with even or uneven spacing. Linear, nearest-neighbor, spline
55 interpolations are supported. After setting up the interpolator object,
56 the interpolation method may be chosen at each evaluation.
58 Parameters
59 ----------
60 points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, )
61 The points defining the regular grid in n dimensions. The points in
62 each dimension (i.e. every elements of the points tuple) must be
63 strictly ascending or descending.
65 values : array_like, shape (m1, ..., mn, ...)
66 The data on the regular grid in n dimensions. Complex data can be
67 acceptable.
69 method : str, optional
70 The method of interpolation to perform. Supported are "linear",
71 "nearest", "slinear", "cubic", "quintic" and "pchip". This
72 parameter will become the default for the object's ``__call__``
73 method. Default is "linear".
75 bounds_error : bool, optional
76 If True, when interpolated values are requested outside of the
77 domain of the input data, a ValueError is raised.
78 If False, then `fill_value` is used.
79 Default is True.
81 fill_value : float or None, optional
82 The value to use for points outside of the interpolation domain.
83 If None, values outside the domain are extrapolated.
84 Default is ``np.nan``.
86 Methods
87 -------
88 __call__
90 Attributes
91 ----------
92 grid : tuple of ndarrays
93 The points defining the regular grid in n dimensions.
94 This tuple defines the full grid via
95 ``np.meshgrid(*grid, indexing='ij')``
96 values : ndarray
97 Data values at the grid.
98 method : str
99 Interpolation method.
100 fill_value : float or ``None``
101 Use this value for out-of-bounds arguments to `__call__`.
102 bounds_error : bool
103 If ``True``, out-of-bounds argument raise a ``ValueError``.
105 Notes
106 -----
107 Contrary to `LinearNDInterpolator` and `NearestNDInterpolator`, this class
108 avoids expensive triangulation of the input data by taking advantage of the
109 regular grid structure.
111 In other words, this class assumes that the data is defined on a
112 *rectilinear* grid.
114 .. versionadded:: 0.14
116 The 'slinear'(k=1), 'cubic'(k=3), and 'quintic'(k=5) methods are
117 tensor-product spline interpolators, where `k` is the spline degree,
118 If any dimension has fewer points than `k` + 1, an error will be raised.
120 .. versionadded:: 1.9
122 If the input data is such that dimensions have incommensurate
123 units and differ by many orders of magnitude, the interpolant may have
124 numerical artifacts. Consider rescaling the data before interpolating.
126 Examples
127 --------
128 **Evaluate a function on the points of a 3-D grid**
130 As a first example, we evaluate a simple example function on the points of
131 a 3-D grid:
133 >>> from scipy.interpolate import RegularGridInterpolator
134 >>> import numpy as np
135 >>> def f(x, y, z):
136 ... return 2 * x**3 + 3 * y**2 - z
137 >>> x = np.linspace(1, 4, 11)
138 >>> y = np.linspace(4, 7, 22)
139 >>> z = np.linspace(7, 9, 33)
140 >>> xg, yg ,zg = np.meshgrid(x, y, z, indexing='ij', sparse=True)
141 >>> data = f(xg, yg, zg)
143 ``data`` is now a 3-D array with ``data[i, j, k] = f(x[i], y[j], z[k])``.
144 Next, define an interpolating function from this data:
146 >>> interp = RegularGridInterpolator((x, y, z), data)
148 Evaluate the interpolating function at the two points
149 ``(x,y,z) = (2.1, 6.2, 8.3)`` and ``(3.3, 5.2, 7.1)``:
151 >>> pts = np.array([[2.1, 6.2, 8.3],
152 ... [3.3, 5.2, 7.1]])
153 >>> interp(pts)
154 array([ 125.80469388, 146.30069388])
156 which is indeed a close approximation to
158 >>> f(2.1, 6.2, 8.3), f(3.3, 5.2, 7.1)
159 (125.54200000000002, 145.894)
161 **Interpolate and extrapolate a 2D dataset**
163 As a second example, we interpolate and extrapolate a 2D data set:
165 >>> x, y = np.array([-2, 0, 4]), np.array([-2, 0, 2, 5])
166 >>> def ff(x, y):
167 ... return x**2 + y**2
169 >>> xg, yg = np.meshgrid(x, y, indexing='ij')
170 >>> data = ff(xg, yg)
171 >>> interp = RegularGridInterpolator((x, y), data,
172 ... bounds_error=False, fill_value=None)
174 >>> import matplotlib.pyplot as plt
175 >>> fig = plt.figure()
176 >>> ax = fig.add_subplot(projection='3d')
177 >>> ax.scatter(xg.ravel(), yg.ravel(), data.ravel(),
178 ... s=60, c='k', label='data')
180 Evaluate and plot the interpolator on a finer grid
182 >>> xx = np.linspace(-4, 9, 31)
183 >>> yy = np.linspace(-4, 9, 31)
184 >>> X, Y = np.meshgrid(xx, yy, indexing='ij')
186 >>> # interpolator
187 >>> ax.plot_wireframe(X, Y, interp((X, Y)), rstride=3, cstride=3,
188 ... alpha=0.4, color='m', label='linear interp')
190 >>> # ground truth
191 >>> ax.plot_wireframe(X, Y, ff(X, Y), rstride=3, cstride=3,
192 ... alpha=0.4, label='ground truth')
193 >>> plt.legend()
194 >>> plt.show()
196 Other examples are given
197 :ref:`in the tutorial <tutorial-interpolate_regular_grid_interpolator>`.
199 See Also
200 --------
201 NearestNDInterpolator : Nearest neighbor interpolation on *unstructured*
202 data in N dimensions
204 LinearNDInterpolator : Piecewise linear interpolant on *unstructured* data
205 in N dimensions
207 interpn : a convenience function which wraps `RegularGridInterpolator`
209 scipy.ndimage.map_coordinates : interpolation on grids with equal spacing
210 (suitable for e.g., N-D image resampling)
212 References
213 ----------
214 .. [1] Python package *regulargrid* by Johannes Buchner, see
215 https://pypi.python.org/pypi/regulargrid/
216 .. [2] Wikipedia, "Trilinear interpolation",
217 https://en.wikipedia.org/wiki/Trilinear_interpolation
218 .. [3] Weiser, Alan, and Sergio E. Zarantonello. "A note on piecewise linear
219 and multilinear table interpolation in many dimensions." MATH.
220 COMPUT. 50.181 (1988): 189-196.
221 https://www.ams.org/journals/mcom/1988-50-181/S0025-5718-1988-0917826-0/S0025-5718-1988-0917826-0.pdf
222 :doi:`10.1090/S0025-5718-1988-0917826-0`
224 """
225 # this class is based on code originally programmed by Johannes Buchner,
226 # see https://github.com/JohannesBuchner/regulargrid
228 _SPLINE_DEGREE_MAP = {"slinear": 1, "cubic": 3, "quintic": 5, 'pchip': 3}
229 _SPLINE_METHODS = list(_SPLINE_DEGREE_MAP.keys())
230 _ALL_METHODS = ["linear", "nearest"] + _SPLINE_METHODS
232 def __init__(self, points, values, method="linear", bounds_error=True,
233 fill_value=np.nan):
234 if method not in self._ALL_METHODS:
235 raise ValueError("Method '%s' is not defined" % method)
236 elif method in self._SPLINE_METHODS:
237 self._validate_grid_dimensions(points, method)
238 self.method = method
239 self.bounds_error = bounds_error
240 self.grid, self._descending_dimensions = _check_points(points)
241 self.values = self._check_values(values)
242 self._check_dimensionality(self.grid, self.values)
243 self.fill_value = self._check_fill_value(self.values, fill_value)
244 if self._descending_dimensions:
245 self.values = np.flip(values, axis=self._descending_dimensions)
247 def _check_dimensionality(self, grid, values):
248 _check_dimensionality(grid, values)
250 def _check_points(self, points):
251 return _check_points(points)
253 def _check_values(self, values):
254 if not hasattr(values, 'ndim'):
255 # allow reasonable duck-typed values
256 values = np.asarray(values)
258 if hasattr(values, 'dtype') and hasattr(values, 'astype'):
259 if not np.issubdtype(values.dtype, np.inexact):
260 values = values.astype(float)
262 return values
264 def _check_fill_value(self, values, fill_value):
265 if fill_value is not None:
266 fill_value_dtype = np.asarray(fill_value).dtype
267 if (hasattr(values, 'dtype') and not
268 np.can_cast(fill_value_dtype, values.dtype,
269 casting='same_kind')):
270 raise ValueError("fill_value must be either 'None' or "
271 "of a type compatible with values")
272 return fill_value
274 def __call__(self, xi, method=None):
275 """
276 Interpolation at coordinates.
278 Parameters
279 ----------
280 xi : ndarray of shape (..., ndim)
281 The coordinates to evaluate the interpolator at.
283 method : str, optional
284 The method of interpolation to perform. Supported are "linear",
285 "nearest", "slinear", "cubic", "quintic" and "pchip". Default is
286 the method chosen when the interpolator was created.
288 Returns
289 -------
290 values_x : ndarray, shape xi.shape[:-1] + values.shape[ndim:]
291 Interpolated values at `xi`. See notes for behaviour when
292 ``xi.ndim == 1``.
294 Notes
295 -----
296 In the case that ``xi.ndim == 1`` a new axis is inserted into
297 the 0 position of the returned array, values_x, so its shape is
298 instead ``(1,) + values.shape[ndim:]``.
300 Examples
301 --------
302 Here we define a nearest-neighbor interpolator of a simple function
304 >>> import numpy as np
305 >>> x, y = np.array([0, 1, 2]), np.array([1, 3, 7])
306 >>> def f(x, y):
307 ... return x**2 + y**2
308 >>> data = f(*np.meshgrid(x, y, indexing='ij', sparse=True))
309 >>> from scipy.interpolate import RegularGridInterpolator
310 >>> interp = RegularGridInterpolator((x, y), data, method='nearest')
312 By construction, the interpolator uses the nearest-neighbor
313 interpolation
315 >>> interp([[1.5, 1.3], [0.3, 4.5]])
316 array([2., 9.])
318 We can however evaluate the linear interpolant by overriding the
319 `method` parameter
321 >>> interp([[1.5, 1.3], [0.3, 4.5]], method='linear')
322 array([ 4.7, 24.3])
323 """
324 is_method_changed = self.method != method
325 method = self.method if method is None else method
326 if method not in self._ALL_METHODS:
327 raise ValueError("Method '%s' is not defined" % method)
329 xi, xi_shape, ndim, nans, out_of_bounds = self._prepare_xi(xi)
331 if method == "linear":
332 indices, norm_distances = self._find_indices(xi.T)
333 if (ndim == 2 and hasattr(self.values, 'dtype') and
334 self.values.ndim == 2 and self.values.flags.writeable and
335 self.values.dtype in (np.float64, np.complex128) and
336 self.values.dtype.byteorder == '='):
337 # until cython supports const fused types, the fast path
338 # cannot support non-writeable values
339 # a fast path
340 out = np.empty(indices.shape[1], dtype=self.values.dtype)
341 result = evaluate_linear_2d(self.values,
342 indices,
343 norm_distances,
344 self.grid,
345 out)
346 else:
347 result = self._evaluate_linear(indices, norm_distances)
348 elif method == "nearest":
349 indices, norm_distances = self._find_indices(xi.T)
350 result = self._evaluate_nearest(indices, norm_distances)
351 elif method in self._SPLINE_METHODS:
352 if is_method_changed:
353 self._validate_grid_dimensions(self.grid, method)
354 result = self._evaluate_spline(xi, method)
356 if not self.bounds_error and self.fill_value is not None:
357 result[out_of_bounds] = self.fill_value
359 # f(nan) = nan, if any
360 if np.any(nans):
361 result[nans] = np.nan
362 return result.reshape(xi_shape[:-1] + self.values.shape[ndim:])
364 def _prepare_xi(self, xi):
365 ndim = len(self.grid)
366 xi = _ndim_coords_from_arrays(xi, ndim=ndim)
367 if xi.shape[-1] != len(self.grid):
368 raise ValueError("The requested sample points xi have dimension "
369 f"{xi.shape[-1]} but this "
370 f"RegularGridInterpolator has dimension {ndim}")
372 xi_shape = xi.shape
373 xi = xi.reshape(-1, xi_shape[-1])
374 xi = np.asarray(xi, dtype=float)
376 # find nans in input
377 nans = np.any(np.isnan(xi), axis=-1)
379 if self.bounds_error:
380 for i, p in enumerate(xi.T):
381 if not np.logical_and(np.all(self.grid[i][0] <= p),
382 np.all(p <= self.grid[i][-1])):
383 raise ValueError("One of the requested xi is out of bounds "
384 "in dimension %d" % i)
385 out_of_bounds = None
386 else:
387 out_of_bounds = self._find_out_of_bounds(xi.T)
389 return xi, xi_shape, ndim, nans, out_of_bounds
391 def _evaluate_linear(self, indices, norm_distances):
392 # slice for broadcasting over trailing dimensions in self.values
393 vslice = (slice(None),) + (None,)*(self.values.ndim - len(indices))
395 # Compute shifting up front before zipping everything together
396 shift_norm_distances = [1 - yi for yi in norm_distances]
397 shift_indices = [i + 1 for i in indices]
399 # The formula for linear interpolation in 2d takes the form:
400 # values = self.values[(i0, i1)] * (1 - y0) * (1 - y1) + \
401 # self.values[(i0, i1 + 1)] * (1 - y0) * y1 + \
402 # self.values[(i0 + 1, i1)] * y0 * (1 - y1) + \
403 # self.values[(i0 + 1, i1 + 1)] * y0 * y1
404 # We pair i with 1 - yi (zipped1) and i + 1 with yi (zipped2)
405 zipped1 = zip(indices, shift_norm_distances)
406 zipped2 = zip(shift_indices, norm_distances)
408 # Take all products of zipped1 and zipped2 and iterate over them
409 # to get the terms in the above formula. This corresponds to iterating
410 # over the vertices of a hypercube.
411 hypercube = itertools.product(*zip(zipped1, zipped2))
412 value = np.array([0.])
413 for h in hypercube:
414 edge_indices, weights = zip(*h)
415 weight = np.array([1.])
416 for w in weights:
417 weight = weight * w
418 term = np.asarray(self.values[edge_indices]) * weight[vslice]
419 value = value + term # cannot use += because broadcasting
420 return value
422 def _evaluate_nearest(self, indices, norm_distances):
423 idx_res = [np.where(yi <= .5, i, i + 1)
424 for i, yi in zip(indices, norm_distances)]
425 return self.values[tuple(idx_res)]
427 def _validate_grid_dimensions(self, points, method):
428 k = self._SPLINE_DEGREE_MAP[method]
429 for i, point in enumerate(points):
430 ndim = len(np.atleast_1d(point))
431 if ndim <= k:
432 raise ValueError(f"There are {ndim} points in dimension {i},"
433 f" but method {method} requires at least "
434 f" {k+1} points per dimension.")
436 def _evaluate_spline(self, xi, method):
437 # ensure xi is 2D list of points to evaluate (`m` is the number of
438 # points and `n` is the number of interpolation dimensions,
439 # ``n == len(self.grid)``.)
440 if xi.ndim == 1:
441 xi = xi.reshape((1, xi.size))
442 m, n = xi.shape
444 # Reorder the axes: n-dimensional process iterates over the
445 # interpolation axes from the last axis downwards: E.g. for a 4D grid
446 # the order of axes is 3, 2, 1, 0. Each 1D interpolation works along
447 # the 0th axis of its argument array (for 1D routine it's its ``y``
448 # array). Thus permute the interpolation axes of `values` *and keep
449 # trailing dimensions trailing*.
450 axes = tuple(range(self.values.ndim))
451 axx = axes[:n][::-1] + axes[n:]
452 values = self.values.transpose(axx)
454 if method == 'pchip':
455 _eval_func = self._do_pchip
456 else:
457 _eval_func = self._do_spline_fit
458 k = self._SPLINE_DEGREE_MAP[method]
460 # Non-stationary procedure: difficult to vectorize this part entirely
461 # into numpy-level operations. Unfortunately this requires explicit
462 # looping over each point in xi.
464 # can at least vectorize the first pass across all points in the
465 # last variable of xi.
466 last_dim = n - 1
467 first_values = _eval_func(self.grid[last_dim],
468 values,
469 xi[:, last_dim],
470 k)
472 # the rest of the dimensions have to be on a per point-in-xi basis
473 shape = (m, *self.values.shape[n:])
474 result = np.empty(shape, dtype=self.values.dtype)
475 for j in range(m):
476 # Main process: Apply 1D interpolate in each dimension
477 # sequentially, starting with the last dimension.
478 # These are then "folded" into the next dimension in-place.
479 folded_values = first_values[j, ...]
480 for i in range(last_dim-1, -1, -1):
481 # Interpolate for each 1D from the last dimensions.
482 # This collapses each 1D sequence into a scalar.
483 folded_values = _eval_func(self.grid[i],
484 folded_values,
485 xi[j, i],
486 k)
487 result[j, ...] = folded_values
489 return result
491 @staticmethod
492 def _do_spline_fit(x, y, pt, k):
493 local_interp = make_interp_spline(x, y, k=k, axis=0)
494 values = local_interp(pt)
495 return values
497 @staticmethod
498 def _do_pchip(x, y, pt, k):
499 local_interp = PchipInterpolator(x, y, axis=0)
500 values = local_interp(pt)
501 return values
503 def _find_indices(self, xi):
504 return find_indices(self.grid, xi)
506 def _find_out_of_bounds(self, xi):
507 # check for out of bounds xi
508 out_of_bounds = np.zeros((xi.shape[1]), dtype=bool)
509 # iterate through dimensions
510 for x, grid in zip(xi, self.grid):
511 out_of_bounds += x < grid[0]
512 out_of_bounds += x > grid[-1]
513 return out_of_bounds
516def interpn(points, values, xi, method="linear", bounds_error=True,
517 fill_value=np.nan):
518 """
519 Multidimensional interpolation on regular or rectilinear grids.
521 Strictly speaking, not all regular grids are supported - this function
522 works on *rectilinear* grids, that is, a rectangular grid with even or
523 uneven spacing.
525 Parameters
526 ----------
527 points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, )
528 The points defining the regular grid in n dimensions. The points in
529 each dimension (i.e. every elements of the points tuple) must be
530 strictly ascending or descending.
532 values : array_like, shape (m1, ..., mn, ...)
533 The data on the regular grid in n dimensions. Complex data can be
534 acceptable.
536 xi : ndarray of shape (..., ndim)
537 The coordinates to sample the gridded data at
539 method : str, optional
540 The method of interpolation to perform. Supported are "linear",
541 "nearest", "slinear", "cubic", "quintic", "pchip", and "splinef2d".
542 "splinef2d" is only supported for 2-dimensional data.
544 bounds_error : bool, optional
545 If True, when interpolated values are requested outside of the
546 domain of the input data, a ValueError is raised.
547 If False, then `fill_value` is used.
549 fill_value : number, optional
550 If provided, the value to use for points outside of the
551 interpolation domain. If None, values outside
552 the domain are extrapolated. Extrapolation is not supported by method
553 "splinef2d".
555 Returns
556 -------
557 values_x : ndarray, shape xi.shape[:-1] + values.shape[ndim:]
558 Interpolated values at `xi`. See notes for behaviour when
559 ``xi.ndim == 1``.
561 Notes
562 -----
564 .. versionadded:: 0.14
566 In the case that ``xi.ndim == 1`` a new axis is inserted into
567 the 0 position of the returned array, values_x, so its shape is
568 instead ``(1,) + values.shape[ndim:]``.
570 If the input data is such that input dimensions have incommensurate
571 units and differ by many orders of magnitude, the interpolant may have
572 numerical artifacts. Consider rescaling the data before interpolation.
574 Examples
575 --------
576 Evaluate a simple example function on the points of a regular 3-D grid:
578 >>> import numpy as np
579 >>> from scipy.interpolate import interpn
580 >>> def value_func_3d(x, y, z):
581 ... return 2 * x + 3 * y - z
582 >>> x = np.linspace(0, 4, 5)
583 >>> y = np.linspace(0, 5, 6)
584 >>> z = np.linspace(0, 6, 7)
585 >>> points = (x, y, z)
586 >>> values = value_func_3d(*np.meshgrid(*points, indexing='ij'))
588 Evaluate the interpolating function at a point
590 >>> point = np.array([2.21, 3.12, 1.15])
591 >>> print(interpn(points, values, point))
592 [12.63]
594 See Also
595 --------
596 NearestNDInterpolator : Nearest neighbor interpolation on unstructured
597 data in N dimensions
599 LinearNDInterpolator : Piecewise linear interpolant on unstructured data
600 in N dimensions
602 RegularGridInterpolator : interpolation on a regular or rectilinear grid
603 in arbitrary dimensions (`interpn` wraps this
604 class).
606 RectBivariateSpline : Bivariate spline approximation over a rectangular mesh
608 scipy.ndimage.map_coordinates : interpolation on grids with equal spacing
609 (suitable for e.g., N-D image resampling)
611 """
612 # sanity check 'method' kwarg
613 if method not in ["linear", "nearest", "cubic", "quintic", "pchip",
614 "splinef2d", "slinear"]:
615 raise ValueError("interpn only understands the methods 'linear', "
616 "'nearest', 'slinear', 'cubic', 'quintic', 'pchip', "
617 f"and 'splinef2d'. You provided {method}.")
619 if not hasattr(values, 'ndim'):
620 values = np.asarray(values)
622 ndim = values.ndim
623 if ndim > 2 and method == "splinef2d":
624 raise ValueError("The method splinef2d can only be used for "
625 "2-dimensional input data")
626 if not bounds_error and fill_value is None and method == "splinef2d":
627 raise ValueError("The method splinef2d does not support extrapolation.")
629 # sanity check consistency of input dimensions
630 if len(points) > ndim:
631 raise ValueError("There are %d point arrays, but values has %d "
632 "dimensions" % (len(points), ndim))
633 if len(points) != ndim and method == 'splinef2d':
634 raise ValueError("The method splinef2d can only be used for "
635 "scalar data with one point per coordinate")
637 grid, descending_dimensions = _check_points(points)
638 _check_dimensionality(grid, values)
640 # sanity check requested xi
641 xi = _ndim_coords_from_arrays(xi, ndim=len(grid))
642 if xi.shape[-1] != len(grid):
643 raise ValueError("The requested sample points xi have dimension "
644 "%d, but this RegularGridInterpolator has "
645 "dimension %d" % (xi.shape[-1], len(grid)))
647 if bounds_error:
648 for i, p in enumerate(xi.T):
649 if not np.logical_and(np.all(grid[i][0] <= p),
650 np.all(p <= grid[i][-1])):
651 raise ValueError("One of the requested xi is out of bounds "
652 "in dimension %d" % i)
654 # perform interpolation
655 if method in ["linear", "nearest", "slinear", "cubic", "quintic", "pchip"]:
656 interp = RegularGridInterpolator(points, values, method=method,
657 bounds_error=bounds_error,
658 fill_value=fill_value)
659 return interp(xi)
660 elif method == "splinef2d":
661 xi_shape = xi.shape
662 xi = xi.reshape(-1, xi.shape[-1])
664 # RectBivariateSpline doesn't support fill_value; we need to wrap here
665 idx_valid = np.all((grid[0][0] <= xi[:, 0], xi[:, 0] <= grid[0][-1],
666 grid[1][0] <= xi[:, 1], xi[:, 1] <= grid[1][-1]),
667 axis=0)
668 result = np.empty_like(xi[:, 0])
670 # make a copy of values for RectBivariateSpline
671 interp = RectBivariateSpline(points[0], points[1], values[:])
672 result[idx_valid] = interp.ev(xi[idx_valid, 0], xi[idx_valid, 1])
673 result[np.logical_not(idx_valid)] = fill_value
675 return result.reshape(xi_shape[:-1])