interpolate.fitpack2

fitpack — curve and surface fitting with splines

fitpack is based on a collection of Fortran routines DIERCKX by P. Dierckx (see http://www.netlib.org/dierckx/) transformed to double routines by Pearu Peterson.

Module Contents

Classes

UnivariateSpline(self,x,y,w=None,bbox=None,k=3,s=None,ext=0,check_finite=False) One-dimensional smoothing spline fit to a given set of data points.
InterpolatedUnivariateSpline(self,x,y,w=None,bbox=None,k=3,ext=0,check_finite=False) One-dimensional interpolating spline for a given set of data points.
LSQUnivariateSpline(self,x,y,t,w=None,bbox=None,k=3,ext=0,check_finite=False) One-dimensional spline with explicit internal knots.
_BivariateSplineBase() Base class for Bivariate spline s(x,y) interpolation on the rectangle
BivariateSpline() Base class for bivariate splines.
SmoothBivariateSpline(self,x,y,z,w=None,bbox=None,kx=3,ky=3,s=None,eps=None) Smooth bivariate spline approximation.
LSQBivariateSpline(self,x,y,z,tx,ty,w=None,bbox=None,kx=3,ky=3,eps=None) Weighted least-squares bivariate spline approximation.
RectBivariateSpline(self,x,y,z,bbox=None,kx=3,ky=3,s=0) Bivariate spline approximation over a rectangular mesh.
SphereBivariateSpline() Bivariate spline s(x,y) of degrees 3 on a sphere, calculated from a
SmoothSphereBivariateSpline(self,theta,phi,r,w=None,s=0.0,eps=1e-16) Smooth bivariate spline approximation in spherical coordinates.
LSQSphereBivariateSpline(self,theta,phi,r,tt,tp,w=None,eps=1e-16) Weighted least-squares bivariate spline approximation in spherical
RectSphereBivariateSpline(self,u,v,r,s=0.0,pole_continuity=False,pole_values=None,pole_exact=False,pole_flat=False) Bivariate spline approximation over a rectangular mesh on a sphere.
class UnivariateSpline(x, y, w=None, bbox=None, k=3, s=None, ext=0, check_finite=False)

One-dimensional smoothing spline fit to a given set of data points.

Fits a spline y = spl(x) of degree k to the provided x, y data. s specifies the number of knots by specifying a smoothing condition.

x : (N,) array_like
1-D array of independent input data. Must be increasing.
y : (N,) array_like
1-D array of dependent input data, of the same length as x.
w : (N,) array_like, optional
Weights for spline fitting. Must be positive. If None (default), weights are all equal.
bbox : (2,) array_like, optional
2-sequence specifying the boundary of the approximation interval. If None (default), bbox=[x[0], x[-1]].
k : int, optional
Degree of the smoothing spline. Must be <= 5. Default is k=3, a cubic spline.
s : float or None, optional

Positive smoothing factor used to choose the number of knots. Number of knots will be increased until the smoothing condition is satisfied:

sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s

If None (default), s = len(w) which should be a good value if 1/w[i] is an estimate of the standard deviation of y[i]. If 0, spline will interpolate through all data points.

ext : int or str, optional

Controls the extrapolation mode for elements not in the interval defined by the knot sequence.

  • if ext=0 or ‘extrapolate’, return the extrapolated value.
  • if ext=1 or ‘zeros’, return 0
  • if ext=2 or ‘raise’, raise a ValueError
  • if ext=3 of ‘const’, return the boundary value.

The default value is 0.

check_finite : bool, optional
Whether to check that the input arrays contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination or non-sensical results) if the inputs do contain infinities or NaNs. Default is False.

InterpolatedUnivariateSpline : Subclass with smoothing forced to 0 LSQUnivariateSpline : Subclass in which knots are user-selected instead of

being set by smoothing condition

splrep : An older, non object-oriented wrapping of FITPACK splev, sproot, splint, spalde BivariateSpline : A similar class for two-dimensional spline interpolation

The number of data points must be larger than the spline degree k.

NaN handling: If the input arrays contain nan values, the result is not useful, since the underlying spline fitting routines cannot deal with nan . A workaround is to use zero weights for not-a-number data points:

