Wrapped spgrid.f from FITPACK into RectSpherBivariateSpline class #117

Closed
wants to merge 7 commits into
from

Projects

None yet

4 participants

@andreas-h
Contributor

According to http://projects.scipy.org/scipy/ticket/1549, here's the pull request. I wrapped FITPACK's spgrid.f function for interpolating on a rectangular grid on a spherical surface for using in scipy.interpolate.RectSpherBivariateSpline.

@rgommers rgommers commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
+
+ See Also from ``scipy.interpolate``
+ -----------------------------------
+ RectBivariateSpline : vivariate spline approximation over a rectangular mesh
+ bisplrep, bisplev : an older wrapping of FITPACK
+
+ """
+ def __init__(self, u, v, r, s=0., iopt=np.array([0, 0, 0], dtype=int),
+ ider=np.array([-1, 0, -1, 0], dtype=int), r0=None, r1=None):
+ 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.min() == u[0]) and (u.max() == u[-1])):
+ raise TypeError('u must be strictly ascending')
@rgommers
rgommers Dec 28, 2011 Member

Aren't this and the next check unnecessary? How can this not be true if np.all(np.diff(u) > 0.0)?

@rgommers rgommers commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
+
+ 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.
+ 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).
+ iopt : nd.array, optional
@rgommers
rgommers Dec 28, 2011 Member

nd.array --> ndarray. Same for ider parameter.

@rgommers rgommers commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
+ data value r1 at the pole u=pi. if ider[2]=1, r1 will be
+ considered to be the right function value, and it will be
+ fitted exactly (s(pi,v)=r1). if ider[2]=0, r1 will be
+ considered to be a data value just like the other data values
+ r[i,j].
+ ider[3]: specify whether (ider[3]=1) or not (ider[3]=0) the
+ approximation has vanishing derivatives dr(5) and dr(6) at the
+ pole u=pi (in case iopt[2]=1).
+ r0 : float, optional
+ specify the data value at the pole u=0
+ r1 : float, optional
+ specify the data value at the pole u=pi
+
+ See Also from ``scipy.interpolate``
+ -----------------------------------
+ RectBivariateSpline : vivariate spline approximation over a rectangular mesh
@rgommers
rgommers Dec 28, 2011 Member

bivariate

@rgommers rgommers commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
+ pole u=0 (in case iopt[1]=1).
+ ider[2]: specify whether (ider[2]=0 or 1) or not (ider[2]=-1) there is a
+ data value r1 at the pole u=pi. if ider[2]=1, r1 will be
+ considered to be the right function value, and it will be
+ fitted exactly (s(pi,v)=r1). if ider[2]=0, r1 will be
+ considered to be a data value just like the other data values
+ r[i,j].
+ ider[3]: specify whether (ider[3]=1) or not (ider[3]=0) the
+ approximation has vanishing derivatives dr(5) and dr(6) at the
+ pole u=pi (in case iopt[2]=1).
+ r0 : float, optional
+ specify the data value at the pole u=0
+ r1 : float, optional
+ specify the data value at the pole u=pi
+
+ See Also from ``scipy.interpolate``
@rgommers
rgommers Dec 28, 2011 Member

Section title can simply be "See Also"

@rgommers rgommers commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
+ '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 not iopt.size == 3:
+ raise TypeError("iopt must be of shape (3,)")
+ if not ider.size == 4:
+ raise TypeError("ider must be of shape (4,)")
+ if ider[0] > -1 and r0 == None:
+ raise TypeError("if ider[0] > -1, you must give r0.")
+ if ider[2] > -1 and r1 == None:
+ raise TypeError("if ider[2] > -1, you must give r1.")
+ r = np.ravel(r)
+ print u
+ print v
+ print r
@rgommers
rgommers Dec 28, 2011 Member

print statements left in

@rgommers rgommers and 1 other commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
+ if not ider.size == 4:
+ raise TypeError("ider must be of shape (4,)")
+ if ider[0] > -1 and r0 == None:
+ raise TypeError("if ider[0] > -1, you must give r0.")
+ if ider[2] > -1 and r1 == None:
+ raise TypeError("if ider[2] > -1, you must give r1.")
+ r = np.ravel(r)
+ print u
+ print v
+ print r
+ nu,tu,nv,tv,c,fp,ier = dfitpack.regrid_smth_spher(iopt,ider,u.copy(),
+ v.copy(),r.copy(),r0,
+ r1,s)
+ if ier in [0,-1,-2]: # normal return
+ pass
+ else:
@rgommers
rgommers Dec 28, 2011 Member

above 3 lines can be written more compact as if not ier in [0,-1,-2]:

@pv
pv Jan 18, 2012 Member

Or even if ier not in [0, -1, -2]:

@rgommers rgommers commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
@@ -735,3 +736,172 @@ 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 = {1:"""
@rgommers
rgommers Dec 28, 2011 Member

The only error message that is substantially different from _surfit_messages is 10. Why not do:

_spfit_messages = _surfit_messages.copy()
-spfit_message[10] = """<msg>"""
@rgommers rgommers commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
+ 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 both smoothing and interpolating data.
+
+ For more information, see the FITPACK_ site about this function.
+
+ .. _FITPACK: http://www.netlib.org/dierckx/spgrid.f
@rgommers
rgommers Dec 28, 2011 Member

Some blank lines can be removed here.

