# Looking for triangulator/interpolator

Discussion in 'Python' started by Grant Edwards, May 27, 2006.

1. ### Grant EdwardsGuest

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?

--
Grant Edwards grante Yow! LOOK!! Sullen
at American teens wearing
visi.com MADRAS shorts and "Flock of
Seagulls" HAIRCUTS!

Grant Edwards, May 27, 2006

2. ### Scott David DanielsGuest

Grant Edwards wrote:
> 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

Scott David Daniels, May 27, 2006

3. ### Grant EdwardsGuest

On 2006-05-26, Scott David Daniels <> wrote:

>> 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.

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

> 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.

--
Grant Edwards grante Yow! Hmmm... an arrogant
at bouquet with a subtle
visi.com suggestion of POLYVINYL
CHLORIDE...

Grant Edwards, May 27, 2006
4. ### Grant EdwardsGuest

On 2006-05-27, Grant Edwards <> wrote:
> On 2006-05-26, Scott David Daniels <> wrote:
>
>>> 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.

>
> 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.

>> 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.

--
Grant Edwards grante Yow! YOW!! Now I
visi.com MICROBIOLOGY and th' new
TAX REFORM laws!!

Grant Edwards, May 27, 2006
5. ### John MachinGuest

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

John Machin, May 27, 2006
6. ### Grant EdwardsGuest

On 2006-05-27, Grant Edwards <> wrote:
> On 2006-05-27, Grant Edwards <> wrote:
>> On 2006-05-26, Scott David Daniels <> wrote:
>>
>>>> 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.

>>
>> 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.

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.

--
Grant Edwards grante Yow! I'm ZIPPY!! Are we
at having FUN yet??
visi.com

Grant Edwards, May 27, 2006
7. ### Travis E. OliphantGuest

Grant Edwards wrote:
> 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

Travis E. Oliphant, May 27, 2006
8. ### Robert KernGuest

Grant Edwards wrote:

> 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
an underlying truth."
-- Umberto Eco

Robert Kern, May 27, 2006
9. ### Grant EdwardsGuest

On 2006-05-27, Grant Edwards <> wrote:

> 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 grante Yow! MMM-MM!! So THIS is
at BIO-NEBULATION!
visi.com

Grant Edwards, May 27, 2006
10. ### Grant EdwardsGuest

On 2006-05-27, Robert Kern <> wrote:
> Grant Edwards wrote:
>
>> 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?

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

--
Grant Edwards grante Yow! Are we on STRIKE yet?
at
visi.com

Grant Edwards, May 27, 2006
11. ### Grant EdwardsGuest

On 2006-05-27, Travis E. Oliphant <> wrote:

>> 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?

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.

--
Grant Edwards grante Yow! I call it a "SARDINE
at ON WHEAT"!
visi.com

Grant Edwards, May 27, 2006
12. ### Travis E. OliphantGuest

Grant Edwards wrote:
> On 2006-05-27, Travis E. Oliphant <> wrote:
>
>>> 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?

>
> 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.
>

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

Travis E. Oliphant, May 27, 2006
13. ### Grant EdwardsGuest

On 2006-05-27, Travis E. Oliphant <> wrote:

> 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))].

--
Grant Edwards grante Yow! Are the STEWED PRUNES
at still in the HAIR DRYER?
visi.com

Grant Edwards, May 27, 2006
14. ### Travis E. OliphantGuest

Grant Edwards wrote:
> On 2006-05-27, Travis E. Oliphant <> wrote:
>
>> 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.

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

Travis E. Oliphant, May 28, 2006