>>> from scipy.interpolate import UnivariateSpline
>>> x, y = np.array([1, 2, 3, 4]), np.array([1, np.nan, 3, 4])
>>> w = np.isnan(y)
>>> y[w] = 0.
>>> spl = UnivariateSpline(x, y, w=~w)

Notice the need to replace a nan by a numerical value (precise value does not matter as long as the corresponding weight is zero.)

>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(-3, 3, 50)
>>> y = np.exp(-x**2) + 0.1 * np.random.randn(50)
>>> plt.plot(x, y, 'ro', ms=5)

Use the default value for the smoothing parameter:

>>> spl = UnivariateSpline(x, y)
>>> xs = np.linspace(-3, 3, 1000)
>>> plt.plot(xs, spl(xs), 'g', lw=3)

Manually change the amount of smoothing:

>>> spl.set_smoothing_factor(0.5)
>>> plt.plot(xs, spl(xs), 'b', lw=3)
>>> plt.show()
__init__(x, y, w=None, bbox=None, k=3, s=None, ext=0, check_finite=False)
_from_tck(tck, ext=0)

Construct a spline object from given tck

_reset_class()
_set_class(cls)
_reset_nest(data, nest=None)
set_smoothing_factor(s)

Continue spline computation with the given smoothing factor s and with the knots found at the last call.

This routine modifies the spline in place.

__call__()

Evaluate spline (or its nu-th derivative) at positions x.

x : array_like
A 1-D array of points at which to return the value of the smoothed spline or its derivatives. Note: x can be unordered but the evaluation is more efficient if x is (partially) ordered.
nu : int
The order of derivative of the spline to compute.
ext : int

Controls the value returned for elements of x not in the interval defined by the knot sequence.

  • if ext=0 or ‘extrapolate’, return the extrapolated value.
  • if ext=1 or ‘zeros’, return 0
  • if ext=2 or ‘raise’, raise a ValueError
  • if ext=3 or ‘const’, return the boundary value.

The default value is 0, passed from the initialization of UnivariateSpline.

get_knots()

Return positions of interior knots of the spline.

Internally, the knot vector contains 2*k additional boundary knots.

get_coeffs()

Return spline coefficients.

get_residual()

Return weighted sum of squared residuals of the spline approximation.

This is equivalent to:

sum((w[i] * (y[i]-spl(x[i])))**2, axis=0)
integral(a, b)

Return definite integral of the spline between two given points.

a : float
Lower limit of integration.
b : float
Upper limit of integration.
integral : float
The value of the definite integral of the spline between limits.
>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(0, 3, 11)
>>> y = x**2
>>> spl = UnivariateSpline(x, y)
>>> spl.integral(0, 3)
9.0

which agrees with between the limits of 0 and 3.

A caveat is that this routine assumes the spline to be zero outside of the data limits:

>>> spl.integral(-1, 4)
9.0
>>> spl.integral(-1, 0)
0.0
derivatives(x)

Return all derivatives of the spline at the point x.

x : float
The point to evaluate the derivatives at.
der : ndarray, shape(k+1,)
Derivatives of the orders 0 to k.
>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(0, 3, 11)
>>> y = x**2
>>> spl = UnivariateSpline(x, y)
>>> spl.derivatives(1.5)
array([2.25, 3.0, 2.0, 0])
roots()

Return the zeros of the spline.

Restriction: only cubic splines are supported by fitpack.

derivative(n=1)

Construct a new spline representing the derivative of this spline.

n : int, optional
Order of derivative to evaluate. Default: 1
spline : UnivariateSpline
Spline of order k2=k-n representing the derivative of this spline.

splder, antiderivative

New in version 0.13.0.

This can be used for finding maxima of a curve:

>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(0, 10, 70)
>>> y = np.sin(x)
>>> spl = UnivariateSpline(x, y, k=4, s=0)

Now, differentiate the spline and find the zeros of the derivative. (NB: sproot only works for order 3 splines, so we fit an order 4 spline):

>>> spl.derivative().roots() / np.pi
array([ 0.50000001,  1.5       ,  2.49999998])

This agrees well with roots of .

antiderivative(n=1)

Construct a new spline representing the antiderivative of this spline.

n : int, optional
Order of antiderivative to evaluate. Default: 1
spline : UnivariateSpline
Spline of order k2=k+n representing the antiderivative of this spline.

New in version 0.13.0.

splantider, derivative

>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(0, np.pi/2, 70)
>>> y = 1 / np.sqrt(1 - 0.8*np.sin(x)**2)
>>> spl = UnivariateSpline(x, y, s=0)

The derivative is the inverse operation of the antiderivative, although some floating point error accumulates:

>>> spl(1.7), spl.antiderivative().derivative()(1.7)
(array(2.1565429877197317), array(2.1565429877201865))

Antiderivative can be used to evaluate definite integrals:

>>> ispl = spl.antiderivative()
>>> ispl(np.pi/2) - ispl(0)
2.2572053588768486

This is indeed an approximation to the complete elliptic integral :

>>> from scipy.special import ellipk
>>> ellipk(0.8)
2.2572053268208538
class InterpolatedUnivariateSpline(x, y, w=None, bbox=None, k=3, ext=0, check_finite=False)

One-dimensional interpolating spline for a given set of data points.

Fits a spline y = spl(x) of degree k to the provided x, y data. Spline function passes through all provided points. Equivalent to UnivariateSpline with s=0.

x : (N,) array_like
Input dimension of data points – must be increasing
y : (N,) array_like
input dimension of data points
w : (N,) array_like, optional
Weights for spline fitting. Must be positive. If None (default), weights are all equal.
bbox : (2,) array_like, optional
2-sequence specifying the boundary of the approximation interval. If None (default), bbox=[x[0], x[-1]].
k : int, optional
Degree of the smoothing spline. Must be 1 <= k <= 5.
ext : int or str, optional

Controls the extrapolation mode for elements not in the interval defined by the knot sequence.

  • if ext=0 or ‘extrapolate’, return the extrapolated value.
  • if ext=1 or ‘zeros’, return 0
  • if ext=2 or ‘raise’, raise a ValueError
  • if ext=3 of ‘const’, return the boundary value.

The default value is 0.

check_finite : bool, optional
Whether to check that the input arrays contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination or non-sensical results) if the inputs do contain infinities or NaNs. Default is False.
UnivariateSpline : Superclass – allows knots to be selected by a
smoothing condition

LSQUnivariateSpline : spline for which knots are user-selected splrep : An older, non object-oriented wrapping of FITPACK splev, sproot, splint, spalde BivariateSpline : A similar class for two-dimensional spline interpolation

The number of data points must be larger than the spline degree k.

>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import InterpolatedUnivariateSpline
>>> x = np.linspace(-3, 3, 50)
>>> y = np.exp(-x**2) + 0.1 * np.random.randn(50)
>>> spl = InterpolatedUnivariateSpline(x, y)
>>> plt.plot(x, y, 'ro', ms=5)
>>> xs = np.linspace(-3, 3, 1000)
>>> plt.plot(xs, spl(xs), 'g', lw=3, alpha=0.7)
>>> plt.show()

Notice that the spl(x) interpolates y:

>>> spl.get_residual()
0.0
__init__(x, y, w=None, bbox=None, k=3, ext=0, check_finite=False)
class LSQUnivariateSpline(x, y, t, w=None, bbox=None, k=3, ext=0, check_finite=False)

One-dimensional spline with explicit internal knots.

Fits a spline y = spl(x) of degree k to the provided x, y data. t specifies the internal knots of the spline

x : (N,) array_like
Input dimension of data points – must be increasing
y : (N,) array_like
Input dimension of data points
t : (M,) array_like

interior knots of the spline. Must be in ascending order and:

bbox[0] < t[0] < ... < t[-1] < bbox[-1]
w : (N,) array_like, optional
weights for spline fitting. Must be positive. If None (default), weights are all equal.
bbox : (2,) array_like, optional
2-sequence specifying the boundary of the approximation interval. If None (default), bbox = [x[0], x[-1]].
k : int, optional
Degree of the smoothing spline. Must be 1 <= k <= 5. Default is k=3, a cubic spline.
ext : int or str, optional

Controls the extrapolation mode for elements not in the interval defined by the knot sequence.

  • if ext=0 or ‘extrapolate’, return the extrapolated value.
  • if ext=1 or ‘zeros’, return 0
  • if ext=2 or ‘raise’, raise a ValueError
  • if ext=3 of ‘const’, return the boundary value.

