[SciPy-dev] Interpolating scattered 2D data

[Robert Kern]

Short post:

  http://svn.scipy.org/svn/scipy/branches/newscipy/Lib/sandbox/delaunay/

Long post:

I finally found some code to do 2-D Delaunay triangulations that was

* efficient and robust
* BSD-ish licensed
* coded to behave nicely as a library rather than a standalone program
* compileable with a modern compiler

  http://mapviewer.skynet.ie/voronoi.html

It's a modest C++-ification of Steve Fortune's classic sweepline code
(which previously failed the last two criterion). It needed some
tweaking to expose all of the information I needed for interpolation,
but that was pretty easy.

I've implemented two interpolation algorithms. One is just a linear
interpolation on each triangle using the three node values to determine
a plane. I've also implemented natural neighbor interpolation via
Sibson's method. Currently, only interpolation to a regular, rectangular
grid is supported, but I'll add general interpolation shortly.

Here's an example:

In [5]: import scipy
In [6]: from matplotlib.pylab import *
In [7]: from scipy.sandbox import delaunay
In [8]: x = scipy.random.uniform(-0.5, 1.5, size=1000)
In [9]: y = scipy.random.uniform(-0.5, 1.5, size=1000)
In [10]: z = 2.0*cos(10.0*x)*sin(10.0*y) + sin(10.0*x*y)
In [11]: tri = delaunay.Triangulation(x, y)
In [12]: nni = tri.nn_interpolator(z)
In [13]: x0 = y0 = 0.0
In [14]: x1 = y1 = 1.0
In [15]: f = nni[y0:y1:101j, x0:x1:101j]
In [16]: imshow(f, extent=(x0,x1,y0,y1))
Out[16]: <matplotlib.image.AxesImage instance at 0x31691e8>

You can see a bunch of images of test functions here:

  http://starship.python.net/crew/kernr/nn-images/

Note that the contours aren't being computed in any intelligent way. I
just interpolated to a 101x101 grid and used pylab's contour() function.

-- 
Robert Kern
rkern at ucsd.edu

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter