Nice interpolation

Halldor_Bjornsson · March 30, 2006, 11:33am

First, I want to thank John Hunter and Jouni Seppanen for their very useful help with my xticklabel problems. I have little experience with python, and am a total newbie wrt. matplotlib (- which is very impressive...).

Few days ago I was looking for a method to interpolate 12hourly data, and was not happy with the python alternatives I found. I wanted something in native python, not wrappers to C or Fortran programs
(in this case speed is not a big issue). It so happened that few weeks ago a coworker gave me a faint photocopy of a photocopy from an article published 26 years ago with an interpolation method I had not seen before.

The main point of the method was that it did not generate spline like excursions. The claim was that the method was "consistantly well behaved".

So, I wrote it up in Python (well, I translated my Matlab code). It seems to work as the author claimed. The interpolation method needs x,y and y' (the slope of y(x)) as input, but the article also provided a simple algorithm to estimate the slope. I am including these functions in an attachment. To use them is very simple:

Example where y' is known
  x = seq(0,2*pi,10); y = sin(x); yp = cos(x)
  xi = seq(0,2*pi,40);
  yi = StinemanInterp(x,y,yp,xi);
  plot(x,y,'o',xi,yi)

Example where y' is calculated using slopes()

  x = seq(0,2*pi,10); y = sin(x);
         yp = slopes(x,y)
  xi = seq(0,2*pi,40);
  yi = StinemanInterp(x,y,yp,xi);
  plot(x,y,'o',xi,yi)

I make no claim that the method is fast, nor that my coding is optimal. - I am sure it can be sped up, but it is fast enough for my uses.
Nor do I claim that this is the best interpolation algorithm, - I have no interest in taking part in a religous debate on the merits of different interpolation methods.

For me this is simply a consise method that works resonably well. The slope y' can be calculated using other methods than slopes(),- using divided differences to get the derivative of a quadratic through any three points is straightforward (but not implemented here). Also if yp is simply the linear slope the method yields linear interpolation.

As I said, I have not seen this method before, but it may well be that it is known under a different guise....My coworker got it from a Swiss glaciologist....

Sincerely,
Halldór

niceinterp.py (4.38 KB)

···

--
Halldor Bjornsson
Bustadavegur 9 IS150 Reykjavik Iceland.
e: halldor()vedur.is tel:+354-522600

Norbert_Nemec1 · March 30, 2006, 9:31pm

Nice work, Halldor!

I've spent a bit of time on data interpolation recently, but this
Stineman interpolation method beats everything I came up with in quality
and simplicity.

I took the freedom of going over your code and putting in all the
experience I gathered before working with Python on data interpolation
and related issues. See attached what I came up with. Compared with your
version this is heavily modified:

* The core change was to get rid of the explicit loop in the
interpolation routine. The method now fully exploits the power of numpy.

* The interface of the interpolation routine is changed. Intention was
to move yp to the back and make it optional. If it is not given, the
"slopes" routine is called automatically.

* The "slopes" routine is changed in its core. Instead of a circular
interpolation (which is problematic if the aspect ratio between the x-
and the y-axis is not known) it now uses a parabola interpolation to
estimate the slope at inner points. For most curves where the original
gave reasonable results, this new version should give something similar.
It should, however, also cover many cases where the original "slopes"
routine would have produced garbage.

* The slopes at the endpoints are now also extrapolated in a much
simpler manner. The intent of the original version was not clear to me,
but it would definitely caused problems in several corner cases.

If you could privately send my a scan of the original paper, I would be
very grateful.

I believe, it would be quite some gain for the matplotlib library to
have this algorithm incorporated. As far as I can see, it should indeed
be very robust for producing nice looking interpolations. One should,
however, warn against its use in scientific context: Interpolation is
always a way of making data look nicer than it actually was measured...

Greetings,
Norbert

stinemaninterp.py (5.01 KB)