Looking for triangulator/interpolator

[Grant Edwards]

I need to interpolate an irregularly spaced set of sampled
points: Given a set of x,y,z points, I need to interpolate z
values for a much finer x,y grid.

I tried using the scipy sandbox delaunay module, but the
interpolators don't work: the natural neighbor interpolator
produces a surface with "holes" in it: the interpolator returns
NaNs for no reason for certain regions within the convex hull
(the convex hull looks right, and the input Z values in that
region don't look any different that regions that work).

The linear interpolator just segfualts no matter what data I
try with it.

In the past i've used the Delny/libqhull module, but it doesn't
do interpolation, so I had to write the interpolator in Python
(which works but is very slow).

Is there any "off the shelf" module that does interpolation of
irregularly spaced data?

 

[Scott David Daniels]

I need to interpolate an irregularly spaced set of sampled
points: Given a set of x,y,z points, I need to interpolate z
values for a much finer x,y grid.

I tried using the scipy sandbox delaunay module, but the
interpolators don't work: the natural neighbor interpolator
produces a surface with "holes" in it: the interpolator returns
NaNs for no reason for certain regions within the convex hull
(the convex hull looks right, and the input Z values in that
region don't look any different that regions that work).

Sounds like a great opportunity for you to contribute. The easiest
contribution would be to find as small a case as you can that
demonstrates the problem, fixing it and add a test case would be,
obviously, a greater contribution.

--Scott David Daniels
(e-mail address removed)

 

[Grant Edwards]

Sounds like a great opportunity for you to contribute.

I'll try, but it looks like I'm going to be working all weekend
as it is.

 

[Grant Edwards]

I'll try, but it looks like I'm going to be working all weekend
as it is.

OTOH, it looks like I'm screwed either way. My python
interpolator is so hopelessly slow it's useless in practice. It
can only process 4 points per second and I need to process
arrays of 10,000 to 50,000 elements.

 

[John Machin]

Grant> OTOH, it looks like I'm screwed either way. My python
interpolator is so hopelessly slow it's useless in practice. It
can only process 4 points per second and I need to process
arrays of 10,000 to 50,000 elements.

Pardon my utter ignorance of scipy, but are neither psyco nor pyrex any
use?
Cheers,
John

 

[Grant Edwards]

OTOH, it looks like I'm screwed either way. My python
interpolator is so hopelessly slow it's useless in practice. It
can only process 4 points per second and I need to process
arrays of 10,000 to 50,000 elements.

I found another module that claims to do what I want

http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/python/griddata.html

But, no matter what data I pass, I get either all zeros or all
NaNs back. :/

I'm 0 for 3 now.

 

[Travis E. Oliphant]

I need to interpolate an irregularly spaced set of sampled
points: Given a set of x,y,z points, I need to interpolate z
values for a much finer x,y grid.

How many x,y,z points do you have?

Did you try the fitpack function bisplrep in scipy? It can work well as
long as you don't have too many starting points.

Here is an example of how to use it


from numpy import rand, exp, ogrid
from scipy import interpolate
x = 2*rand(20)-1
y = 2*rand(20)-1
z = (x+y)*exp(-6.0*(x*x+y*y))

tck = interpolate.bisplrep(x,y,z,s=0,xb=-1,xe=1,yb=-1,ye=1)

xnew = r_[-1:1:70j]
ynew = r_[-1:1:70j]
znew = interpolate.bisplev(xnew,ynew,tck)


There is a buglet that is fixed in SVN scipy that means you need to
enter xb, xe, yb, and ye manually.



-Travis

 

[Robert Kern]

I found another module that claims to do what I want

http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/python/griddata.html

But, no matter what data I pass, I get either all zeros or all
NaNs back. :/

I'm 0 for 3 now.

I pointed you to

http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data

earlier. Does that not do what you want?

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

 

[Grant Edwards]

I found another module that claims to do what I want

http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/python/griddata.html

But, no matter what data I pass, I get either all zeros or all
NaNs back. :/

Aaarrrggh. After some more sweating and swearing, it looks like
both the griddata module and the scipy.delaunay.nn_interpolator
do work as long as you pass the preferred brand of arrays to
them and specify the mesh/grid using the right scheme.

My problems seem to have been caused by the interaction of a
number of factors:

1) Gnuplot.py seems to like convert _some_ floating-point
arrays to integer values before plotting them -- this only
seems to happen when passing 2D arrays to splot(). That
was breaking some of my data.

2) Converting the arrays to nested lists prevents the rounding
to an integer problem but it apparently transposed the x/y
axis without my noticing. Then the "holes" in some of the
interpolated surfaces showed up. It turns out there were
NaNs in the input data that were causing the holes, but
because of the transposed x/y axis, I was looking in the
wrong place in the data.