The default value is 0.

check_finite : bool, optional
Whether to check that the input arrays contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination or non-sensical results) if the inputs do contain infinities or NaNs. Default is False.
ValueError
If the interior knots do not satisfy the Schoenberg-Whitney conditions
UnivariateSpline : Superclass – knots are specified by setting a
smoothing condition

InterpolatedUnivariateSpline : spline passing through all points splrep : An older, non object-oriented wrapping of FITPACK splev, sproot, splint, spalde BivariateSpline : A similar class for two-dimensional spline interpolation

The number of data points must be larger than the spline degree k.

Knots t must satisfy the Schoenberg-Whitney conditions, i.e., there must be a subset of data points x[j] such that t[j] < x[j] < t[j+k+1], for j=0, 1,...,n-k-2.

>>> from scipy.interpolate import LSQUnivariateSpline, UnivariateSpline
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(-3, 3, 50)
>>> y = np.exp(-x**2) + 0.1 * np.random.randn(50)

Fit a smoothing spline with a pre-defined internal knots:

>>> t = [-1, 0, 1]
>>> spl = LSQUnivariateSpline(x, y, t)
>>> xs = np.linspace(-3, 3, 1000)
>>> plt.plot(x, y, 'ro', ms=5)
>>> plt.plot(xs, spl(xs), 'g-', lw=3)
>>> plt.show()

Check the knot vector:

>>> spl.get_knots()
array([-3., -1., 0., 1., 3.])

Constructing lsq spline using the knots from another spline:

>>> x = np.arange(10)
>>> s = UnivariateSpline(x, x, s=0)
>>> s.get_knots()
array([ 0.,  2.,  3.,  4.,  5.,  6.,  7.,  9.])
>>> knt = s.get_knots()
>>> s1 = LSQUnivariateSpline(x, x, knt[1:-1])    # Chop 1st and last knot
>>> s1.get_knots()
array([ 0.,  2.,  3.,  4.,  5.,  6.,  7.,  9.])
__init__(x, y, t, w=None, bbox=None, k=3, ext=0, check_finite=False)
class _BivariateSplineBase

Base class for Bivariate spline s(x,y) interpolation on the rectangle [xb,xe] x [yb, ye] calculated from a given set of data points (x,y,z).

bisplrep, bisplev : an older wrapping of FITPACK BivariateSpline :

implementation of bivariate spline interpolation on a plane grid
SphereBivariateSpline :
implementation of bivariate spline interpolation on a spherical grid
get_residual()

