Permalink
Browse files

Merge branch 'pull-117-spgrid' into master.

This was reviewed in #117.

Conflicts:
	scipy/interpolate/tests/test_fitpack.py
  • Loading branch information...
2 parents 438e7cc + 3b79227 commit b9081a8cc1258fc62d5fa3fffad2233e8fe17644 @rgommers rgommers committed Mar 10, 2012
View
@@ -98,6 +98,8 @@ Christoph Gohlke for many bug fixes and help with Windows specific issues.
Jacob Silterra for cwt-based peak finding in scipy.signal.
Denis Laxalde for the unified interface to minimizers in scipy.optimize.
David Fong for the sparse LSMR solver.
+Andreas Hilboll for wrapping FITPACK's spgrid.f
+
Institutions
------------
@@ -100,6 +100,7 @@
:toctree: generated/
RectBivariateSpline
+ RectSpherBivariateSpline
For unstructured data:
@@ -15,7 +15,10 @@
'BivariateSpline',
'LSQBivariateSpline',
'SmoothBivariateSpline',
- 'RectBivariateSpline']
+ 'RectBivariateSpline',
+ 'RectSpherBivariateSpline']
+
+from types import NoneType
import warnings
from numpy import zeros, concatenate, alltrue, ravel, all, diff, array
@@ -24,6 +27,7 @@
import fitpack
import dfitpack
+
################ Univariate spline ####################
_curfit_messages = {1:"""
@@ -53,6 +57,7 @@
xb<t[k+1]<t[k+2]<...<t[n-k-2]<xe"""
}
+
class UnivariateSpline(object):
"""
One-dimensional smoothing spline fit to a given set of data points.
@@ -268,6 +273,7 @@ def roots(self):
raise NotImplementedError('finding roots unsupported for '
'non-cubic splines')
+
class InterpolatedUnivariateSpline(UnivariateSpline):
"""
One-dimensional interpolating spline for a given set of data points.
@@ -337,6 +343,7 @@ def __init__(self, x, y, w=None, bbox = [None]*2, k=3):
xb=bbox[0],xe=bbox[1],s=0)
self._reset_class()
+
class LSQUnivariateSpline(UnivariateSpline):
"""
One-dimensional spline with explicit internal knots.
@@ -554,6 +561,7 @@ def integral(self, xa, xb, ya, yb):
kx,ky = self.degrees
return dfitpack.dblint(tx,ty,c,kx,ky,xa,xb,ya,yb)
+
class SmoothBivariateSpline(BivariateSpline):
""" Smooth bivariate spline approximation.
@@ -608,6 +616,7 @@ def __init__(self, x, y, z, w=None, bbox = [None]*4, kx=3, ky=3, s=None,
self.tck = tx[:nx],ty[:ny],c[:(nx-kx-1)*(ny-ky-1)]
self.degrees = kx,ky
+
class LSQBivariateSpline(BivariateSpline):
""" Weighted least-squares bivariate spline approximation.
@@ -644,6 +653,7 @@ class LSQBivariateSpline(BivariateSpline):
bisplrep, bisplev : an older wrapping of FITPACK
UnivariateSpline : a similar class for univariate spline interpolation
SmoothUnivariateSpline : To create a BivariateSpline through the given points
+
"""
def __init__(self, x, y, z, tx, ty, w=None, bbox=[None]*4, kx=3, ky=3,
@@ -676,6 +686,7 @@ def __init__(self, x, y, z, tx, ty, w=None, bbox=[None]*4, kx=3, ky=3,
self.tck = tx1,ty1,c
self.degrees = kx,ky
+
class RectBivariateSpline(BivariateSpline):
""" Bivariate spline approximation over a rectangular mesh.
@@ -703,10 +714,10 @@ class RectBivariateSpline(BivariateSpline):
SmoothBivariateSpline : a smoothing bivariate spline for scattered data
bisplrep, bisplev : an older wrapping of FITPACK
UnivariateSpline : a similar class for univariate spline interpolation
- """
+ """
def __init__(self, x, y, z, bbox = [None]*4, kx=3, ky=3, s=0):
- x,y = ravel(x),ravel(y)
+ x,y = ravel(x), ravel(y)
if not all(diff(x) > 0.0):
raise TypeError('x must be strictly increasing')
if not all(diff(y) > 0.0):
@@ -726,12 +737,227 @@ def __init__(self, x, y, z, bbox = [None]*4, kx=3, ky=3, s=0):
nx,tx,ny,ty,c,fp,ier = dfitpack.regrid_smth(x,y,z,
xb,xe,yb,ye,
kx,ky,s)
- if ier in [0,-1,-2]: # normal return
- pass
- else:
- message = _surfit_messages.get(ier,'ier=%s' % (ier))
- warnings.warn(message)
+
+ if not ier in [0,-1,-2]:
+ msg = _surfit_messages.get(ier, 'ier=%s' % (ier))
+ raise ValueError(msg)
self.fp = fp
self.tck = tx[:nx],ty[:ny],c[:(nx-kx-1)*(ny-ky-1)]
self.degrees = kx,ky
+
+
+_spfit_messages = _surfit_messages.copy()
+_spfit_messages[10] = """
+ERROR: on entry, the input data are controlled on validity
+ the following restrictions must be satisfied.
+ -1<=iopt(1)<=1, 0<=iopt(2)<=1, 0<=iopt(3)<=1,
+ -1<=ider(1)<=1, 0<=ider(2)<=1, ider(2)=0 if iopt(2)=0.
+ -1<=ider(3)<=1, 0<=ider(4)<=1, ider(4)=0 if iopt(3)=0.
+ mu >= mumin (see above), mv >= 4, nuest >=8, nvest >= 8,
+ kwrk>=5+mu+mv+nuest+nvest,
+ lwrk >= 12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+max(nuest,mv+nvest)
+ 0< u(i-1)<u(i)< pi,i=2,..,mu,
+ -pi<=v(1)< pi, v(1)<v(i-1)<v(i)<v(1)+2*pi, i=3,...,mv
+ if iopt(1)=-1: 8<=nu<=min(nuest,mu+6+iopt(2)+iopt(3))
+ 0<tu(5)<tu(6)<...<tu(nu-4)< pi
+ 8<=nv<=min(nvest,mv+7)
+ v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(1)+2*pi
+ the schoenberg-whitney conditions, i.e. there must be
+ subset of grid co-ordinates uu(p) and vv(q) such that
+ tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4
+ (iopt(2)=1 and iopt(3)=1 also count for a uu-value
+ tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4
+ (vv(q) is either a value v(j) or v(j)+2*pi)
+ if iopt(1)>=0: s>=0
+ if s=0: nuest>=mu+6+iopt(2)+iopt(3), nvest>=mv+7
+ if one of these conditions is found to be violated,control is
+ immediately repassed to the calling program. in that case there is no
+ approximation returned."""
+
+
+class RectSpherBivariateSpline(BivariateSpline):
+ """
+ Bivariate spline approximation over a rectangular mesh on a sphere.
+
+ Can be used for smoothing data.
+
+ Parameters
+ ----------
+ 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, and must lie within an interval of
+ length 2pi which in turn lies within the interval (-pi, 2pi).
+ 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.
+
+ See Also
+ --------
+ RectBivariateSpline : bivariate spline approximation over a rectangular mesh
+
+ Notes
+ -----
+ 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.
+
+ .. _FITPACK: http://www.netlib.org/dierckx/spgrid.f
+
+ Examples
+ --------
+ 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 RectSpherBivariateSpline
+ >>> lut = RectSpherBivariateSpline(lats, lons, data)
+
+ Finally we interpolate the data. The `RectSpherBivariateSpline` 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:
+
+ >>> 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 xrange(len(s)):
+ >>> lut = RectSpherBivariateSpline(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()
+
+ """
+ def __init__(self, u, v, r, s=0., pole_continuity=False, pole_values=None,
+ pole_exact=False, pole_flat=False):
+ iopt = np.array([0, 0, 0], dtype=int)
+ ider = np.array([-1, 0, -1, 0], dtype=int)
+ if isinstance(pole_values, NoneType):
+ pole_values = (None, None)
+ elif isinstance(pole_values, (float, np.float32, np.float64)):
+ pole_values = (pole_values, pole_values)
+ if isinstance(pole_continuity, bool):
+ pole_continuity = (pole_continuity, pole_continuity)
+ if isinstance(pole_exact, bool):
+ pole_exact = (pole_exact, pole_exact)
+ if isinstance(pole_flat, bool):
+ pole_flat = (pole_flat, pole_flat)
+
+ r0, r1 = pole_values
+ iopt[1:] = pole_continuity
+ if r0 is None:
+ ider[0] = -1
+ else:
+ ider[0] = pole_exact[0]
+
+ if r1 is None:
+ ider[2] = -1
+ else:
+ ider[2] = pole_exact[1]
+
+ ider[1], ider[3] = pole_flat
+
+ u, v = np.ravel(u), np.ravel(v)
+ if not np.all(np.diff(u) > 0.0):
+ raise TypeError('u must be strictly increasing')
+ if not np.all(np.diff(v) > 0.0):
+ raise TypeError('v must be strictly increasing')
+
+ if not u.size == r.shape[0]:
+ raise TypeError('u dimension of r must have same number of '
+ 'elements as u')
+ if not v.size == r.shape[1]:
+ raise TypeError('v dimension of r must have same number of '
+ 'elements as v')
+
+ if pole_continuity[1] is False and pole_flat[1] is True:
+ raise TypeError('if pole_continuity is False, so must be pole_flat')
+ if pole_continuity[0] is False and pole_flat[0] is True:
+ raise TypeError('if pole_continuity is False, so must be pole_flat')
+
+ r = np.ravel(r)
+ nu, tu, nv, tv, c, fp, ier = dfitpack.regrid_smth_spher(iopt, ider,
+ u.copy(), v.copy(), r.copy(), r0, r1, s)
+
+ if not ier in [0,-1,-2]:
+ msg = _spfit_messages.get(ier, 'ier=%s' % (ier))
+ raise ValueError(msg)
+
+ self.fp = fp
+ self.tck = tu[:nu], tv[:nv], c[:(nu - 4) * (nv-4)]
+ self.degrees = (3, 3)
+
@@ -478,6 +478,42 @@ static int calc_regrid_lwrk(int mx, int my, int kx, int ky,
integer intent(out) :: ier
end subroutine regrid_smth
+ subroutine regrid_smth_spher(iopt,ider,mu,u,mv,v,r,r0,r1,s,nuest,nvest,&
+ nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
+ ! nu,tu,nv,tv,c,fp,ier = regrid_smth_spher(u,v,r,[r0,r1,s])
+
+ fortranname spgrid
+
+ integer dimension(3) :: iopt
+ integer dimension(4) :: ider
+ integer intent(hide) :: mu=len(u)
+ real*8 dimension(mu) :: u
+ integer intent(hide) :: mv=len(v)
+ real*8 dimension(mv) :: v
+ real*8 dimension(mu*mv),depend(mu,mv),check(len(r)==mu*mv) :: r
+ real*8 optional :: r0
+ real*8 optional :: r1
+ real*8 optional,check(0.0<=s) :: s = 0.0
+ integer intent(hide),depend(mu),check(nuest>=8) &
+ :: nuest = mu+6
+ integer intent(hide),depend(mv),check(nvest>=8) &
+ :: nvest = mv+7
+ integer intent(out) :: nu
+ real*8 dimension(nuest),intent(out),depend(nuest) :: tu
+ integer intent(out) :: nv
+ real*8 dimension(nvest),intent(out),depend(nvest) :: tv
+ real*8 dimension((nuest-4)*(nvest-4)), &
+ depend(nuest,nvest),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk),intent(cache,hide),depend(lwrk) :: wrk
+ integer intent(hide),depend(mu,mv,nuest,nvest) &
+ :: lwrk=12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+max(mv+nvest,nuest)
+ integer dimension(kwrk),intent(cache,hide),depend(kwrk) :: iwrk
+ integer intent(hide),depend(mu,mv,nuest,nvest) &
+ :: kwrk=5+mu+mv+nuest+nvest
+ integer intent(out) :: ier
+ end subroutine regrid_smth_spher
+
function dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
! iy = dblint(tx,ty,c,kx,ky,xb,xe,yb,ye)
real*8 dimension(nx),intent(in) :: tx
@@ -218,6 +218,7 @@ def test_smoke_bisplrep_bisplev(self):
put("***************** bisplev")
self.check_5()
+
if __name__ == "__main__":
__put_prints = True
import nose
Oops, something went wrong.

0 comments on commit b9081a8

Please sign in to comment.