3) Attempting to use the griddata module resulted in mixing
array objects from pylab, numpy, numeric, and scipy (some
of which may or may not be the same -- I can't keep track).
Mixing array types seems to have tripped up some
extensions. AFAICT, python code is happy with any of the
array types since they're pretty much the same if you go by
"duck" typing. But C/Fortran extensions only seem to work
with one sort or the other, and some python modules that
wrap those extensions will pass anything that quacks on
down to C/Fortran code, when then gets confused. Maybe.

4) Even when the arrays were OK, there are a couple
incompatible ways to specify a mesh/grid, and I picked the
wrong one in at least one case.

 

[Grant Edwards]

I pointed you to

http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data

earlier. Does that not do what you want?

Yes, after installer newer versions of things and straightening out
some other issues.

 

[Grant Edwards]

How many x,y,z points do you have?

I've got about 700 data points. They're arranged irregularly
along arcs arranged in sort of truncated fan-shape.

Did you try the fitpack function bisplrep in scipy?

Not recently. I spent a month or so wrestling with bisplrep a
while back, but was unable to get stable results (either with
real world data or with the example from the tutorial):

http://www.visi.com/~grante/scipy/interpolate.html
http://www.scipy.net/pipermail/scipy-user/2004-December/003921.html

I switched to a trianguation scheme shortly after that posting,
and haven't tired bisplrep since then.

 

[Travis E. Oliphant]

I've got about 700 data points. They're arranged irregularly
along arcs arranged in sort of truncated fan-shape.


Not recently. I spent a month or so wrestling with bisplrep a
while back, but was unable to get stable results (either with
real world data or with the example from the tutorial):

http://www.visi.com/~grante/scipy/interpolate.html
http://www.scipy.net/pipermail/scipy-user/2004-December/003921.html

I switched to a trianguation scheme shortly after that posting,
and haven't tired bisplrep since then.

Not that you made a bad choice. I do wonder, how much of your
difficulty was with the interface to the underlying fitpack routines.
The interface is pretty involved as there are lots of parameter choices
to the underlying routine. It's very possible the interface is broken
in some strange way.

Given how many people want to do interpolation. I'm surprised nobody
has looked into fixing up the fitpack routines that do it. I suppose
it's possible that the underlying fitpack routines are faulty in 2-d.

Your examples will help in that study.

Pearu made some nice class-based interfaces to fitpack in 2003. They
are in SciPy, but, it appears that people aren't making much use of
them. Did you try them?

They are in fitpack2.py

There are new wrappers for it

 

[Grant Edwards]

Not that you made a bad choice. I do wonder, how much of your
difficulty was with the interface to the underlying fitpack
routines.

I've no idea. I had never done anything with splines before,
so it's quite possible I just wasn't doing things right. I
never got it to work at all for non-gridded data (which is what
I needed). Since it wasn't stable even for gridded data, I
more or less gave up.

The interface is pretty involved as there are lots of
parameter choices to the underlying routine. It's very
possible the interface is broken in some strange way.

I took a look at underlying Fortran code, but that made my
dizzy.

Given how many people want to do interpolation. I'm surprised
nobody has looked into fixing up the fitpack routines that do
it. I suppose it's possible that the underlying fitpack
routines are faulty in 2-d.

Your examples will help in that study.

Pearu made some nice class-based interfaces to fitpack in
2003. They are in SciPy, but, it appears that people aren't
making much use of them. Did you try them?

Not that I remember.

For that particular application, the interpolation had to be
done in real time, so a triangulation scheme turned out to be a
lot faster [with some constraints on the input data to make the
triangle search O(sqrt(N))].

 

[Travis E. Oliphant]

I've no idea. I had never done anything with splines before,
so it's quite possible I just wasn't doing things right. I
never got it to work at all for non-gridded data (which is what
I needed). Since it wasn't stable even for gridded data, I
more or less gave up.


I took a look at underlying Fortran code, but that made my
dizzy.

Just to finish this thread, I should mention I spent some time yesterday
looking into the details of the underlying surfit fucntion used by
bisplrep. It appears to have difficulties when s=0 is used. It seems
to be geared toward smoothing applications where you are trying to fit
"noisy" scattered data to a smooth function. There are many warnings
about choosing s to be too low.

Thus, trying to use bisplrep for interpolation (instead of data
smoothing) is probably going to be difficult with bisplrep.

We still need a good N-d re-gridding algorithm in SciPy.

-Travis