Return weighted sum of squared residuals of the spline approximation: sum ((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0)

get_knots()

Return a tuple (tx,ty) where tx,ty contain knots positions of the spline with respect to x-, y-variable, respectively. The position of interior and additional knots are given as t[k+1:-k-1] and t[:k+1]=b, t[-k-1:]=e, respectively.

get_coeffs()

Return spline coefficients.

__call__(x, y, dx=0, dy=0, grid=True)

Evaluate the spline or its derivatives at given positions.

x, y : array_like

Input coordinates.

If grid is False, evaluate the spline at points (x[i], y[i]), i=0, ..., len(x)-1. Standard Numpy broadcasting is obeyed.

If grid is True: evaluate spline at the grid points defined by the coordinate arrays x, y. The arrays must be sorted to increasing order.

dx : int

Order of x-derivative

New in version 0.14.0.

dy : int

Order of y-derivative

New in version 0.14.0.

grid : bool

Whether to evaluate the results on a grid spanned by the input arrays, or at points specified by the input arrays.

New in version 0.14.0.

class BivariateSpline

Base class for bivariate splines.

This describes a spline s(x, y) of degrees kx and ky on the rectangle [xb, xe] * [yb, ye] calculated from a given set of data points (x, y, z).

This class is meant to be subclassed, not instantiated directly. To construct these splines, call either SmoothBivariateSpline or LSQBivariateSpline.

UnivariateSpline : a similar class for univariate spline interpolation SmoothBivariateSpline :

to create a BivariateSpline through the given points
LSQBivariateSpline :
to create a BivariateSpline using weighted least-squares fitting
SphereBivariateSpline :
bivariate spline interpolation in spherical cooridinates

bisplrep : older wrapping of FITPACK bisplev : older wrapping of FITPACK

_from_tck(tck)

Construct a spline object from given tck and degree

ev(xi, yi, dx=0, dy=0)

Evaluate the spline at points

Returns the interpolated value at (xi[i], yi[i]), i=0,...,len(xi)-1.

xi, yi : array_like
Input coordinates. Standard Numpy broadcasting is obeyed.
dx : int, optional

Order of x-derivative

New in version 0.14.0.

dy : int, optional

Order of y-derivative

New in version 0.14.0.

integral(xa, xb, ya, yb)

Evaluate the integral of the spline over area [xa,xb] x [ya,yb].

xa, xb : float
The end-points of the x integration interval.
ya, yb : float
The end-points of the y integration interval.
integ : float
The value of the resulting integral.
class SmoothBivariateSpline(x, y, z, w=None, bbox=None, kx=3, ky=3, s=None, eps=None)

Smooth bivariate spline approximation.

x, y, z : array_like
1-D sequences of data points (order is not important).
w : array_like, optional
Positive 1-D sequence of weights, of same length as x, y and z.
bbox : array_like, optional
Sequence of length 4 specifying the boundary of the rectangular approximation domain. By default, bbox=[min(x,tx),max(x,tx), min(y,ty),max(y,ty)].
kx, ky : ints, optional
Degrees of the bivariate spline. Default is 3.
s : float, optional
Positive smoothing factor defined for estimation condition: sum((w[i]*(z[i]-s(x[i], y[i])))**2, axis=0) <= s Default s=len(w) which should be a good value if 1/w[i] is an estimate of the standard deviation of z[i].
eps : float, optional
A threshold for determining the effective rank of an over-determined linear system of equations. eps should have a value between 0 and 1, the default is 1e-16.

bisplrep : an older wrapping of FITPACK bisplev : an older wrapping of FITPACK UnivariateSpline : a similar class for univariate spline interpolation LSQUnivariateSpline : to create a BivariateSpline using weighted

The length of x, y and z should be at least (kx+1) * (ky+1).

__init__(x, y, z, w=None, bbox=None, kx=3, ky=3, s=None, eps=None)
class LSQBivariateSpline(x, y, z, tx, ty, w=None, bbox=None, kx=3, ky=3, eps=None)

Weighted least-squares bivariate spline approximation.

x, y, z : array_like
1-D sequences of data points (order is not important).
tx, ty : array_like
Strictly ordered 1-D sequences of knots coordinates.
w : array_like, optional
Positive 1-D array of weights, of the same length as x, y and z.
bbox : (4,) array_like, optional
Sequence of length 4 specifying the boundary of the rectangular approximation domain. By default, bbox=[min(x,tx),max(x,tx), min(y,ty),max(y,ty)].
kx, ky : ints, optional
Degrees of the bivariate spline. Default is 3.
eps : float, optional
A threshold for determining the effective rank of an over-determined linear system of equations. eps should have a value between 0 and 1, the default is 1e-16.

bisplrep : an older wrapping of FITPACK bisplev : an older wrapping of FITPACK UnivariateSpline : a similar class for univariate spline interpolation SmoothBivariateSpline : create a smoothing BivariateSpline

The length of x, y and z should be at least (kx+1) * (ky+1).

__init__(x, y, z, tx, ty, w=None, bbox=None, kx=3, ky=3, eps=None)
class RectBivariateSpline(x, y, z, bbox=None, kx=3, ky=3, s=0)

Bivariate spline approximation over a rectangular mesh.

Can be used for both smoothing and interpolating data.

x,y : array_like
1-D arrays of coordinates in strictly ascending order.
z : array_like
2-D array of data with shape (x.size,y.size).
bbox : array_like, optional
Sequence of length 4 specifying the boundary of the rectangular approximation domain. By default, bbox=[min(x,tx),max(x,tx), min(y,ty),max(y,ty)].
kx, ky : ints, optional
Degrees of the bivariate spline. Default is 3.
s : float, optional
Positive smoothing factor defined for estimation condition: sum((w[i]*(z[i]-s(x[i], y[i])))**2, axis=0) <= s Default is s=0, which is for interpolation.

SmoothBivariateSpline : a smoothing bivariate spline for scattered data bisplrep : an older wrapping of FITPACK bisplev : an older wrapping of FITPACK UnivariateSpline : a similar class for univariate spline interpolation

__init__(x, y, z, bbox=None, kx=3, ky=3, s=0)
class SphereBivariateSpline

Bivariate spline s(x,y) of degrees 3 on a sphere, calculated from a given set of data points (theta,phi,r).

New in version 0.11.0.

bisplrep, bisplev : an older wrapping of FITPACK UnivariateSpline : a similar class for univariate spline interpolation SmoothUnivariateSpline :

to create a BivariateSpline through the given points
LSQUnivariateSpline :
to create a BivariateSpline using weighted least-squares fitting
__call__(theta, phi, dtheta=0, dphi=0, grid=True)

Evaluate the spline or its derivatives at given positions.

theta, phi : array_like

Input coordinates.

If grid is False, evaluate the spline at points (theta[i], phi[i]), i=0, ..., len(x)-1. Standard Numpy broadcasting is obeyed.

If grid is True: evaluate spline at the grid points defined by the coordinate arrays theta, phi. The arrays must be sorted to increasing order.

dtheta : int, optional

Order of theta-derivative

New in version 0.14.0.

dphi : int

Order of phi-derivative

New in version 0.14.0.

grid : bool

Whether to evaluate the results on a grid spanned by the input arrays, or at points specified by the input arrays.

New in version 0.14.0.

ev(theta, phi, dtheta=0, dphi=0)

Evaluate the spline at points

Returns the interpolated value at (theta[i], phi[i]), i=0,...,len(theta)-1.

theta, phi : array_like
Input coordinates. Standard Numpy broadcasting is obeyed.
dtheta : int, optional

Order of theta-derivative

New in version 0.14.0.

dphi : int, optional

Order of phi-derivative

New in version 0.14.0.

class SmoothSphereBivariateSpline(theta, phi, r, w=None, s=0.0, eps=1e-16)

Smooth bivariate spline approximation in spherical coordinates.

New in version 0.11.0.

theta, phi, r : array_like
1-D sequences of data points (order is not important). Coordinates must be given in radians. Theta must lie within the interval (0, pi), and phi must lie within the interval (0, 2pi).
w : array_like, optional
Positive 1-D sequence of weights.
s : float, optional
Positive smoothing factor defined for estimation condition: sum((w(i)*(r(i) - s(theta(i), phi(i))))**2, axis=0) <= s Default s=len(w) which should be a good value if 1/w[i] is an estimate of the standard deviation of r[i].
eps : float, optional
A threshold for determining the effective rank of an over-determined linear system of equations. eps should have a value between 0 and 1, the default is 1e-16.

For more information, see the FITPACK_ site about this function.

Suppose we have global data on a coarse grid (the input data does not have to be on a grid):

>>> theta = np.linspace(0., np.pi, 7)
>>> phi = np.linspace(0., 2*np.pi, 9)
>>> data = np.empty((theta.shape[0], phi.shape[0]))
>>> data[:,0], data[0,:], data[-1,:] = 0., 0., 0.
>>> data[1:-1,1], data[1:-1,-1] = 1., 1.
>>> data[1,1:-1], data[-2,1:-1] = 1., 1.
>>> data[2:-2,2], data[2:-2,-2] = 2., 2.
>>> data[2,2:-2], data[-3,2:-2] = 2., 2.
>>> data[3,3:-2] = 3.
>>> data = np.roll(data, 4, 1)

We need to set up the interpolator object

>>> lats, lons = np.meshgrid(theta, phi)
>>> from scipy.interpolate import SmoothSphereBivariateSpline
>>> lut = SmoothSphereBivariateSpline(lats.ravel(), lons.ravel(),
...                                   data.T.ravel(), s=3.5)

As a first test, we’ll see what the algorithm returns when run on the input coordinates

>>> data_orig = lut(theta, phi)

Finally we interpolate the data to a finer grid

>>> fine_lats = np.linspace(0., np.pi, 70)
>>> fine_lons = np.linspace(0., 2 * np.pi, 90)
>>> data_smth = lut(fine_lats, fine_lons)
>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(131)
>>> ax1.imshow(data, interpolation='nearest')
>>> ax2 = fig.add_subplot(132)
>>> ax2.imshow(data_orig, interpolation='nearest')
>>> ax3 = fig.add_subplot(133)
>>> ax3.imshow(data_smth, interpolation='nearest')
>>> plt.show()
__init__(theta, phi, r, w=None, s=0.0, eps=1e-16)
class LSQSphereBivariateSpline(theta, phi, r, tt, tp, w=None, eps=1e-16)

Weighted least-squares bivariate spline approximation in spherical coordinates.

New in version 0.11.0.

theta, phi, r : array_like
1-D sequences of data points (order is not important). Coordinates must be given in radians. Theta must lie within the interval (0, pi), and phi must lie within the interval (0, 2pi).
tt, tp : array_like
Strictly ordered 1-D sequences of knots coordinates. Coordinates must satisfy 0 < tt[i] < pi, 0 < tp[i] < 2*pi.
w : array_like, optional
Positive 1-D sequence of weights, of the same length as theta, phi and r.
eps : float, optional
A threshold for determining the effective rank of an over-determined linear system of equations. eps should have a value between 0 and 1, the default is 1e-16.

For more information, see the FITPACK_ site about this function.

Suppose we have global data on a coarse grid (the input data does not have to be on a grid):

>>> theta = np.linspace(0., np.pi, 7)
>>> phi = np.linspace(0., 2*np.pi, 9)
>>> data = np.empty((theta.shape[0], phi.shape[0]))
>>> data[:,0], data[0,:], data[-1,:] = 0., 0., 0.
>>> data[1:-1,1], data[1:-1,-1] = 1., 1.
>>> data[1,1:-1], data[-2,1:-1] = 1., 1.
>>> data[2:-2,2], data[2:-2,-2] = 2., 2.
>>> data[2,2:-2], data[-3,2:-2] = 2., 2.
>>> data[3,3:-2] = 3.
>>> data = np.roll(data, 4, 1)

We need to set up the interpolator object. Here, we must also specify the coordinates of the knots to use.

>>> lats, lons = np.meshgrid(theta, phi)
>>> knotst, knotsp = theta.copy(), phi.copy()
>>> knotst[0] += .0001
>>> knotst[-1] -= .0001
>>> knotsp[0] += .0001
>>> knotsp[-1] -= .0001
>>> from scipy.interpolate import LSQSphereBivariateSpline
>>> lut = LSQSphereBivariateSpline(lats.ravel(), lons.ravel(),
...                                data.T.ravel(), knotst, knotsp)

As a first test, we’ll see what the algorithm returns when run on the input coordinates

>>> data_orig = lut(theta, phi)

Finally we interpolate the data to a finer grid

>>> fine_lats = np.linspace(0., np.pi, 70)
>>> fine_lons = np.linspace(0., 2*np.pi, 90)
>>> data_lsq = lut(fine_lats, fine_lons)
>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(131)
>>> ax1.imshow(data, interpolation='nearest')
>>> ax2 = fig.add_subplot(132)
>>> ax2.imshow(data_orig, interpolation='nearest')
>>> ax3 = fig.add_subplot(133)
>>> ax3.imshow(data_lsq, interpolation='nearest')
>>> plt.show()
__init__(theta, phi, r, tt, tp, w=None, eps=1e-16)
class RectSphereBivariateSpline(u, v, r, s=0.0, pole_continuity=False, pole_values=None, pole_exact=False, pole_flat=False)

Bivariate spline approximation over a rectangular mesh on a sphere.

Can be used for smoothing data.

New in version 0.11.0.

u : array_like
1-D array of latitude coordinates in strictly ascending order. Coordinates must be given in radians and lie within the interval (0, pi).
v : array_like
1-D array of longitude coordinates in strictly ascending order. Coordinates must be given in radians. First element (v[0]) must lie within the interval [-pi, pi). Last element (v[-1]) must satisfy v[-1] <= v[0] + 2*pi.
r : array_like
2-D array of data with shape (u.size, v.size).
s : float, optional
Positive smoothing factor defined for estimation condition (s=0 is for interpolation).
pole_continuity : bool or (bool, bool), optional
Order of continuity at the poles u=0 (pole_continuity[0]) and u=pi (pole_continuity[1]). The order of continuity at the pole will be 1 or 0 when this is True or False, respectively. Defaults to False.
pole_values : float or (float, float), optional
Data values at the poles u=0 and u=pi. Either the whole parameter or each individual element can be None. Defaults to None.
pole_exact : bool or (bool, bool), optional
Data value exactness at the poles u=0 and u=pi. If True, the value is considered to be the right function value, and it will be fitted exactly. If False, the value will be considered to be a data value just like the other data values. Defaults to False.
pole_flat : bool or (bool, bool), optional
For the poles at u=0 and u=pi, specify whether or not the approximation has vanishing derivatives. Defaults to False.
RectBivariateSpline : bivariate spline approximation over a rectangular
mesh

Currently, only the smoothing spline approximation (iopt[0] = 0 and iopt[0] = 1 in the FITPACK routine) is supported. The exact least-squares spline approximation is not implemented yet.

When actually performing the interpolation, the requested v values must lie within the same length 2pi interval that the original v values were chosen from.

For more information, see the FITPACK_ site about this function.

Suppose we have global data on a coarse grid

>>> lats = np.linspace(10, 170, 9) * np.pi / 180.
>>> lons = np.linspace(0, 350, 18) * np.pi / 180.
>>> data = np.dot(np.atleast_2d(90. - np.linspace(-80., 80., 18)).T,
...               np.atleast_2d(180. - np.abs(np.linspace(0., 350., 9)))).T

We want to interpolate it to a global one-degree grid

>>> new_lats = np.linspace(1, 180, 180) * np.pi / 180
>>> new_lons = np.linspace(1, 360, 360) * np.pi / 180
>>> new_lats, new_lons = np.meshgrid(new_lats, new_lons)

We need to set up the interpolator object

>>> from scipy.interpolate import RectSphereBivariateSpline
>>> lut = RectSphereBivariateSpline(lats, lons, data)

Finally we interpolate the data. The RectSphereBivariateSpline object only takes 1-D arrays as input, therefore we need to do some reshaping.

>>> data_interp = lut.ev(new_lats.ravel(),
...                      new_lons.ravel()).reshape((360, 180)).T

Looking at the original and the interpolated data, one can see that the interpolant reproduces the original data very well:

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(211)
>>> ax1.imshow(data, interpolation='nearest')
>>> ax2 = fig.add_subplot(212)
>>> ax2.imshow(data_interp, interpolation='nearest')
>>> plt.show()

Chosing the optimal value of s can be a delicate task. Recommended values for s depend on the accuracy of the data values. If the user has an idea of the statistical errors on the data, she can also find a proper estimate for s. By assuming that, if she specifies the right s, the interpolator will use a spline f(u,v) which exactly reproduces the function underlying the data, she can evaluate sum((r(i,j)-s(u(i),v(j)))**2) to find a good estimate for this s. For example, if she knows that the statistical errors on her r(i,j)-values are not greater than 0.1, she may expect that a good s should have a value not larger than u.size * v.size * (0.1)**2.

If nothing is known about the statistical error in r(i,j), s must be determined by trial and error. The best is then to start with a very large value of s (to determine the least-squares polynomial and the corresponding upper bound fp0 for s) and then to progressively decrease the value of s (say by a factor 10 in the beginning, i.e. s = fp0 / 10, fp0 / 100, ... and more carefully as the approximation shows more detail) to obtain closer fits.

The interpolation results for different values of s give some insight into this process:

>>> fig2 = plt.figure()
>>> s = [3e9, 2e9, 1e9, 1e8]
>>> for ii in range(len(s)):
...     lut = RectSphereBivariateSpline(lats, lons, data, s=s[ii])
...     data_interp = lut.ev(new_lats.ravel(),
...                          new_lons.ravel()).reshape((360, 180)).T
...     ax = fig2.add_subplot(2, 2, ii+1)
...     ax.imshow(data_interp, interpolation='nearest')
...     ax.set_title("s = %g" % s[ii])
>>> plt.show()
__init__(u, v, r, s=0.0, pole_continuity=False, pole_values=None, pole_exact=False, pole_flat=False)