[SciPy-dev] Interpolating scattered 2D data
Robert Kern
rkern at ucsd.edu
Mon Oct 31 19:32:49 CST 2005
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
More information about the Scipy-dev
mailing list