Wrapped spgrid.f from FITPACK into RectSpherBivariateSpline class #117

Closed
wants to merge 7 commits into
from
View
@@ -95,6 +95,7 @@ Fazlul Shahriar for fixes to the NetCDF3 I/O.
Chris Jordan-Squire for bug fixes, documentation improvements and
scipy.special.logit & expit.
Christoph Gohlke for many bug fixes and help with Windows specific issues.
+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
@@ -644,6 +647,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,
@@ -703,8 +707,8 @@ 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)
if not all(diff(x) > 0.0):
@@ -735,3 +739,202 @@ def __init__(self, x, y, z, bbox = [None]*4, kx=3, ky=3, s=0):
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.
+ For more information, see the FITPACK_ site about this function.
+
+ .. _FITPACK: http://www.netlib.org/dierckx/spgrid.f
+
+ 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``.
+
+ Notice
+ ------
+ 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.
+
+ See Also
+ --------
+ RectBivariateSpline : bivariate spline approximation over a rectangular mesh
+
+ Examples
+ --------
+ Suppose we have global data on a coarse grid
+
+ >>> from numpy import abs, atleast_2d, dot, linspace, meshgrid, pi
+ >>> lats = (linspace(-80., 80., 9) + 90.) * pi / 180.
+ >>> lons = linspace(0., 350., 18) * pi / 180.
+ >>> data = dot(atleast_2d(90. - linspace(-80., 80., 18)).T,
+ atleast_2d(180. - abs(linspace(0., 350., 9)))).T
+
+ We want to interpolate it to a global one-degree grid
+
+ >>> lats = (linspace(-89., 89., 179) + 90.) * pi / 180.
+ >>> lons = linspace(0., 360., 361) * pi / 180.
+ >>> lats, lons = meshgrid(lats, 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 1d arrays as input, therefore we need to do some reshaping.
+
+ >>> data = lut.ev(new_lats.ravel(), new_lons.ravel()).reshape((361, 179)).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()
+ >>> orig = fig.add_subplot(211)
+ >>> orig.pcolor(data_orig)
+ >>> interp = fig.add_subplot(212)
+ >>> interp.pcolor(data)
+
+ 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 the
+ 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:
+
+ >>> figs = plt.figure()
+ >>> s = [3E9, 2E9, 1E9, 1E8]
+ >>> for i in xrange(len(s)):
+ >>> lut = RectSpherBivariateSpline(lats_orig, lons_orig, data_orig, s=s)
+ >>> data = lut.ev(lats.ravel(), lons.ravel()).reshape((361, 179)).T
+ >>> ax = figs.add_subplot(2, 2, i+1)
+ >>> ax.pcolor(data)
+ >>> ax.set_title("s = %g" % s)
+
+ """
+ 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
+ ider[0] = -1 if r0 is None else pole_exact[0]
@pv

pv Feb 25, 2012

Owner

We'd like to retain Python 2.4 compatibility for now, so conditional expressions shouldn't be used.

+ ider[2] = -1 if r1 is None else 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]:
+ message = _spfit_messages.get(ier, 'ier=%s' % (ier))
+ warnings.warn(message)
+ raise ValueError()
@pv

pv Feb 25, 2012

Owner

Here the message should be given to ValueError, rather than in the warning.

@andreas-h

andreas-h Feb 25, 2012

Contributor

That's a copy-past from the RectBivariateSpline class.

Am 25.02.2012 15:13, schrieb Pauli Virtanen:

  •                        '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]:
    
  •        message = _spfit_messages.get(ier, 'ier=%s' % (ier))
    
  •        warnings.warn(message)
    
  •        raise ValueError()
    

Here the message should be given to ValueError, rather than in the warning.


Reply to this email directly or view it on GitHub:
https://github.com/scipy/scipy/pull/117/files#r488632

@andreas-h

andreas-h Feb 25, 2012

Contributor

Sorry, my mistake. It came as a copy-paste, but then I decided it would
be good to raise the ValueError, because the results aren't
meaningful in any way. So yes, the message should go to the exception.

I propose to change the behaviour in the RectBivariateSpline as well
to raise an error, because also there, the results are meaningless if
ier > 0.

Am Sa 25 Feb 2012 18:01:52 CET schrieb Andreas H.:

That's a copy-past from the RectBivariateSpline class.

Am 25.02.2012 15:13, schrieb Pauli Virtanen:

  •                        '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]:
    
  •        message = _spfit_messages.get(ier, 'ier=%s' % (ier))
    
  •        warnings.warn(message)
    
  •        raise ValueError()
    

Here the message should be given to ValueError, rather than in the warning.


Reply to this email directly or view it on GitHub:
https://github.com/scipy/scipy/pull/117/files#r488632

+
+ 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
@@ -14,9 +14,9 @@
from numpy.testing import assert_equal, assert_almost_equal, assert_array_equal, \
assert_array_almost_equal, assert_allclose, TestCase, run_module_suite
-from numpy import array, diff, shape
+from numpy import array, diff, linspace, pi, shape
from scipy.interpolate.fitpack2 import UnivariateSpline, LSQBivariateSpline, \
- SmoothBivariateSpline, RectBivariateSpline
+ SmoothBivariateSpline, RectBivariateSpline, RectSpherBivariateSpline
class TestUnivariateSpline(TestCase):
def test_linear_constant(self):
@@ -211,5 +211,28 @@ def test_evaluate(self):
assert_almost_equal(zi, zi2)
+class TestRectSpherBivariateSpline(TestCase):
+ def test_defaults(self):
+ y = linspace(0.01, 2*pi-0.01, 7)
+ x = linspace(0.01, pi-0.01, 7)
+ z = array([[1,2,1,2,1,2,1],[1,2,1,2,1,2,1],[1,2,3,2,1,2,1],
+ [1,2,2,2,1,2,1],[1,2,1,2,1,2,1],[1,2,2,2,1,2,1],
+ [1,2,1,2,1,2,1]])
+ lut = RectSpherBivariateSpline(x,y,z)
+ assert_array_almost_equal(lut(x,y),z)
+
+ def test_evaluate(self):
+ y = linspace(0.01, 2*pi-0.01, 7)
+ x = linspace(0.01, pi-0.01, 7)
+ z = array([[1,2,1,2,1,2,1],[1,2,1,2,1,2,1],[1,2,3,2,1,2,1],
+ [1,2,2,2,1,2,1],[1,2,1,2,1,2,1],[1,2,2,2,1,2,1],
+ [1,2,1,2,1,2,1]])
+ lut = RectSpherBivariateSpline(x,y,z)
+ yi = [0.2, 1, 2.3, 2.35, 3.0, 3.99, 5.25]
+ xi = [1.5, 0.4, 1.1, 0.45, 0.2345, 1., 0.0001]
+ zi = lut.ev(xi, yi)
+ zi2 = array([lut(xp, yp)[0,0] for xp, yp in zip(xi, yi)])
+ assert_almost_equal(zi, zi2)
+
if __name__ == "__main__":
run_module_suite()