Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/interpolate/_fitpack2.py: 16%
470 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2fitpack --- curve and surface fitting with splines
4fitpack is based on a collection of Fortran routines DIERCKX
5by P. Dierckx (see http://www.netlib.org/dierckx/) transformed
6to double routines by Pearu Peterson.
7"""
8# Created by Pearu Peterson, June,August 2003
9__all__ = [
10 'UnivariateSpline',
11 'InterpolatedUnivariateSpline',
12 'LSQUnivariateSpline',
13 'BivariateSpline',
14 'LSQBivariateSpline',
15 'SmoothBivariateSpline',
16 'LSQSphereBivariateSpline',
17 'SmoothSphereBivariateSpline',
18 'RectBivariateSpline',
19 'RectSphereBivariateSpline']
22import warnings
24from numpy import zeros, concatenate, ravel, diff, array, ones
25import numpy as np
27from . import _fitpack_impl
28from . import dfitpack
31dfitpack_int = dfitpack.types.intvar.dtype
34# ############### Univariate spline ####################
36_curfit_messages = {1: """
37The required storage space exceeds the available storage space, as
38specified by the parameter nest: nest too small. If nest is already
39large (say nest > m/2), it may also indicate that s is too small.
40The approximation returned is the weighted least-squares spline
41according to the knots t[0],t[1],...,t[n-1]. (n=nest) the parameter fp
42gives the corresponding weighted sum of squared residuals (fp>s).
43""",
44 2: """
45A theoretically impossible result was found during the iteration
46process for finding a smoothing spline with fp = s: s too small.
47There is an approximation returned but the corresponding weighted sum
48of squared residuals does not satisfy the condition abs(fp-s)/s < tol.""",
49 3: """
50The maximal number of iterations maxit (set to 20 by the program)
51allowed for finding a smoothing spline with fp=s has been reached: s
52too small.
53There is an approximation returned but the corresponding weighted sum
54of squared residuals does not satisfy the condition abs(fp-s)/s < tol.""",
55 10: """
56Error on entry, no approximation returned. The following conditions
57must hold:
58xb<=x[0]<x[1]<...<x[m-1]<=xe, w[i]>0, i=0..m-1
59if iopt=-1:
60 xb<t[k+1]<t[k+2]<...<t[n-k-2]<xe"""
61 }
64# UnivariateSpline, ext parameter can be an int or a string
65_extrap_modes = {0: 0, 'extrapolate': 0,
66 1: 1, 'zeros': 1,
67 2: 2, 'raise': 2,
68 3: 3, 'const': 3}
71class UnivariateSpline:
72 """
73 1-D smoothing spline fit to a given set of data points.
75 Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data. `s`
76 specifies the number of knots by specifying a smoothing condition.
78 Parameters
79 ----------
80 x : (N,) array_like
81 1-D array of independent input data. Must be increasing;
82 must be strictly increasing if `s` is 0.
83 y : (N,) array_like
84 1-D array of dependent input data, of the same length as `x`.
85 w : (N,) array_like, optional
86 Weights for spline fitting. Must be positive. If `w` is None,
87 weights are all 1. Default is None.
88 bbox : (2,) array_like, optional
89 2-sequence specifying the boundary of the approximation interval. If
90 `bbox` is None, ``bbox=[x[0], x[-1]]``. Default is None.
91 k : int, optional
92 Degree of the smoothing spline. Must be 1 <= `k` <= 5.
93 ``k = 3`` is a cubic spline. Default is 3.
94 s : float or None, optional
95 Positive smoothing factor used to choose the number of knots. Number
96 of knots will be increased until the smoothing condition is satisfied::
98 sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s
100 However, because of numerical issues, the actual condition is::
102 abs(sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) - s) < 0.001 * s
104 If `s` is None, `s` will be set as `len(w)` for a smoothing spline
105 that uses all data points.
106 If 0, spline will interpolate through all data points. This is
107 equivalent to `InterpolatedUnivariateSpline`.
108 Default is None.
109 The user can use the `s` to control the tradeoff between closeness
110 and smoothness of fit. Larger `s` means more smoothing while smaller
111 values of `s` indicate less smoothing.
112 Recommended values of `s` depend on the weights, `w`. If the weights
113 represent the inverse of the standard-deviation of `y`, then a good
114 `s` value should be found in the range (m-sqrt(2*m),m+sqrt(2*m))
115 where m is the number of datapoints in `x`, `y`, and `w`. This means
116 ``s = len(w)`` should be a good value if ``1/w[i]`` is an
117 estimate of the standard deviation of ``y[i]``.
118 ext : int or str, optional
119 Controls the extrapolation mode for elements
120 not in the interval defined by the knot sequence.
122 * if ext=0 or 'extrapolate', return the extrapolated value.
123 * if ext=1 or 'zeros', return 0
124 * if ext=2 or 'raise', raise a ValueError
125 * if ext=3 of 'const', return the boundary value.
127 Default is 0.
129 check_finite : bool, optional
130 Whether to check that the input arrays contain only finite numbers.
131 Disabling may give a performance gain, but may result in problems
132 (crashes, non-termination or non-sensical results) if the inputs
133 do contain infinities or NaNs.
134 Default is False.
136 See Also
137 --------
138 BivariateSpline :
139 a base class for bivariate splines.
140 SmoothBivariateSpline :
141 a smoothing bivariate spline through the given points
142 LSQBivariateSpline :
143 a bivariate spline using weighted least-squares fitting
144 RectSphereBivariateSpline :
145 a bivariate spline over a rectangular mesh on a sphere
146 SmoothSphereBivariateSpline :
147 a smoothing bivariate spline in spherical coordinates
148 LSQSphereBivariateSpline :
149 a bivariate spline in spherical coordinates using weighted
150 least-squares fitting
151 RectBivariateSpline :
152 a bivariate spline over a rectangular mesh
153 InterpolatedUnivariateSpline :
154 a interpolating univariate spline for a given set of data points.
155 bisplrep :
156 a function to find a bivariate B-spline representation of a surface
157 bisplev :
158 a function to evaluate a bivariate B-spline and its derivatives
159 splrep :
160 a function to find the B-spline representation of a 1-D curve
161 splev :
162 a function to evaluate a B-spline or its derivatives
163 sproot :
164 a function to find the roots of a cubic B-spline
165 splint :
166 a function to evaluate the definite integral of a B-spline between two
167 given points
168 spalde :
169 a function to evaluate all derivatives of a B-spline
171 Notes
172 -----
173 The number of data points must be larger than the spline degree `k`.
175 **NaN handling**: If the input arrays contain ``nan`` values, the result
176 is not useful, since the underlying spline fitting routines cannot deal
177 with ``nan``. A workaround is to use zero weights for not-a-number
178 data points:
180 >>> import numpy as np
181 >>> from scipy.interpolate import UnivariateSpline
182 >>> x, y = np.array([1, 2, 3, 4]), np.array([1, np.nan, 3, 4])
183 >>> w = np.isnan(y)
184 >>> y[w] = 0.
185 >>> spl = UnivariateSpline(x, y, w=~w)
187 Notice the need to replace a ``nan`` by a numerical value (precise value
188 does not matter as long as the corresponding weight is zero.)
190 Examples
191 --------
192 >>> import numpy as np
193 >>> import matplotlib.pyplot as plt
194 >>> from scipy.interpolate import UnivariateSpline
195 >>> rng = np.random.default_rng()
196 >>> x = np.linspace(-3, 3, 50)
197 >>> y = np.exp(-x**2) + 0.1 * rng.standard_normal(50)
198 >>> plt.plot(x, y, 'ro', ms=5)
200 Use the default value for the smoothing parameter:
202 >>> spl = UnivariateSpline(x, y)
203 >>> xs = np.linspace(-3, 3, 1000)
204 >>> plt.plot(xs, spl(xs), 'g', lw=3)
206 Manually change the amount of smoothing:
208 >>> spl.set_smoothing_factor(0.5)
209 >>> plt.plot(xs, spl(xs), 'b', lw=3)
210 >>> plt.show()
212 """
214 def __init__(self, x, y, w=None, bbox=[None]*2, k=3, s=None,
215 ext=0, check_finite=False):
217 x, y, w, bbox, self.ext = self.validate_input(x, y, w, bbox, k, s, ext,
218 check_finite)
220 # _data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
221 data = dfitpack.fpcurf0(x, y, k, w=w, xb=bbox[0],
222 xe=bbox[1], s=s)
223 if data[-1] == 1:
224 # nest too small, setting to maximum bound
225 data = self._reset_nest(data)
226 self._data = data
227 self._reset_class()
229 @staticmethod
230 def validate_input(x, y, w, bbox, k, s, ext, check_finite):
231 x, y, bbox = np.asarray(x), np.asarray(y), np.asarray(bbox)
232 if w is not None:
233 w = np.asarray(w)
234 if check_finite:
235 w_finite = np.isfinite(w).all() if w is not None else True
236 if (not np.isfinite(x).all() or not np.isfinite(y).all() or
237 not w_finite):
238 raise ValueError("x and y array must not contain "
239 "NaNs or infs.")
240 if s is None or s > 0:
241 if not np.all(diff(x) >= 0.0):
242 raise ValueError("x must be increasing if s > 0")
243 else:
244 if not np.all(diff(x) > 0.0):
245 raise ValueError("x must be strictly increasing if s = 0")
246 if x.size != y.size:
247 raise ValueError("x and y should have a same length")
248 elif w is not None and not x.size == y.size == w.size:
249 raise ValueError("x, y, and w should have a same length")
250 elif bbox.shape != (2,):
251 raise ValueError("bbox shape should be (2,)")
252 elif not (1 <= k <= 5):
253 raise ValueError("k should be 1 <= k <= 5")
254 elif s is not None and not s >= 0.0:
255 raise ValueError("s should be s >= 0.0")
257 try:
258 ext = _extrap_modes[ext]
259 except KeyError as e:
260 raise ValueError("Unknown extrapolation mode %s." % ext) from e
262 return x, y, w, bbox, ext
264 @classmethod
265 def _from_tck(cls, tck, ext=0):
266 """Construct a spline object from given tck"""
267 self = cls.__new__(cls)
268 t, c, k = tck
269 self._eval_args = tck
270 # _data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
271 self._data = (None, None, None, None, None, k, None, len(t), t,
272 c, None, None, None, None)
273 self.ext = ext
274 return self
276 def _reset_class(self):
277 data = self._data
278 n, t, c, k, ier = data[7], data[8], data[9], data[5], data[-1]
279 self._eval_args = t[:n], c[:n], k
280 if ier == 0:
281 # the spline returned has a residual sum of squares fp
282 # such that abs(fp-s)/s <= tol with tol a relative
283 # tolerance set to 0.001 by the program
284 pass
285 elif ier == -1:
286 # the spline returned is an interpolating spline
287 self._set_class(InterpolatedUnivariateSpline)
288 elif ier == -2:
289 # the spline returned is the weighted least-squares
290 # polynomial of degree k. In this extreme case fp gives
291 # the upper bound fp0 for the smoothing factor s.
292 self._set_class(LSQUnivariateSpline)
293 else:
294 # error
295 if ier == 1:
296 self._set_class(LSQUnivariateSpline)
297 message = _curfit_messages.get(ier, 'ier=%s' % (ier))
298 warnings.warn(message)
300 def _set_class(self, cls):
301 self._spline_class = cls
302 if self.__class__ in (UnivariateSpline, InterpolatedUnivariateSpline,
303 LSQUnivariateSpline):
304 self.__class__ = cls
305 else:
306 # It's an unknown subclass -- don't change class. cf. #731
307 pass
309 def _reset_nest(self, data, nest=None):
310 n = data[10]
311 if nest is None:
312 k, m = data[5], len(data[0])
313 nest = m+k+1 # this is the maximum bound for nest
314 else:
315 if not n <= nest:
316 raise ValueError("`nest` can only be increased")
317 t, c, fpint, nrdata = [np.resize(data[j], nest) for j in
318 [8, 9, 11, 12]]
320 args = data[:8] + (t, c, n, fpint, nrdata, data[13])
321 data = dfitpack.fpcurf1(*args)
322 return data
324 def set_smoothing_factor(self, s):
325 """ Continue spline computation with the given smoothing
326 factor s and with the knots found at the last call.
328 This routine modifies the spline in place.
330 """
331 data = self._data
332 if data[6] == -1:
333 warnings.warn('smoothing factor unchanged for'
334 'LSQ spline with fixed knots')
335 return
336 args = data[:6] + (s,) + data[7:]
337 data = dfitpack.fpcurf1(*args)
338 if data[-1] == 1:
339 # nest too small, setting to maximum bound
340 data = self._reset_nest(data)
341 self._data = data
342 self._reset_class()
344 def __call__(self, x, nu=0, ext=None):
345 """
346 Evaluate spline (or its nu-th derivative) at positions x.
348 Parameters
349 ----------
350 x : array_like
351 A 1-D array of points at which to return the value of the smoothed
352 spline or its derivatives. Note: `x` can be unordered but the
353 evaluation is more efficient if `x` is (partially) ordered.
354 nu : int
355 The order of derivative of the spline to compute.
356 ext : int
357 Controls the value returned for elements of `x` not in the
358 interval defined by the knot sequence.
360 * if ext=0 or 'extrapolate', return the extrapolated value.
361 * if ext=1 or 'zeros', return 0
362 * if ext=2 or 'raise', raise a ValueError
363 * if ext=3 or 'const', return the boundary value.
365 The default value is 0, passed from the initialization of
366 UnivariateSpline.
368 """
369 x = np.asarray(x)
370 # empty input yields empty output
371 if x.size == 0:
372 return array([])
373 if ext is None:
374 ext = self.ext
375 else:
376 try:
377 ext = _extrap_modes[ext]
378 except KeyError as e:
379 raise ValueError("Unknown extrapolation mode %s." % ext) from e
380 return _fitpack_impl.splev(x, self._eval_args, der=nu, ext=ext)
382 def get_knots(self):
383 """ Return positions of interior knots of the spline.
385 Internally, the knot vector contains ``2*k`` additional boundary knots.
386 """
387 data = self._data
388 k, n = data[5], data[7]
389 return data[8][k:n-k]
391 def get_coeffs(self):
392 """Return spline coefficients."""
393 data = self._data
394 k, n = data[5], data[7]
395 return data[9][:n-k-1]
397 def get_residual(self):
398 """Return weighted sum of squared residuals of the spline approximation.
400 This is equivalent to::
402 sum((w[i] * (y[i]-spl(x[i])))**2, axis=0)
404 """
405 return self._data[10]
407 def integral(self, a, b):
408 """ Return definite integral of the spline between two given points.
410 Parameters
411 ----------
412 a : float
413 Lower limit of integration.
414 b : float
415 Upper limit of integration.
417 Returns
418 -------
419 integral : float
420 The value of the definite integral of the spline between limits.
422 Examples
423 --------
424 >>> import numpy as np
425 >>> from scipy.interpolate import UnivariateSpline
426 >>> x = np.linspace(0, 3, 11)
427 >>> y = x**2
428 >>> spl = UnivariateSpline(x, y)
429 >>> spl.integral(0, 3)
430 9.0
432 which agrees with :math:`\\int x^2 dx = x^3 / 3` between the limits
433 of 0 and 3.
435 A caveat is that this routine assumes the spline to be zero outside of
436 the data limits:
438 >>> spl.integral(-1, 4)
439 9.0
440 >>> spl.integral(-1, 0)
441 0.0
443 """
444 return _fitpack_impl.splint(a, b, self._eval_args)
446 def derivatives(self, x):
447 """ Return all derivatives of the spline at the point x.
449 Parameters
450 ----------
451 x : float
452 The point to evaluate the derivatives at.
454 Returns
455 -------
456 der : ndarray, shape(k+1,)
457 Derivatives of the orders 0 to k.
459 Examples
460 --------
461 >>> import numpy as np
462 >>> from scipy.interpolate import UnivariateSpline
463 >>> x = np.linspace(0, 3, 11)
464 >>> y = x**2
465 >>> spl = UnivariateSpline(x, y)
466 >>> spl.derivatives(1.5)
467 array([2.25, 3.0, 2.0, 0])
469 """
470 return _fitpack_impl.spalde(x, self._eval_args)
472 def roots(self):
473 """ Return the zeros of the spline.
475 Notes
476 -----
477 Restriction: only cubic splines are supported by FITPACK. For non-cubic
478 splines, use `PPoly.root` (see below for an example).
480 Examples
481 --------
483 For some data, this method may miss a root. This happens when one of
484 the spline knots (which FITPACK places automatically) happens to
485 coincide with the true root. A workaround is to convert to `PPoly`,
486 which uses a different root-finding algorithm.
488 For example,
490 >>> x = [1.96, 1.97, 1.98, 1.99, 2.00, 2.01, 2.02, 2.03, 2.04, 2.05]
491 >>> y = [-6.365470e-03, -4.790580e-03, -3.204320e-03, -1.607270e-03,
492 ... 4.440892e-16, 1.616930e-03, 3.243000e-03, 4.877670e-03,
493 ... 6.520430e-03, 8.170770e-03]
494 >>> from scipy.interpolate import UnivariateSpline
495 >>> spl = UnivariateSpline(x, y, s=0)
496 >>> spl.roots()
497 array([], dtype=float64)
499 Converting to a PPoly object does find the roots at `x=2`:
501 >>> from scipy.interpolate import splrep, PPoly
502 >>> tck = splrep(x, y, s=0)
503 >>> ppoly = PPoly.from_spline(tck)
504 >>> ppoly.roots(extrapolate=False)
505 array([2.])
507 See Also
508 --------
509 sproot
510 PPoly.roots
512 """
513 k = self._data[5]
514 if k == 3:
515 return _fitpack_impl.sproot(self._eval_args)
516 raise NotImplementedError('finding roots unsupported for '
517 'non-cubic splines')
519 def derivative(self, n=1):
520 """
521 Construct a new spline representing the derivative of this spline.
523 Parameters
524 ----------
525 n : int, optional
526 Order of derivative to evaluate. Default: 1
528 Returns
529 -------
530 spline : UnivariateSpline
531 Spline of order k2=k-n representing the derivative of this
532 spline.
534 See Also
535 --------
536 splder, antiderivative
538 Notes
539 -----
541 .. versionadded:: 0.13.0
543 Examples
544 --------
545 This can be used for finding maxima of a curve:
547 >>> import numpy as np
548 >>> from scipy.interpolate import UnivariateSpline
549 >>> x = np.linspace(0, 10, 70)
550 >>> y = np.sin(x)
551 >>> spl = UnivariateSpline(x, y, k=4, s=0)
553 Now, differentiate the spline and find the zeros of the
554 derivative. (NB: `sproot` only works for order 3 splines, so we
555 fit an order 4 spline):
557 >>> spl.derivative().roots() / np.pi
558 array([ 0.50000001, 1.5 , 2.49999998])
560 This agrees well with roots :math:`\\pi/2 + n\\pi` of
561 :math:`\\cos(x) = \\sin'(x)`.
563 """
564 tck = _fitpack_impl.splder(self._eval_args, n)
565 # if self.ext is 'const', derivative.ext will be 'zeros'
566 ext = 1 if self.ext == 3 else self.ext
567 return UnivariateSpline._from_tck(tck, ext=ext)
569 def antiderivative(self, n=1):
570 """
571 Construct a new spline representing the antiderivative of this spline.
573 Parameters
574 ----------
575 n : int, optional
576 Order of antiderivative to evaluate. Default: 1
578 Returns
579 -------
580 spline : UnivariateSpline
581 Spline of order k2=k+n representing the antiderivative of this
582 spline.
584 Notes
585 -----
587 .. versionadded:: 0.13.0
589 See Also
590 --------
591 splantider, derivative
593 Examples
594 --------
595 >>> import numpy as np
596 >>> from scipy.interpolate import UnivariateSpline
597 >>> x = np.linspace(0, np.pi/2, 70)
598 >>> y = 1 / np.sqrt(1 - 0.8*np.sin(x)**2)
599 >>> spl = UnivariateSpline(x, y, s=0)
601 The derivative is the inverse operation of the antiderivative,
602 although some floating point error accumulates:
604 >>> spl(1.7), spl.antiderivative().derivative()(1.7)
605 (array(2.1565429877197317), array(2.1565429877201865))
607 Antiderivative can be used to evaluate definite integrals:
609 >>> ispl = spl.antiderivative()
610 >>> ispl(np.pi/2) - ispl(0)
611 2.2572053588768486
613 This is indeed an approximation to the complete elliptic integral
614 :math:`K(m) = \\int_0^{\\pi/2} [1 - m\\sin^2 x]^{-1/2} dx`:
616 >>> from scipy.special import ellipk
617 >>> ellipk(0.8)
618 2.2572053268208538
620 """
621 tck = _fitpack_impl.splantider(self._eval_args, n)
622 return UnivariateSpline._from_tck(tck, self.ext)
625class InterpolatedUnivariateSpline(UnivariateSpline):
626 """
627 1-D interpolating spline for a given set of data points.
629 Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data.
630 Spline function passes through all provided points. Equivalent to
631 `UnivariateSpline` with `s` = 0.
633 Parameters
634 ----------
635 x : (N,) array_like
636 Input dimension of data points -- must be strictly increasing
637 y : (N,) array_like
638 input dimension of data points
639 w : (N,) array_like, optional
640 Weights for spline fitting. Must be positive. If None (default),
641 weights are all 1.
642 bbox : (2,) array_like, optional
643 2-sequence specifying the boundary of the approximation interval. If
644 None (default), ``bbox=[x[0], x[-1]]``.
645 k : int, optional
646 Degree of the smoothing spline. Must be ``1 <= k <= 5``. Default is
647 ``k = 3``, a cubic spline.
648 ext : int or str, optional
649 Controls the extrapolation mode for elements
650 not in the interval defined by the knot sequence.
652 * if ext=0 or 'extrapolate', return the extrapolated value.
653 * if ext=1 or 'zeros', return 0
654 * if ext=2 or 'raise', raise a ValueError
655 * if ext=3 of 'const', return the boundary value.
657 The default value is 0.
659 check_finite : bool, optional
660 Whether to check that the input arrays contain only finite numbers.
661 Disabling may give a performance gain, but may result in problems
662 (crashes, non-termination or non-sensical results) if the inputs
663 do contain infinities or NaNs.
664 Default is False.
666 See Also
667 --------
668 UnivariateSpline :
669 a smooth univariate spline to fit a given set of data points.
670 LSQUnivariateSpline :
671 a spline for which knots are user-selected
672 SmoothBivariateSpline :
673 a smoothing bivariate spline through the given points
674 LSQBivariateSpline :
675 a bivariate spline using weighted least-squares fitting
676 splrep :
677 a function to find the B-spline representation of a 1-D curve
678 splev :
679 a function to evaluate a B-spline or its derivatives
680 sproot :
681 a function to find the roots of a cubic B-spline
682 splint :
683 a function to evaluate the definite integral of a B-spline between two
684 given points
685 spalde :
686 a function to evaluate all derivatives of a B-spline
688 Notes
689 -----
690 The number of data points must be larger than the spline degree `k`.
692 Examples
693 --------
694 >>> import numpy as np
695 >>> import matplotlib.pyplot as plt
696 >>> from scipy.interpolate import InterpolatedUnivariateSpline
697 >>> rng = np.random.default_rng()
698 >>> x = np.linspace(-3, 3, 50)
699 >>> y = np.exp(-x**2) + 0.1 * rng.standard_normal(50)
700 >>> spl = InterpolatedUnivariateSpline(x, y)
701 >>> plt.plot(x, y, 'ro', ms=5)
702 >>> xs = np.linspace(-3, 3, 1000)
703 >>> plt.plot(xs, spl(xs), 'g', lw=3, alpha=0.7)
704 >>> plt.show()
706 Notice that the ``spl(x)`` interpolates `y`:
708 >>> spl.get_residual()
709 0.0
711 """
713 def __init__(self, x, y, w=None, bbox=[None]*2, k=3,
714 ext=0, check_finite=False):
716 x, y, w, bbox, self.ext = self.validate_input(x, y, w, bbox, k, None,
717 ext, check_finite)
718 if not np.all(diff(x) > 0.0):
719 raise ValueError('x must be strictly increasing')
721 # _data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
722 self._data = dfitpack.fpcurf0(x, y, k, w=w, xb=bbox[0],
723 xe=bbox[1], s=0)
724 self._reset_class()
727_fpchec_error_string = """The input parameters have been rejected by fpchec. \
728This means that at least one of the following conditions is violated:
7301) k+1 <= n-k-1 <= m
7312) t(1) <= t(2) <= ... <= t(k+1)
732 t(n-k) <= t(n-k+1) <= ... <= t(n)
7333) t(k+1) < t(k+2) < ... < t(n-k)
7344) t(k+1) <= x(i) <= t(n-k)
7355) The conditions specified by Schoenberg and Whitney must hold
736 for at least one subset of data points, i.e., there must be a
737 subset of data points y(j) such that
738 t(j) < y(j) < t(j+k+1), j=1,2,...,n-k-1
739"""
742class LSQUnivariateSpline(UnivariateSpline):
743 """
744 1-D spline with explicit internal knots.
746 Fits a spline y = spl(x) of degree `k` to the provided `x`, `y` data. `t`
747 specifies the internal knots of the spline
749 Parameters
750 ----------
751 x : (N,) array_like
752 Input dimension of data points -- must be increasing
753 y : (N,) array_like
754 Input dimension of data points
755 t : (M,) array_like
756 interior knots of the spline. Must be in ascending order and::
758 bbox[0] < t[0] < ... < t[-1] < bbox[-1]
760 w : (N,) array_like, optional
761 weights for spline fitting. Must be positive. If None (default),
762 weights are all 1.
763 bbox : (2,) array_like, optional
764 2-sequence specifying the boundary of the approximation interval. If
765 None (default), ``bbox = [x[0], x[-1]]``.
766 k : int, optional
767 Degree of the smoothing spline. Must be 1 <= `k` <= 5.
768 Default is `k` = 3, a cubic spline.
769 ext : int or str, optional
770 Controls the extrapolation mode for elements
771 not in the interval defined by the knot sequence.
773 * if ext=0 or 'extrapolate', return the extrapolated value.
774 * if ext=1 or 'zeros', return 0
775 * if ext=2 or 'raise', raise a ValueError
776 * if ext=3 of 'const', return the boundary value.
778 The default value is 0.
780 check_finite : bool, optional
781 Whether to check that the input arrays contain only finite numbers.
782 Disabling may give a performance gain, but may result in problems
783 (crashes, non-termination or non-sensical results) if the inputs
784 do contain infinities or NaNs.
785 Default is False.
787 Raises
788 ------
789 ValueError
790 If the interior knots do not satisfy the Schoenberg-Whitney conditions
792 See Also
793 --------
794 UnivariateSpline :
795 a smooth univariate spline to fit a given set of data points.
796 InterpolatedUnivariateSpline :
797 a interpolating univariate spline for a given set of data points.
798 splrep :
799 a function to find the B-spline representation of a 1-D curve
800 splev :
801 a function to evaluate a B-spline or its derivatives
802 sproot :
803 a function to find the roots of a cubic B-spline
804 splint :
805 a function to evaluate the definite integral of a B-spline between two
806 given points
807 spalde :
808 a function to evaluate all derivatives of a B-spline
810 Notes
811 -----
812 The number of data points must be larger than the spline degree `k`.
814 Knots `t` must satisfy the Schoenberg-Whitney conditions,
815 i.e., there must be a subset of data points ``x[j]`` such that
816 ``t[j] < x[j] < t[j+k+1]``, for ``j=0, 1,...,n-k-2``.
818 Examples
819 --------
820 >>> import numpy as np
821 >>> from scipy.interpolate import LSQUnivariateSpline, UnivariateSpline
822 >>> import matplotlib.pyplot as plt
823 >>> rng = np.random.default_rng()
824 >>> x = np.linspace(-3, 3, 50)
825 >>> y = np.exp(-x**2) + 0.1 * rng.standard_normal(50)
827 Fit a smoothing spline with a pre-defined internal knots:
829 >>> t = [-1, 0, 1]
830 >>> spl = LSQUnivariateSpline(x, y, t)
832 >>> xs = np.linspace(-3, 3, 1000)
833 >>> plt.plot(x, y, 'ro', ms=5)
834 >>> plt.plot(xs, spl(xs), 'g-', lw=3)
835 >>> plt.show()
837 Check the knot vector:
839 >>> spl.get_knots()
840 array([-3., -1., 0., 1., 3.])
842 Constructing lsq spline using the knots from another spline:
844 >>> x = np.arange(10)
845 >>> s = UnivariateSpline(x, x, s=0)
846 >>> s.get_knots()
847 array([ 0., 2., 3., 4., 5., 6., 7., 9.])
848 >>> knt = s.get_knots()
849 >>> s1 = LSQUnivariateSpline(x, x, knt[1:-1]) # Chop 1st and last knot
850 >>> s1.get_knots()
851 array([ 0., 2., 3., 4., 5., 6., 7., 9.])
853 """
855 def __init__(self, x, y, t, w=None, bbox=[None]*2, k=3,
856 ext=0, check_finite=False):
858 x, y, w, bbox, self.ext = self.validate_input(x, y, w, bbox, k, None,
859 ext, check_finite)
860 if not np.all(diff(x) >= 0.0):
861 raise ValueError('x must be increasing')
863 # _data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
864 xb = bbox[0]
865 xe = bbox[1]
866 if xb is None:
867 xb = x[0]
868 if xe is None:
869 xe = x[-1]
870 t = concatenate(([xb]*(k+1), t, [xe]*(k+1)))
871 n = len(t)
872 if not np.all(t[k+1:n-k]-t[k:n-k-1] > 0, axis=0):
873 raise ValueError('Interior knots t must satisfy '
874 'Schoenberg-Whitney conditions')
875 if not dfitpack.fpchec(x, t, k) == 0:
876 raise ValueError(_fpchec_error_string)
877 data = dfitpack.fpcurfm1(x, y, k, t, w=w, xb=xb, xe=xe)
878 self._data = data[:-3] + (None, None, data[-1])
879 self._reset_class()
882# ############### Bivariate spline ####################
884class _BivariateSplineBase:
885 """ Base class for Bivariate spline s(x,y) interpolation on the rectangle
886 [xb,xe] x [yb, ye] calculated from a given set of data points
887 (x,y,z).
889 See Also
890 --------
891 bisplrep :
892 a function to find a bivariate B-spline representation of a surface
893 bisplev :
894 a function to evaluate a bivariate B-spline and its derivatives
895 BivariateSpline :
896 a base class for bivariate splines.
897 SphereBivariateSpline :
898 a bivariate spline on a spherical grid
899 """
901 @classmethod
902 def _from_tck(cls, tck):
903 """Construct a spline object from given tck and degree"""
904 self = cls.__new__(cls)
905 if len(tck) != 5:
906 raise ValueError("tck should be a 5 element tuple of tx,"
907 " ty, c, kx, ky")
908 self.tck = tck[:3]
909 self.degrees = tck[3:]
910 return self
912 def get_residual(self):
913 """ Return weighted sum of squared residuals of the spline
914 approximation: sum ((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0)
915 """
916 return self.fp
918 def get_knots(self):
919 """ Return a tuple (tx,ty) where tx,ty contain knots positions
920 of the spline with respect to x-, y-variable, respectively.
921 The position of interior and additional knots are given as
922 t[k+1:-k-1] and t[:k+1]=b, t[-k-1:]=e, respectively.
923 """
924 return self.tck[:2]
926 def get_coeffs(self):
927 """ Return spline coefficients."""
928 return self.tck[2]
930 def __call__(self, x, y, dx=0, dy=0, grid=True):
931 """
932 Evaluate the spline or its derivatives at given positions.
934 Parameters
935 ----------
936 x, y : array_like
937 Input coordinates.
939 If `grid` is False, evaluate the spline at points ``(x[i],
940 y[i]), i=0, ..., len(x)-1``. Standard Numpy broadcasting
941 is obeyed.
943 If `grid` is True: evaluate spline at the grid points
944 defined by the coordinate arrays x, y. The arrays must be
945 sorted to increasing order.
947 Note that the axis ordering is inverted relative to
948 the output of meshgrid.
949 dx : int
950 Order of x-derivative
952 .. versionadded:: 0.14.0
953 dy : int
954 Order of y-derivative
956 .. versionadded:: 0.14.0
957 grid : bool
958 Whether to evaluate the results on a grid spanned by the
959 input arrays, or at points specified by the input arrays.
961 .. versionadded:: 0.14.0
963 """
964 x = np.asarray(x)
965 y = np.asarray(y)
967 tx, ty, c = self.tck[:3]
968 kx, ky = self.degrees
969 if grid:
970 if x.size == 0 or y.size == 0:
971 return np.zeros((x.size, y.size), dtype=self.tck[2].dtype)
973 if (x.size >= 2) and (not np.all(np.diff(x) >= 0.0)):
974 raise ValueError("x must be strictly increasing when `grid` is True")
975 if (y.size >= 2) and (not np.all(np.diff(y) >= 0.0)):
976 raise ValueError("y must be strictly increasing when `grid` is True")
978 if dx or dy:
979 z, ier = dfitpack.parder(tx, ty, c, kx, ky, dx, dy, x, y)
980 if not ier == 0:
981 raise ValueError("Error code returned by parder: %s" % ier)
982 else:
983 z, ier = dfitpack.bispev(tx, ty, c, kx, ky, x, y)
984 if not ier == 0:
985 raise ValueError("Error code returned by bispev: %s" % ier)
986 else:
987 # standard Numpy broadcasting
988 if x.shape != y.shape:
989 x, y = np.broadcast_arrays(x, y)
991 shape = x.shape
992 x = x.ravel()
993 y = y.ravel()
995 if x.size == 0 or y.size == 0:
996 return np.zeros(shape, dtype=self.tck[2].dtype)
998 if dx or dy:
999 z, ier = dfitpack.pardeu(tx, ty, c, kx, ky, dx, dy, x, y)
1000 if not ier == 0:
1001 raise ValueError("Error code returned by pardeu: %s" % ier)
1002 else:
1003 z, ier = dfitpack.bispeu(tx, ty, c, kx, ky, x, y)
1004 if not ier == 0:
1005 raise ValueError("Error code returned by bispeu: %s" % ier)
1007 z = z.reshape(shape)
1008 return z
1010 def partial_derivative(self, dx, dy):
1011 """Construct a new spline representing a partial derivative of this
1012 spline.
1014 Parameters
1015 ----------
1016 dx, dy : int
1017 Orders of the derivative in x and y respectively. They must be
1018 non-negative integers and less than the respective degree of the
1019 original spline (self) in that direction (``kx``, ``ky``).
1021 Returns
1022 -------
1023 spline :
1024 A new spline of degrees (``kx - dx``, ``ky - dy``) representing the
1025 derivative of this spline.
1027 Notes
1028 -----
1030 .. versionadded:: 1.9.0
1032 """
1033 if dx == 0 and dy == 0:
1034 return self
1035 else:
1036 kx, ky = self.degrees
1037 if not (dx >= 0 and dy >= 0):
1038 raise ValueError("order of derivative must be positive or"
1039 " zero")
1040 if not (dx < kx and dy < ky):
1041 raise ValueError("order of derivative must be less than"
1042 " degree of spline")
1043 tx, ty, c = self.tck[:3]
1044 newc, ier = dfitpack.pardtc(tx, ty, c, kx, ky, dx, dy)
1045 if ier != 0:
1046 # This should not happen under normal conditions.
1047 raise ValueError("Unexpected error code returned by"
1048 " pardtc: %d" % ier)
1049 nx = len(tx)
1050 ny = len(ty)
1051 newtx = tx[dx:nx - dx]
1052 newty = ty[dy:ny - dy]
1053 newkx, newky = kx - dx, ky - dy
1054 newclen = (nx - dx - kx - 1) * (ny - dy - ky - 1)
1055 return _DerivedBivariateSpline._from_tck((newtx, newty,
1056 newc[:newclen],
1057 newkx, newky))
1060_surfit_messages = {1: """
1061The required storage space exceeds the available storage space: nxest
1062or nyest too small, or s too small.
1063The weighted least-squares spline corresponds to the current set of
1064knots.""",
1065 2: """
1066A theoretically impossible result was found during the iteration
1067process for finding a smoothing spline with fp = s: s too small or
1068badly chosen eps.
1069Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
1070 3: """
1071the maximal number of iterations maxit (set to 20 by the program)
1072allowed for finding a smoothing spline with fp=s has been reached:
1073s too small.
1074Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
1075 4: """
1076No more knots can be added because the number of b-spline coefficients
1077(nx-kx-1)*(ny-ky-1) already exceeds the number of data points m:
1078either s or m too small.
1079The weighted least-squares spline corresponds to the current set of
1080knots.""",
1081 5: """
1082No more knots can be added because the additional knot would (quasi)
1083coincide with an old one: s too small or too large a weight to an
1084inaccurate data point.
1085The weighted least-squares spline corresponds to the current set of
1086knots.""",
1087 10: """
1088Error on entry, no approximation returned. The following conditions
1089must hold:
1090xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
1091If iopt==-1, then
1092 xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
1093 yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""",
1094 -3: """
1095The coefficients of the spline returned have been computed as the
1096minimal norm least-squares solution of a (numerically) rank deficient
1097system (deficiency=%i). If deficiency is large, the results may be
1098inaccurate. Deficiency may strongly depend on the value of eps."""
1099 }
1102class BivariateSpline(_BivariateSplineBase):
1103 """
1104 Base class for bivariate splines.
1106 This describes a spline ``s(x, y)`` of degrees ``kx`` and ``ky`` on
1107 the rectangle ``[xb, xe] * [yb, ye]`` calculated from a given set
1108 of data points ``(x, y, z)``.
1110 This class is meant to be subclassed, not instantiated directly.
1111 To construct these splines, call either `SmoothBivariateSpline` or
1112 `LSQBivariateSpline` or `RectBivariateSpline`.
1114 See Also
1115 --------
1116 UnivariateSpline :
1117 a smooth univariate spline to fit a given set of data points.
1118 SmoothBivariateSpline :
1119 a smoothing bivariate spline through the given points
1120 LSQBivariateSpline :
1121 a bivariate spline using weighted least-squares fitting
1122 RectSphereBivariateSpline :
1123 a bivariate spline over a rectangular mesh on a sphere
1124 SmoothSphereBivariateSpline :
1125 a smoothing bivariate spline in spherical coordinates
1126 LSQSphereBivariateSpline :
1127 a bivariate spline in spherical coordinates using weighted
1128 least-squares fitting
1129 RectBivariateSpline :
1130 a bivariate spline over a rectangular mesh.
1131 bisplrep :
1132 a function to find a bivariate B-spline representation of a surface
1133 bisplev :
1134 a function to evaluate a bivariate B-spline and its derivatives
1135 """
1137 def ev(self, xi, yi, dx=0, dy=0):
1138 """
1139 Evaluate the spline at points
1141 Returns the interpolated value at ``(xi[i], yi[i]),
1142 i=0,...,len(xi)-1``.
1144 Parameters
1145 ----------
1146 xi, yi : array_like
1147 Input coordinates. Standard Numpy broadcasting is obeyed.
1148 dx : int, optional
1149 Order of x-derivative
1151 .. versionadded:: 0.14.0
1152 dy : int, optional
1153 Order of y-derivative
1155 .. versionadded:: 0.14.0
1156 """
1157 return self.__call__(xi, yi, dx=dx, dy=dy, grid=False)
1159 def integral(self, xa, xb, ya, yb):
1160 """
1161 Evaluate the integral of the spline over area [xa,xb] x [ya,yb].
1163 Parameters
1164 ----------
1165 xa, xb : float
1166 The end-points of the x integration interval.
1167 ya, yb : float
1168 The end-points of the y integration interval.
1170 Returns
1171 -------
1172 integ : float
1173 The value of the resulting integral.
1175 """
1176 tx, ty, c = self.tck[:3]
1177 kx, ky = self.degrees
1178 return dfitpack.dblint(tx, ty, c, kx, ky, xa, xb, ya, yb)
1180 @staticmethod
1181 def _validate_input(x, y, z, w, kx, ky, eps):
1182 x, y, z = np.asarray(x), np.asarray(y), np.asarray(z)
1183 if not x.size == y.size == z.size:
1184 raise ValueError('x, y, and z should have a same length')
1186 if w is not None:
1187 w = np.asarray(w)
1188 if x.size != w.size:
1189 raise ValueError('x, y, z, and w should have a same length')
1190 elif not np.all(w >= 0.0):
1191 raise ValueError('w should be positive')
1192 if (eps is not None) and (not 0.0 < eps < 1.0):
1193 raise ValueError('eps should be between (0, 1)')
1194 if not x.size >= (kx + 1) * (ky + 1):
1195 raise ValueError('The length of x, y and z should be at least'
1196 ' (kx+1) * (ky+1)')
1197 return x, y, z, w
1200class _DerivedBivariateSpline(_BivariateSplineBase):
1201 """Bivariate spline constructed from the coefficients and knots of another
1202 spline.
1204 Notes
1205 -----
1206 The class is not meant to be instantiated directly from the data to be
1207 interpolated or smoothed. As a result, its ``fp`` attribute and
1208 ``get_residual`` method are inherited but overriden; ``AttributeError`` is
1209 raised when they are accessed.
1211 The other inherited attributes can be used as usual.
1212 """
1213 _invalid_why = ("is unavailable, because _DerivedBivariateSpline"
1214 " instance is not constructed from data that are to be"
1215 " interpolated or smoothed, but derived from the"
1216 " underlying knots and coefficients of another spline"
1217 " object")
1219 @property
1220 def fp(self):
1221 raise AttributeError("attribute \"fp\" %s" % self._invalid_why)
1223 def get_residual(self):
1224 raise AttributeError("method \"get_residual\" %s" % self._invalid_why)
1227class SmoothBivariateSpline(BivariateSpline):
1228 """
1229 Smooth bivariate spline approximation.
1231 Parameters
1232 ----------
1233 x, y, z : array_like
1234 1-D sequences of data points (order is not important).
1235 w : array_like, optional
1236 Positive 1-D sequence of weights, of same length as `x`, `y` and `z`.
1237 bbox : array_like, optional
1238 Sequence of length 4 specifying the boundary of the rectangular
1239 approximation domain. By default,
1240 ``bbox=[min(x), max(x), min(y), max(y)]``.
1241 kx, ky : ints, optional
1242 Degrees of the bivariate spline. Default is 3.
1243 s : float, optional
1244 Positive smoothing factor defined for estimation condition:
1245 ``sum((w[i]*(z[i]-s(x[i], y[i])))**2, axis=0) <= s``
1246 Default ``s=len(w)`` which should be a good value if ``1/w[i]`` is an
1247 estimate of the standard deviation of ``z[i]``.
1248 eps : float, optional
1249 A threshold for determining the effective rank of an over-determined
1250 linear system of equations. `eps` should have a value within the open
1251 interval ``(0, 1)``, the default is 1e-16.
1253 See Also
1254 --------
1255 BivariateSpline :
1256 a base class for bivariate splines.
1257 UnivariateSpline :
1258 a smooth univariate spline to fit a given set of data points.
1259 LSQBivariateSpline :
1260 a bivariate spline using weighted least-squares fitting
1261 RectSphereBivariateSpline :
1262 a bivariate spline over a rectangular mesh on a sphere
1263 SmoothSphereBivariateSpline :
1264 a smoothing bivariate spline in spherical coordinates
1265 LSQSphereBivariateSpline :
1266 a bivariate spline in spherical coordinates using weighted
1267 least-squares fitting
1268 RectBivariateSpline :
1269 a bivariate spline over a rectangular mesh
1270 bisplrep :
1271 a function to find a bivariate B-spline representation of a surface
1272 bisplev :
1273 a function to evaluate a bivariate B-spline and its derivatives
1275 Notes
1276 -----
1277 The length of `x`, `y` and `z` should be at least ``(kx+1) * (ky+1)``.
1279 If the input data is such that input dimensions have incommensurate
1280 units and differ by many orders of magnitude, the interpolant may have
1281 numerical artifacts. Consider rescaling the data before interpolating.
1283 This routine constructs spline knot vectors automatically via the FITPACK
1284 algorithm. The spline knots may be placed away from the data points. For
1285 some data sets, this routine may fail to construct an interpolating spline,
1286 even if one is requested via ``s=0`` parameter. In such situations, it is
1287 recommended to use `bisplrep` / `bisplev` directly instead of this routine
1288 and, if needed, increase the values of ``nxest`` and ``nyest`` parameters
1289 of `bisplrep`.
1291 For linear interpolation, prefer `LinearNDInterpolator`.
1292 See ``https://gist.github.com/ev-br/8544371b40f414b7eaf3fe6217209bff``
1293 for discussion.
1295 """
1297 def __init__(self, x, y, z, w=None, bbox=[None] * 4, kx=3, ky=3, s=None,
1298 eps=1e-16):
1300 x, y, z, w = self._validate_input(x, y, z, w, kx, ky, eps)
1301 bbox = ravel(bbox)
1302 if not bbox.shape == (4,):
1303 raise ValueError('bbox shape should be (4,)')
1304 if s is not None and not s >= 0.0:
1305 raise ValueError("s should be s >= 0.0")
1307 xb, xe, yb, ye = bbox
1308 nx, tx, ny, ty, c, fp, wrk1, ier = dfitpack.surfit_smth(x, y, z, w,
1309 xb, xe, yb,
1310 ye, kx, ky,
1311 s=s, eps=eps,
1312 lwrk2=1)
1313 if ier > 10: # lwrk2 was to small, re-run
1314 nx, tx, ny, ty, c, fp, wrk1, ier = dfitpack.surfit_smth(x, y, z, w,
1315 xb, xe, yb,
1316 ye, kx, ky,
1317 s=s,
1318 eps=eps,
1319 lwrk2=ier)
1320 if ier in [0, -1, -2]: # normal return
1321 pass
1322 else:
1323 message = _surfit_messages.get(ier, 'ier=%s' % (ier))
1324 warnings.warn(message)
1326 self.fp = fp
1327 self.tck = tx[:nx], ty[:ny], c[:(nx-kx-1)*(ny-ky-1)]
1328 self.degrees = kx, ky
1331class LSQBivariateSpline(BivariateSpline):
1332 """
1333 Weighted least-squares bivariate spline approximation.
1335 Parameters
1336 ----------
1337 x, y, z : array_like
1338 1-D sequences of data points (order is not important).
1339 tx, ty : array_like
1340 Strictly ordered 1-D sequences of knots coordinates.
1341 w : array_like, optional
1342 Positive 1-D array of weights, of the same length as `x`, `y` and `z`.
1343 bbox : (4,) array_like, optional
1344 Sequence of length 4 specifying the boundary of the rectangular
1345 approximation domain. By default,
1346 ``bbox=[min(x,tx),max(x,tx), min(y,ty),max(y,ty)]``.
1347 kx, ky : ints, optional
1348 Degrees of the bivariate spline. Default is 3.
1349 eps : float, optional
1350 A threshold for determining the effective rank of an over-determined
1351 linear system of equations. `eps` should have a value within the open
1352 interval ``(0, 1)``, the default is 1e-16.
1354 See Also
1355 --------
1356 BivariateSpline :
1357 a base class for bivariate splines.
1358 UnivariateSpline :
1359 a smooth univariate spline to fit a given set of data points.
1360 SmoothBivariateSpline :
1361 a smoothing bivariate spline through the given points
1362 RectSphereBivariateSpline :
1363 a bivariate spline over a rectangular mesh on a sphere
1364 SmoothSphereBivariateSpline :
1365 a smoothing bivariate spline in spherical coordinates
1366 LSQSphereBivariateSpline :
1367 a bivariate spline in spherical coordinates using weighted
1368 least-squares fitting
1369 RectBivariateSpline :
1370 a bivariate spline over a rectangular mesh.
1371 bisplrep :
1372 a function to find a bivariate B-spline representation of a surface
1373 bisplev :
1374 a function to evaluate a bivariate B-spline and its derivatives
1376 Notes
1377 -----
1378 The length of `x`, `y` and `z` should be at least ``(kx+1) * (ky+1)``.
1380 If the input data is such that input dimensions have incommensurate
1381 units and differ by many orders of magnitude, the interpolant may have
1382 numerical artifacts. Consider rescaling the data before interpolating.
1384 """
1386 def __init__(self, x, y, z, tx, ty, w=None, bbox=[None]*4, kx=3, ky=3,
1387 eps=None):
1389 x, y, z, w = self._validate_input(x, y, z, w, kx, ky, eps)
1390 bbox = ravel(bbox)
1391 if not bbox.shape == (4,):
1392 raise ValueError('bbox shape should be (4,)')
1394 nx = 2*kx+2+len(tx)
1395 ny = 2*ky+2+len(ty)
1396 # The Fortran subroutine "surfit" (called as dfitpack.surfit_lsq)
1397 # requires that the knot arrays passed as input should be "real
1398 # array(s) of dimension nmax" where "nmax" refers to the greater of nx
1399 # and ny. We pad the tx1/ty1 arrays here so that this is satisfied, and
1400 # slice them to the desired sizes upon return.
1401 nmax = max(nx, ny)
1402 tx1 = zeros((nmax,), float)
1403 ty1 = zeros((nmax,), float)
1404 tx1[kx+1:nx-kx-1] = tx
1405 ty1[ky+1:ny-ky-1] = ty
1407 xb, xe, yb, ye = bbox
1408 tx1, ty1, c, fp, ier = dfitpack.surfit_lsq(x, y, z, nx, tx1, ny, ty1,
1409 w, xb, xe, yb, ye,
1410 kx, ky, eps, lwrk2=1)
1411 if ier > 10:
1412 tx1, ty1, c, fp, ier = dfitpack.surfit_lsq(x, y, z,
1413 nx, tx1, ny, ty1, w,
1414 xb, xe, yb, ye,
1415 kx, ky, eps, lwrk2=ier)
1416 if ier in [0, -1, -2]: # normal return
1417 pass
1418 else:
1419 if ier < -2:
1420 deficiency = (nx-kx-1)*(ny-ky-1)+ier
1421 message = _surfit_messages.get(-3) % (deficiency)
1422 else:
1423 message = _surfit_messages.get(ier, 'ier=%s' % (ier))
1424 warnings.warn(message)
1425 self.fp = fp
1426 self.tck = tx1[:nx], ty1[:ny], c
1427 self.degrees = kx, ky
1430class RectBivariateSpline(BivariateSpline):
1431 """
1432 Bivariate spline approximation over a rectangular mesh.
1434 Can be used for both smoothing and interpolating data.
1436 Parameters
1437 ----------
1438 x,y : array_like
1439 1-D arrays of coordinates in strictly ascending order.
1440 Evaluated points outside the data range will be extrapolated.
1441 z : array_like
1442 2-D array of data with shape (x.size,y.size).
1443 bbox : array_like, optional
1444 Sequence of length 4 specifying the boundary of the rectangular
1445 approximation domain, which means the start and end spline knots of
1446 each dimension are set by these values. By default,
1447 ``bbox=[min(x), max(x), min(y), max(y)]``.
1448 kx, ky : ints, optional
1449 Degrees of the bivariate spline. Default is 3.
1450 s : float, optional
1451 Positive smoothing factor defined for estimation condition:
1452 ``sum((z[i]-f(x[i], y[i]))**2, axis=0) <= s`` where f is a spline
1453 function. Default is ``s=0``, which is for interpolation.
1455 See Also
1456 --------
1457 BivariateSpline :
1458 a base class for bivariate splines.
1459 UnivariateSpline :
1460 a smooth univariate spline to fit a given set of data points.
1461 SmoothBivariateSpline :
1462 a smoothing bivariate spline through the given points
1463 LSQBivariateSpline :
1464 a bivariate spline using weighted least-squares fitting
1465 RectSphereBivariateSpline :
1466 a bivariate spline over a rectangular mesh on a sphere
1467 SmoothSphereBivariateSpline :
1468 a smoothing bivariate spline in spherical coordinates
1469 LSQSphereBivariateSpline :
1470 a bivariate spline in spherical coordinates using weighted
1471 least-squares fitting
1472 bisplrep :
1473 a function to find a bivariate B-spline representation of a surface
1474 bisplev :
1475 a function to evaluate a bivariate B-spline and its derivatives
1477 Notes
1478 -----
1480 If the input data is such that input dimensions have incommensurate
1481 units and differ by many orders of magnitude, the interpolant may have
1482 numerical artifacts. Consider rescaling the data before interpolating.
1484 """
1486 def __init__(self, x, y, z, bbox=[None] * 4, kx=3, ky=3, s=0):
1487 x, y, bbox = ravel(x), ravel(y), ravel(bbox)
1488 z = np.asarray(z)
1489 if not np.all(diff(x) > 0.0):
1490 raise ValueError('x must be strictly increasing')
1491 if not np.all(diff(y) > 0.0):
1492 raise ValueError('y must be strictly increasing')
1493 if not x.size == z.shape[0]:
1494 raise ValueError('x dimension of z must have same number of '
1495 'elements as x')
1496 if not y.size == z.shape[1]:
1497 raise ValueError('y dimension of z must have same number of '
1498 'elements as y')
1499 if not bbox.shape == (4,):
1500 raise ValueError('bbox shape should be (4,)')
1501 if s is not None and not s >= 0.0:
1502 raise ValueError("s should be s >= 0.0")
1504 z = ravel(z)
1505 xb, xe, yb, ye = bbox
1506 nx, tx, ny, ty, c, fp, ier = dfitpack.regrid_smth(x, y, z, xb, xe, yb,
1507 ye, kx, ky, s)
1509 if ier not in [0, -1, -2]:
1510 msg = _surfit_messages.get(ier, 'ier=%s' % (ier))
1511 raise ValueError(msg)
1513 self.fp = fp
1514 self.tck = tx[:nx], ty[:ny], c[:(nx - kx - 1) * (ny - ky - 1)]
1515 self.degrees = kx, ky
1518_spherefit_messages = _surfit_messages.copy()
1519_spherefit_messages[10] = """
1520ERROR. On entry, the input data are controlled on validity. The following
1521 restrictions must be satisfied:
1522 -1<=iopt<=1, m>=2, ntest>=8 ,npest >=8, 0<eps<1,
1523 0<=teta(i)<=pi, 0<=phi(i)<=2*pi, w(i)>0, i=1,...,m
1524 lwrk1 >= 185+52*v+10*u+14*u*v+8*(u-1)*v**2+8*m
1525 kwrk >= m+(ntest-7)*(npest-7)
1526 if iopt=-1: 8<=nt<=ntest , 9<=np<=npest
1527 0<tt(5)<tt(6)<...<tt(nt-4)<pi
1528 0<tp(5)<tp(6)<...<tp(np-4)<2*pi
1529 if iopt>=0: s>=0
1530 if one of these conditions is found to be violated,control
1531 is immediately repassed to the calling program. in that
1532 case there is no approximation returned."""
1533_spherefit_messages[-3] = """
1534WARNING. The coefficients of the spline returned have been computed as the
1535 minimal norm least-squares solution of a (numerically) rank
1536 deficient system (deficiency=%i, rank=%i). Especially if the rank
1537 deficiency, which is computed by 6+(nt-8)*(np-7)+ier, is large,
1538 the results may be inaccurate. They could also seriously depend on
1539 the value of eps."""
1542class SphereBivariateSpline(_BivariateSplineBase):
1543 """
1544 Bivariate spline s(x,y) of degrees 3 on a sphere, calculated from a
1545 given set of data points (theta,phi,r).
1547 .. versionadded:: 0.11.0
1549 See Also
1550 --------
1551 bisplrep :
1552 a function to find a bivariate B-spline representation of a surface
1553 bisplev :
1554 a function to evaluate a bivariate B-spline and its derivatives
1555 UnivariateSpline :
1556 a smooth univariate spline to fit a given set of data points.
1557 SmoothBivariateSpline :
1558 a smoothing bivariate spline through the given points
1559 LSQUnivariateSpline :
1560 a univariate spline using weighted least-squares fitting
1561 """
1563 def __call__(self, theta, phi, dtheta=0, dphi=0, grid=True):
1564 """
1565 Evaluate the spline or its derivatives at given positions.
1567 Parameters
1568 ----------
1569 theta, phi : array_like
1570 Input coordinates.
1572 If `grid` is False, evaluate the spline at points
1573 ``(theta[i], phi[i]), i=0, ..., len(x)-1``. Standard
1574 Numpy broadcasting is obeyed.
1576 If `grid` is True: evaluate spline at the grid points
1577 defined by the coordinate arrays theta, phi. The arrays
1578 must be sorted to increasing order.
1579 dtheta : int, optional
1580 Order of theta-derivative
1582 .. versionadded:: 0.14.0
1583 dphi : int
1584 Order of phi-derivative
1586 .. versionadded:: 0.14.0
1587 grid : bool
1588 Whether to evaluate the results on a grid spanned by the
1589 input arrays, or at points specified by the input arrays.
1591 .. versionadded:: 0.14.0
1593 """
1594 theta = np.asarray(theta)
1595 phi = np.asarray(phi)
1597 if theta.size > 0 and (theta.min() < 0. or theta.max() > np.pi):
1598 raise ValueError("requested theta out of bounds.")
1600 return _BivariateSplineBase.__call__(self, theta, phi,
1601 dx=dtheta, dy=dphi, grid=grid)
1603 def ev(self, theta, phi, dtheta=0, dphi=0):
1604 """
1605 Evaluate the spline at points
1607 Returns the interpolated value at ``(theta[i], phi[i]),
1608 i=0,...,len(theta)-1``.
1610 Parameters
1611 ----------
1612 theta, phi : array_like
1613 Input coordinates. Standard Numpy broadcasting is obeyed.
1614 dtheta : int, optional
1615 Order of theta-derivative
1617 .. versionadded:: 0.14.0
1618 dphi : int, optional
1619 Order of phi-derivative
1621 .. versionadded:: 0.14.0
1622 """
1623 return self.__call__(theta, phi, dtheta=dtheta, dphi=dphi, grid=False)
1626class SmoothSphereBivariateSpline(SphereBivariateSpline):
1627 """
1628 Smooth bivariate spline approximation in spherical coordinates.
1630 .. versionadded:: 0.11.0
1632 Parameters
1633 ----------
1634 theta, phi, r : array_like
1635 1-D sequences of data points (order is not important). Coordinates
1636 must be given in radians. Theta must lie within the interval
1637 ``[0, pi]``, and phi must lie within the interval ``[0, 2pi]``.
1638 w : array_like, optional
1639 Positive 1-D sequence of weights.
1640 s : float, optional
1641 Positive smoothing factor defined for estimation condition:
1642 ``sum((w(i)*(r(i) - s(theta(i), phi(i))))**2, axis=0) <= s``
1643 Default ``s=len(w)`` which should be a good value if ``1/w[i]`` is an
1644 estimate of the standard deviation of ``r[i]``.
1645 eps : float, optional
1646 A threshold for determining the effective rank of an over-determined
1647 linear system of equations. `eps` should have a value within the open
1648 interval ``(0, 1)``, the default is 1e-16.
1650 See Also
1651 --------
1652 BivariateSpline :
1653 a base class for bivariate splines.
1654 UnivariateSpline :
1655 a smooth univariate spline to fit a given set of data points.
1656 SmoothBivariateSpline :
1657 a smoothing bivariate spline through the given points
1658 LSQBivariateSpline :
1659 a bivariate spline using weighted least-squares fitting
1660 RectSphereBivariateSpline :
1661 a bivariate spline over a rectangular mesh on a sphere
1662 LSQSphereBivariateSpline :
1663 a bivariate spline in spherical coordinates using weighted
1664 least-squares fitting
1665 RectBivariateSpline :
1666 a bivariate spline over a rectangular mesh.
1667 bisplrep :
1668 a function to find a bivariate B-spline representation of a surface
1669 bisplev :
1670 a function to evaluate a bivariate B-spline and its derivatives
1672 Notes
1673 -----
1674 For more information, see the FITPACK_ site about this function.
1676 .. _FITPACK: http://www.netlib.org/dierckx/sphere.f
1678 Examples
1679 --------
1680 Suppose we have global data on a coarse grid (the input data does not
1681 have to be on a grid):
1683 >>> import numpy as np
1684 >>> theta = np.linspace(0., np.pi, 7)
1685 >>> phi = np.linspace(0., 2*np.pi, 9)
1686 >>> data = np.empty((theta.shape[0], phi.shape[0]))
1687 >>> data[:,0], data[0,:], data[-1,:] = 0., 0., 0.
1688 >>> data[1:-1,1], data[1:-1,-1] = 1., 1.
1689 >>> data[1,1:-1], data[-2,1:-1] = 1., 1.
1690 >>> data[2:-2,2], data[2:-2,-2] = 2., 2.
1691 >>> data[2,2:-2], data[-3,2:-2] = 2., 2.
1692 >>> data[3,3:-2] = 3.
1693 >>> data = np.roll(data, 4, 1)
1695 We need to set up the interpolator object
1697 >>> lats, lons = np.meshgrid(theta, phi)
1698 >>> from scipy.interpolate import SmoothSphereBivariateSpline
1699 >>> lut = SmoothSphereBivariateSpline(lats.ravel(), lons.ravel(),
1700 ... data.T.ravel(), s=3.5)
1702 As a first test, we'll see what the algorithm returns when run on the
1703 input coordinates
1705 >>> data_orig = lut(theta, phi)
1707 Finally we interpolate the data to a finer grid
1709 >>> fine_lats = np.linspace(0., np.pi, 70)
1710 >>> fine_lons = np.linspace(0., 2 * np.pi, 90)
1712 >>> data_smth = lut(fine_lats, fine_lons)
1714 >>> import matplotlib.pyplot as plt
1715 >>> fig = plt.figure()
1716 >>> ax1 = fig.add_subplot(131)
1717 >>> ax1.imshow(data, interpolation='nearest')
1718 >>> ax2 = fig.add_subplot(132)
1719 >>> ax2.imshow(data_orig, interpolation='nearest')
1720 >>> ax3 = fig.add_subplot(133)
1721 >>> ax3.imshow(data_smth, interpolation='nearest')
1722 >>> plt.show()
1724 """
1726 def __init__(self, theta, phi, r, w=None, s=0., eps=1E-16):
1728 theta, phi, r = np.asarray(theta), np.asarray(phi), np.asarray(r)
1730 # input validation
1731 if not ((0.0 <= theta).all() and (theta <= np.pi).all()):
1732 raise ValueError('theta should be between [0, pi]')
1733 if not ((0.0 <= phi).all() and (phi <= 2.0 * np.pi).all()):
1734 raise ValueError('phi should be between [0, 2pi]')
1735 if w is not None:
1736 w = np.asarray(w)
1737 if not (w >= 0.0).all():
1738 raise ValueError('w should be positive')
1739 if not s >= 0.0:
1740 raise ValueError('s should be positive')
1741 if not 0.0 < eps < 1.0:
1742 raise ValueError('eps should be between (0, 1)')
1744 if np.issubclass_(w, float):
1745 w = ones(len(theta)) * w
1746 nt_, tt_, np_, tp_, c, fp, ier = dfitpack.spherfit_smth(theta, phi,
1747 r, w=w, s=s,
1748 eps=eps)
1749 if ier not in [0, -1, -2]:
1750 message = _spherefit_messages.get(ier, 'ier=%s' % (ier))
1751 raise ValueError(message)
1753 self.fp = fp
1754 self.tck = tt_[:nt_], tp_[:np_], c[:(nt_ - 4) * (np_ - 4)]
1755 self.degrees = (3, 3)
1757 def __call__(self, theta, phi, dtheta=0, dphi=0, grid=True):
1759 theta = np.asarray(theta)
1760 phi = np.asarray(phi)
1762 if phi.size > 0 and (phi.min() < 0. or phi.max() > 2. * np.pi):
1763 raise ValueError("requested phi out of bounds.")
1765 return SphereBivariateSpline.__call__(self, theta, phi, dtheta=dtheta,
1766 dphi=dphi, grid=grid)
1769class LSQSphereBivariateSpline(SphereBivariateSpline):
1770 """
1771 Weighted least-squares bivariate spline approximation in spherical
1772 coordinates.
1774 Determines a smoothing bicubic spline according to a given
1775 set of knots in the `theta` and `phi` directions.
1777 .. versionadded:: 0.11.0
1779 Parameters
1780 ----------
1781 theta, phi, r : array_like
1782 1-D sequences of data points (order is not important). Coordinates
1783 must be given in radians. Theta must lie within the interval
1784 ``[0, pi]``, and phi must lie within the interval ``[0, 2pi]``.
1785 tt, tp : array_like
1786 Strictly ordered 1-D sequences of knots coordinates.
1787 Coordinates must satisfy ``0 < tt[i] < pi``, ``0 < tp[i] < 2*pi``.
1788 w : array_like, optional
1789 Positive 1-D sequence of weights, of the same length as `theta`, `phi`
1790 and `r`.
1791 eps : float, optional
1792 A threshold for determining the effective rank of an over-determined
1793 linear system of equations. `eps` should have a value within the
1794 open interval ``(0, 1)``, the default is 1e-16.
1796 See Also
1797 --------
1798 BivariateSpline :
1799 a base class for bivariate splines.
1800 UnivariateSpline :
1801 a smooth univariate spline to fit a given set of data points.
1802 SmoothBivariateSpline :
1803 a smoothing bivariate spline through the given points
1804 LSQBivariateSpline :
1805 a bivariate spline using weighted least-squares fitting
1806 RectSphereBivariateSpline :
1807 a bivariate spline over a rectangular mesh on a sphere
1808 SmoothSphereBivariateSpline :
1809 a smoothing bivariate spline in spherical coordinates
1810 RectBivariateSpline :
1811 a bivariate spline over a rectangular mesh.
1812 bisplrep :
1813 a function to find a bivariate B-spline representation of a surface
1814 bisplev :
1815 a function to evaluate a bivariate B-spline and its derivatives
1817 Notes
1818 -----
1819 For more information, see the FITPACK_ site about this function.
1821 .. _FITPACK: http://www.netlib.org/dierckx/sphere.f
1823 Examples
1824 --------
1825 Suppose we have global data on a coarse grid (the input data does not
1826 have to be on a grid):
1828 >>> from scipy.interpolate import LSQSphereBivariateSpline
1829 >>> import numpy as np
1830 >>> import matplotlib.pyplot as plt
1832 >>> theta = np.linspace(0, np.pi, num=7)
1833 >>> phi = np.linspace(0, 2*np.pi, num=9)
1834 >>> data = np.empty((theta.shape[0], phi.shape[0]))
1835 >>> data[:,0], data[0,:], data[-1,:] = 0., 0., 0.
1836 >>> data[1:-1,1], data[1:-1,-1] = 1., 1.
1837 >>> data[1,1:-1], data[-2,1:-1] = 1., 1.
1838 >>> data[2:-2,2], data[2:-2,-2] = 2., 2.
1839 >>> data[2,2:-2], data[-3,2:-2] = 2., 2.
1840 >>> data[3,3:-2] = 3.
1841 >>> data = np.roll(data, 4, 1)
1843 We need to set up the interpolator object. Here, we must also specify the
1844 coordinates of the knots to use.
1846 >>> lats, lons = np.meshgrid(theta, phi)
1847 >>> knotst, knotsp = theta.copy(), phi.copy()
1848 >>> knotst[0] += .0001
1849 >>> knotst[-1] -= .0001
1850 >>> knotsp[0] += .0001
1851 >>> knotsp[-1] -= .0001
1852 >>> lut = LSQSphereBivariateSpline(lats.ravel(), lons.ravel(),
1853 ... data.T.ravel(), knotst, knotsp)
1855 As a first test, we'll see what the algorithm returns when run on the
1856 input coordinates
1858 >>> data_orig = lut(theta, phi)
1860 Finally we interpolate the data to a finer grid
1862 >>> fine_lats = np.linspace(0., np.pi, 70)
1863 >>> fine_lons = np.linspace(0., 2*np.pi, 90)
1864 >>> data_lsq = lut(fine_lats, fine_lons)
1866 >>> fig = plt.figure()
1867 >>> ax1 = fig.add_subplot(131)
1868 >>> ax1.imshow(data, interpolation='nearest')
1869 >>> ax2 = fig.add_subplot(132)
1870 >>> ax2.imshow(data_orig, interpolation='nearest')
1871 >>> ax3 = fig.add_subplot(133)
1872 >>> ax3.imshow(data_lsq, interpolation='nearest')
1873 >>> plt.show()
1875 """
1877 def __init__(self, theta, phi, r, tt, tp, w=None, eps=1E-16):
1879 theta, phi, r = np.asarray(theta), np.asarray(phi), np.asarray(r)
1880 tt, tp = np.asarray(tt), np.asarray(tp)
1882 if not ((0.0 <= theta).all() and (theta <= np.pi).all()):
1883 raise ValueError('theta should be between [0, pi]')
1884 if not ((0.0 <= phi).all() and (phi <= 2*np.pi).all()):
1885 raise ValueError('phi should be between [0, 2pi]')
1886 if not ((0.0 < tt).all() and (tt < np.pi).all()):
1887 raise ValueError('tt should be between (0, pi)')
1888 if not ((0.0 < tp).all() and (tp < 2*np.pi).all()):
1889 raise ValueError('tp should be between (0, 2pi)')
1890 if w is not None:
1891 w = np.asarray(w)
1892 if not (w >= 0.0).all():
1893 raise ValueError('w should be positive')
1894 if not 0.0 < eps < 1.0:
1895 raise ValueError('eps should be between (0, 1)')
1897 if np.issubclass_(w, float):
1898 w = ones(len(theta)) * w
1899 nt_, np_ = 8 + len(tt), 8 + len(tp)
1900 tt_, tp_ = zeros((nt_,), float), zeros((np_,), float)
1901 tt_[4:-4], tp_[4:-4] = tt, tp
1902 tt_[-4:], tp_[-4:] = np.pi, 2. * np.pi
1903 tt_, tp_, c, fp, ier = dfitpack.spherfit_lsq(theta, phi, r, tt_, tp_,
1904 w=w, eps=eps)
1905 if ier > 0:
1906 message = _spherefit_messages.get(ier, 'ier=%s' % (ier))
1907 raise ValueError(message)
1909 self.fp = fp
1910 self.tck = tt_, tp_, c
1911 self.degrees = (3, 3)
1913 def __call__(self, theta, phi, dtheta=0, dphi=0, grid=True):
1915 theta = np.asarray(theta)
1916 phi = np.asarray(phi)
1918 if phi.size > 0 and (phi.min() < 0. or phi.max() > 2. * np.pi):
1919 raise ValueError("requested phi out of bounds.")
1921 return SphereBivariateSpline.__call__(self, theta, phi, dtheta=dtheta,
1922 dphi=dphi, grid=grid)
1925_spfit_messages = _surfit_messages.copy()
1926_spfit_messages[10] = """
1927ERROR: on entry, the input data are controlled on validity
1928 the following restrictions must be satisfied.
1929 -1<=iopt(1)<=1, 0<=iopt(2)<=1, 0<=iopt(3)<=1,
1930 -1<=ider(1)<=1, 0<=ider(2)<=1, ider(2)=0 if iopt(2)=0.
1931 -1<=ider(3)<=1, 0<=ider(4)<=1, ider(4)=0 if iopt(3)=0.
1932 mu >= mumin (see above), mv >= 4, nuest >=8, nvest >= 8,
1933 kwrk>=5+mu+mv+nuest+nvest,
1934 lwrk >= 12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+max(nuest,mv+nvest)
1935 0< u(i-1)<u(i)< pi,i=2,..,mu,
1936 -pi<=v(1)< pi, v(1)<v(i-1)<v(i)<v(1)+2*pi, i=3,...,mv
1937 if iopt(1)=-1: 8<=nu<=min(nuest,mu+6+iopt(2)+iopt(3))
1938 0<tu(5)<tu(6)<...<tu(nu-4)< pi
1939 8<=nv<=min(nvest,mv+7)
1940 v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(1)+2*pi
1941 the schoenberg-whitney conditions, i.e. there must be
1942 subset of grid co-ordinates uu(p) and vv(q) such that
1943 tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4
1944 (iopt(2)=1 and iopt(3)=1 also count for a uu-value
1945 tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4
1946 (vv(q) is either a value v(j) or v(j)+2*pi)
1947 if iopt(1)>=0: s>=0
1948 if s=0: nuest>=mu+6+iopt(2)+iopt(3), nvest>=mv+7
1949 if one of these conditions is found to be violated,control is
1950 immediately repassed to the calling program. in that case there is no
1951 approximation returned."""
1954class RectSphereBivariateSpline(SphereBivariateSpline):
1955 """
1956 Bivariate spline approximation over a rectangular mesh on a sphere.
1958 Can be used for smoothing data.
1960 .. versionadded:: 0.11.0
1962 Parameters
1963 ----------
1964 u : array_like
1965 1-D array of colatitude coordinates in strictly ascending order.
1966 Coordinates must be given in radians and lie within the open interval
1967 ``(0, pi)``.
1968 v : array_like
1969 1-D array of longitude coordinates in strictly ascending order.
1970 Coordinates must be given in radians. First element (``v[0]``) must lie
1971 within the interval ``[-pi, pi)``. Last element (``v[-1]``) must satisfy
1972 ``v[-1] <= v[0] + 2*pi``.
1973 r : array_like
1974 2-D array of data with shape ``(u.size, v.size)``.
1975 s : float, optional
1976 Positive smoothing factor defined for estimation condition
1977 (``s=0`` is for interpolation).
1978 pole_continuity : bool or (bool, bool), optional
1979 Order of continuity at the poles ``u=0`` (``pole_continuity[0]``) and
1980 ``u=pi`` (``pole_continuity[1]``). The order of continuity at the pole
1981 will be 1 or 0 when this is True or False, respectively.
1982 Defaults to False.
1983 pole_values : float or (float, float), optional
1984 Data values at the poles ``u=0`` and ``u=pi``. Either the whole
1985 parameter or each individual element can be None. Defaults to None.
1986 pole_exact : bool or (bool, bool), optional
1987 Data value exactness at the poles ``u=0`` and ``u=pi``. If True, the
1988 value is considered to be the right function value, and it will be
1989 fitted exactly. If False, the value will be considered to be a data
1990 value just like the other data values. Defaults to False.
1991 pole_flat : bool or (bool, bool), optional
1992 For the poles at ``u=0`` and ``u=pi``, specify whether or not the
1993 approximation has vanishing derivatives. Defaults to False.
1995 See Also
1996 --------
1997 BivariateSpline :
1998 a base class for bivariate splines.
1999 UnivariateSpline :
2000 a smooth univariate spline to fit a given set of data points.
2001 SmoothBivariateSpline :
2002 a smoothing bivariate spline through the given points
2003 LSQBivariateSpline :
2004 a bivariate spline using weighted least-squares fitting
2005 SmoothSphereBivariateSpline :
2006 a smoothing bivariate spline in spherical coordinates
2007 LSQSphereBivariateSpline :
2008 a bivariate spline in spherical coordinates using weighted
2009 least-squares fitting
2010 RectBivariateSpline :
2011 a bivariate spline over a rectangular mesh.
2012 bisplrep :
2013 a function to find a bivariate B-spline representation of a surface
2014 bisplev :
2015 a function to evaluate a bivariate B-spline and its derivatives
2017 Notes
2018 -----
2019 Currently, only the smoothing spline approximation (``iopt[0] = 0`` and
2020 ``iopt[0] = 1`` in the FITPACK routine) is supported. The exact
2021 least-squares spline approximation is not implemented yet.
2023 When actually performing the interpolation, the requested `v` values must
2024 lie within the same length 2pi interval that the original `v` values were
2025 chosen from.
2027 For more information, see the FITPACK_ site about this function.
2029 .. _FITPACK: http://www.netlib.org/dierckx/spgrid.f
2031 Examples
2032 --------
2033 Suppose we have global data on a coarse grid
2035 >>> import numpy as np
2036 >>> lats = np.linspace(10, 170, 9) * np.pi / 180.
2037 >>> lons = np.linspace(0, 350, 18) * np.pi / 180.
2038 >>> data = np.dot(np.atleast_2d(90. - np.linspace(-80., 80., 18)).T,
2039 ... np.atleast_2d(180. - np.abs(np.linspace(0., 350., 9)))).T
2041 We want to interpolate it to a global one-degree grid
2043 >>> new_lats = np.linspace(1, 180, 180) * np.pi / 180
2044 >>> new_lons = np.linspace(1, 360, 360) * np.pi / 180
2045 >>> new_lats, new_lons = np.meshgrid(new_lats, new_lons)
2047 We need to set up the interpolator object
2049 >>> from scipy.interpolate import RectSphereBivariateSpline
2050 >>> lut = RectSphereBivariateSpline(lats, lons, data)
2052 Finally we interpolate the data. The `RectSphereBivariateSpline` object
2053 only takes 1-D arrays as input, therefore we need to do some reshaping.
2055 >>> data_interp = lut.ev(new_lats.ravel(),
2056 ... new_lons.ravel()).reshape((360, 180)).T
2058 Looking at the original and the interpolated data, one can see that the
2059 interpolant reproduces the original data very well:
2061 >>> import matplotlib.pyplot as plt
2062 >>> fig = plt.figure()
2063 >>> ax1 = fig.add_subplot(211)
2064 >>> ax1.imshow(data, interpolation='nearest')
2065 >>> ax2 = fig.add_subplot(212)
2066 >>> ax2.imshow(data_interp, interpolation='nearest')
2067 >>> plt.show()
2069 Choosing the optimal value of ``s`` can be a delicate task. Recommended
2070 values for ``s`` depend on the accuracy of the data values. If the user
2071 has an idea of the statistical errors on the data, she can also find a
2072 proper estimate for ``s``. By assuming that, if she specifies the
2073 right ``s``, the interpolator will use a spline ``f(u,v)`` which exactly
2074 reproduces the function underlying the data, she can evaluate
2075 ``sum((r(i,j)-s(u(i),v(j)))**2)`` to find a good estimate for this ``s``.
2076 For example, if she knows that the statistical errors on her
2077 ``r(i,j)``-values are not greater than 0.1, she may expect that a good
2078 ``s`` should have a value not larger than ``u.size * v.size * (0.1)**2``.
2080 If nothing is known about the statistical error in ``r(i,j)``, ``s`` must
2081 be determined by trial and error. The best is then to start with a very
2082 large value of ``s`` (to determine the least-squares polynomial and the
2083 corresponding upper bound ``fp0`` for ``s``) and then to progressively
2084 decrease the value of ``s`` (say by a factor 10 in the beginning, i.e.
2085 ``s = fp0 / 10, fp0 / 100, ...`` and more carefully as the approximation
2086 shows more detail) to obtain closer fits.
2088 The interpolation results for different values of ``s`` give some insight
2089 into this process:
2091 >>> fig2 = plt.figure()
2092 >>> s = [3e9, 2e9, 1e9, 1e8]
2093 >>> for idx, sval in enumerate(s, 1):
2094 ... lut = RectSphereBivariateSpline(lats, lons, data, s=sval)
2095 ... data_interp = lut.ev(new_lats.ravel(),
2096 ... new_lons.ravel()).reshape((360, 180)).T
2097 ... ax = fig2.add_subplot(2, 2, idx)
2098 ... ax.imshow(data_interp, interpolation='nearest')
2099 ... ax.set_title(f"s = {sval:g}")
2100 >>> plt.show()
2102 """
2104 def __init__(self, u, v, r, s=0., pole_continuity=False, pole_values=None,
2105 pole_exact=False, pole_flat=False):
2106 iopt = np.array([0, 0, 0], dtype=dfitpack_int)
2107 ider = np.array([-1, 0, -1, 0], dtype=dfitpack_int)
2108 if pole_values is None:
2109 pole_values = (None, None)
2110 elif isinstance(pole_values, (float, np.float32, np.float64)):
2111 pole_values = (pole_values, pole_values)
2112 if isinstance(pole_continuity, bool):
2113 pole_continuity = (pole_continuity, pole_continuity)
2114 if isinstance(pole_exact, bool):
2115 pole_exact = (pole_exact, pole_exact)
2116 if isinstance(pole_flat, bool):
2117 pole_flat = (pole_flat, pole_flat)
2119 r0, r1 = pole_values
2120 iopt[1:] = pole_continuity
2121 if r0 is None:
2122 ider[0] = -1
2123 else:
2124 ider[0] = pole_exact[0]
2126 if r1 is None:
2127 ider[2] = -1
2128 else:
2129 ider[2] = pole_exact[1]
2131 ider[1], ider[3] = pole_flat
2133 u, v = np.ravel(u), np.ravel(v)
2134 r = np.asarray(r)
2136 if not (0.0 < u[0] and u[-1] < np.pi):
2137 raise ValueError('u should be between (0, pi)')
2138 if not -np.pi <= v[0] < np.pi:
2139 raise ValueError('v[0] should be between [-pi, pi)')
2140 if not v[-1] <= v[0] + 2*np.pi:
2141 raise ValueError('v[-1] should be v[0] + 2pi or less ')
2143 if not np.all(np.diff(u) > 0.0):
2144 raise ValueError('u must be strictly increasing')
2145 if not np.all(np.diff(v) > 0.0):
2146 raise ValueError('v must be strictly increasing')
2148 if not u.size == r.shape[0]:
2149 raise ValueError('u dimension of r must have same number of '
2150 'elements as u')
2151 if not v.size == r.shape[1]:
2152 raise ValueError('v dimension of r must have same number of '
2153 'elements as v')
2155 if pole_continuity[1] is False and pole_flat[1] is True:
2156 raise ValueError('if pole_continuity is False, so must be '
2157 'pole_flat')
2158 if pole_continuity[0] is False and pole_flat[0] is True:
2159 raise ValueError('if pole_continuity is False, so must be '
2160 'pole_flat')
2162 if not s >= 0.0:
2163 raise ValueError('s should be positive')
2165 r = np.ravel(r)
2166 nu, tu, nv, tv, c, fp, ier = dfitpack.regrid_smth_spher(iopt, ider,
2167 u.copy(),
2168 v.copy(),
2169 r.copy(),
2170 r0, r1, s)
2172 if ier not in [0, -1, -2]:
2173 msg = _spfit_messages.get(ier, 'ier=%s' % (ier))
2174 raise ValueError(msg)
2176 self.fp = fp
2177 self.tck = tu[:nu], tv[:nv], c[:(nu - 4) * (nv-4)]
2178 self.degrees = (3, 3)
2179 self.v0 = v[0]
2181 def __call__(self, theta, phi, dtheta=0, dphi=0, grid=True):
2183 theta = np.asarray(theta)
2184 phi = np.asarray(phi)
2186 return SphereBivariateSpline.__call__(self, theta, phi, dtheta=dtheta,
2187 dphi=dphi, grid=grid)