@rgommers rgommers commented on an outdated diff Dec 28, 2011
scipy/interpolate/fitpack2.py
+ considered to be the right function value, and it will be
+ fitted exactly (s(pi,v)=r1). if ider[2]=0, r1 will be
+ considered to be a data value just like the other data values
+ r[i,j].
+ ider[3]: specify whether (ider[3]=1) or not (ider[3]=0) the
+ approximation has vanishing derivatives dr(5) and dr(6) at the
+ pole u=pi (in case iopt[2]=1).
+ r0 : float, optional
+ specify the data value at the pole u=0
+ r1 : float, optional
+ specify the data value at the pole u=pi
+
+ See Also from ``scipy.interpolate``
+ -----------------------------------
+ RectBivariateSpline : vivariate spline approximation over a rectangular mesh
+ bisplrep, bisplev : an older wrapping of FITPACK
@rgommers
rgommers Dec 28, 2011 Member

This line can be removed, because the older wrappers don't wrap this routine.

@rgommers
Member

Looks good to me. Did you call the class RectSpher instead of RectSphere on purpose? The latter seems more natural.

An example would be very useful.

@pv
Member
pv commented Jan 18, 2012

This looks good, overall, and would be an useful addition. There are some things to address, aside what Ralf mentions above: IMO the ider / iopt interface is not really optimal -- passing input parameters in integer arrays is too Fortran-ish.

@andreas-h
Contributor

If you prefer, I can give each of iopt/ider's elements its own input
parameter. I just thought it'd be easier to handle this way ...
Whichever you prefer -- I'm just the newbie ;)

On 2012-01-18 12:23, Pauli Virtanen wrote:

This looks good, overall, and would be an useful addition. There are some things to address, aside what Ralf mentions above: IMO the ider / iopt interface is not really optimal -- passing input parameters in integer arrays is too Fortran-ish.


Reply to this email directly or view it on GitHub:
#117 (comment)

@pv
Member
pv commented Jan 18, 2012

I'd maybe do like this:

  • iopt[0]: could be determined from the value of s: least squares for s = 0 or s = None?
  • pole_continuity = (iopt[1], iopt[2]), as boolean flags.
  • pole_values = (r0, r1). The whole parameter, or r0/r1 separately can be None
  • pole_exact = (ider[0], ider[2]). Also these can be boolean flags -- the -1 is put in if r0/r1 is omitted (None)
  • pole_flat = (ider[1], ider[3]), also boolean flags.
@rgommers
Member

Hi Andreas, are you intending to work on this in the near future? I'd like to get this merged sometime soon, to give it some time to settle before the 0.11.0 release.

@andreas-h
Contributor

@pv as for iopt[0], spgrid.f does not return the same result for s=0, iopt[0]=0 and s=0, iopt[0]=-1. Therefore, I add an additional least_squares parameter. Agreed for the others.

@andreas-h
Contributor

I included everything you two suggested, except for the iopt[0] as stated above.
Only the usage example is missing from the docstring, will do that later, hopefully today.

@andreas-h
Contributor

Yes, I did call it RectSpher on purpose (instead of RectSphere), as I was thinking of spherical.

@andreas-h
Contributor

While working on the usage example I found some problems with my implementation, which I won't be able to fix within the next week or so. So don't wait for me with 0.11.0, please.

@rgommers
Member

Within a few weeks works too. I just wanted to avoid merging everything at the last minute.

@andreas-h
Contributor

Okay, this should be it. I decided to restrain the interface a bit, i.e. only allow for iopt[0] = 0. Both other choices would require a more complicated interface, and since I don't need those cases, I didn't implement them (yet?).
The class is still useful and working.

@andreas-h
Contributor

PS: For some reason, I cannot get the docs to build on my system, so I cannot guarantee the examples actually work correctly. I copy-pasted from an interactive session, where it worked, though.
It would be good if one of you guys who can build the docs verifies that I didn't screw up on that part and correct what's necessary.

@pv pv commented on the diff Feb 25, 2012
scipy/interpolate/fitpack2.py
+ '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 Member

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

@pv pv commented on the diff Feb 25, 2012
scipy/interpolate/fitpack2.py
+ 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 Member

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

@pv
Member
pv commented Feb 25, 2012

Ok, looks good. Personally, I'd add tests checking that all the provided options work. However, the remaining issues can be addressed by the one doing the merge.

@rgommers
Member
rgommers commented Mar 3, 2012

I rebased/squashed this PR, adapted the docstring to our standard, addressed the remaining issues and fixed the example: https://github.com/rgommers/scipy/tree/pull-117-spgrid. Please have a look.

@rgommers
Member
rgommers commented Mar 3, 2012

I didn;t add any new tests.

@rgommers
Member

OK, I've gone ahead and merged this, because my changes were mostly cosmetic.

Thanks Andreas!

@rgommers rgommers closed this Mar 10, 2012
@charris
Member
charris commented Mar 10, 2012

Now that this has gone in, I agree with Ralf that RectSphereBivariateSpline is the more natural name for the class. Bit late, I know.

@rgommers
Member

We could still change it.

@charris
Member
charris commented Mar 10, 2012

I think that would be the way to go. Sphere without the 'e' just doesn't seem right and might be hard to remember.

@andreas-h
Contributor

OK, I've gone ahead and merged this, because my changes were mostly cosmetic.

Thanks Andreas!


Reply to this email directly or view it on GitHub:
#117 (comment)

Cool, thanks! As for the renaming, go ahead, I'm not dogmatic in that
respect ;)

And I apologize for the git mess I left for you to clean up; I just
yesterday stumbled upon

http://docs.scipy.org/doc/numpy/dev/gitwash/development_workflow.html

and will adhere to it in my next pull request.

@rgommers
Member

Rename done in e7d3e33.

No apologies needed Andreas. Keep them coming:)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment