griddata performance

Armin_Moser · February 20, 2009, 2:11pm

Hi,

I would like to interpolate an array of shape (801,676) to regularily
spaced datapoints using griddata. This interpolation is quick if the
(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this
condition is not fullfilled the delaunay triangulation is extremely
slow, i.e. not useable. Is this a known property of the used
triangulation? The triangulation can be performed with matlab without
any problems.

Armin

_Jeff_Whitaker1 · February 20, 2009, 3:52pm

Armin Moser wrote:

Hi,

I would like to interpolate an array of shape (801,676) to regularily
spaced datapoints using griddata. This interpolation is quick if the
(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this
condition is not fullfilled the delaunay triangulation is extremely
slow, i.e. not useable. Is this a known property of the used
triangulation? The triangulation can be performed with matlab without
any problems.

Armin

Armin: You could try installing the natgrid toolkit and see if that speeds up griddata at all. If not, please post a test script with data and maybe we can figure out what is going on.

-Jeff

···

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory

_Ryan_May1 · February 20, 2009, 5:18pm

If you’re not using meshgrid, how do you compute x,y? A small complete example would be helpful.

Ryan

···

On Fri, Feb 20, 2009 at 8:11 AM, Armin Moser <armin.moser@…878…2495…> wrote:

Hi,

I would like to interpolate an array of shape (801,676) to regularily

spaced datapoints using griddata. This interpolation is quick if the

(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this

condition is not fullfilled the delaunay triangulation is extremely

slow, i.e. not useable. Is this a known property of the used

triangulation? The triangulation can be performed with matlab without

any problems.

–
Ryan May
Graduate Research Assistant
School of Meteorology
University of Oklahoma

Sent from: Norman Oklahoma United States.

Armin_Moser · February 20, 2009, 7:05pm

Jeff Whitaker wrote:

Armin Moser wrote:

Hi,

I would like to interpolate an array of shape (801,676) to regularily
spaced datapoints using griddata. This interpolation is quick if the
(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this
condition is not fullfilled the delaunay triangulation is extremely
slow, i.e. not useable. Is this a known property of the used
triangulation? The triangulation can be performed with matlab without
any problems.

Armin

Armin: You could try installing the natgrid toolkit and see if that
speeds up griddata at all. If not, please post a test script with data
and maybe we can figure out what is going on.

I have already tried natgrid and it didn't improve the situation. As
suggested I append a script demonstrating the problem.

Thanks
Armin

------8<-------------
from numpy import *
from pylab import *
import time

deg2rad = pi/180.0
ai = 0.12*deg2rad
x = linspace(13,40,676)
y = linspace(10,22,801)

x = x*deg2rad
y = y*deg2rad
[x,y] = meshgrid(x,y)
z = (x**2+y**2)

xi = linspace(x.min(),x.max(),x.shape[1])
yi = linspace(y.min(),y.max(),y.shape[0])
tic= time.time()
zi = griddata(x.flatten(),y.flatten(),z.flatten(),xi,yi)
toc = time.time()
print toc-tic

fac = 2*pi/1.2681
nx = fac * (cos(y)*cos(x) - cos(ai))
ny = fac * (cos(y)*sin(x))
nz = fac * (sin(y) + sin(ai))
np = sqrt(nx**2 + ny**2)

z = (np**2+nz**2)*exp(-0.001*nz)

xi = linspace(np.min(),np.max(),x.shape[1])
yi = linspace(nz.min(),nz.max(),y.shape[0])
tic = time.time()
zi = griddata(np.flatten(),nz.flatten(),z.flatten(),xi,yi)
toc = time.time()
print toc-tic

_Jeff_Whitaker1 · February 20, 2009, 8:04pm

Armin Moser wrote:

Jeff Whitaker wrote:


Armin Moser wrote:


Hi,

I would like to interpolate an array of shape (801,676) to regularily
spaced datapoints using griddata. This interpolation is quick if the
(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this
condition is not fullfilled the delaunay triangulation is extremely
slow, i.e. not useable. Is this a known property of the used
triangulation? The triangulation can be performed with matlab without
any problems.

Armin


Armin: You could try installing the natgrid toolkit and see if that
speeds up griddata at all. If not, please post a test script with data
and maybe we can figure out what is going on.


I have already tried natgrid and it didn't improve the situation. As
suggested I append a script demonstrating the problem.

Thanks
Armin

Armin: On my mac, your two benchmarks take 15 and 14 seconds. Do you consider that too slow?

Perhaps this is just a toy example to test griddata, but I assume you realize that you wouldn't normally use griddata to interpolate data on one regular grid to another regular grid. griddata is strictly for interpolating scatter data (not on a regular mesh) to a regular mesh.

-Jeff

···

------8<-------------
from numpy import *
from pylab import *
import time

deg2rad = pi/180.0
ai = 0.12*deg2rad
x = linspace(13,40,676)
y = linspace(10,22,801)

x = x*deg2rad
y = y*deg2rad
[x,y] = meshgrid(x,y)
z = (x**2+y**2)

xi = linspace(x.min(),x.max(),x.shape[1])
yi = linspace(y.min(),y.max(),y.shape[0])
tic= time.time()
zi = griddata(x.flatten(),y.flatten(),z.flatten(),xi,yi)
toc = time.time()
print toc-tic

fac = 2*pi/1.2681
nx = fac * (cos(y)*cos(x) - cos(ai))
ny = fac * (cos(y)*sin(x))
nz = fac * (sin(y) + sin(ai))
np = sqrt(nx**2 + ny**2)

z = (np**2+nz**2)*exp(-0.001*nz)

xi = linspace(np.min(),np.max(),x.shape[1])
yi = linspace(nz.min(),nz.max(),y.shape[0])
tic = time.time()
zi = griddata(np.flatten(),nz.flatten(),z.flatten(),xi,yi)
toc = time.time()
print toc-tic

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory

Eric_Firing1 · February 20, 2009, 9:10pm

Jeff Whitaker wrote:

Armin Moser wrote:

Jeff Whitaker wrote:


Armin Moser wrote:


Hi,

I would like to interpolate an array of shape (801,676) to regularily
spaced datapoints using griddata. This interpolation is quick if the
(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this
condition is not fullfilled the delaunay triangulation is extremely
slow, i.e. not useable. Is this a known property of the used
triangulation? The triangulation can be performed with matlab without
any problems.

Armin


Armin: You could try installing the natgrid toolkit and see if that
speeds up griddata at all. If not, please post a test script with data
and maybe we can figure out what is going on.


I have already tried natgrid and it didn't improve the situation. As
suggested I append a script demonstrating the problem.

Thanks
Armin

Armin: On my mac, your two benchmarks take 15 and 14 seconds. Do you consider that too slow?

Jeff,

I got 10 s and 8 s on my Thinkpad--but that was without natgrid. When I installed natgrid, the benchmark was taking so long I killed the window.

I'm running 32-bit Ubuntu.

Eric

···

Perhaps this is just a toy example to test griddata, but I assume you realize that you wouldn't normally use griddata to interpolate data on one regular grid to another regular grid. griddata is strictly for interpolating scatter data (not on a regular mesh) to a regular mesh.

-Jeff

------8<-------------
from numpy import *
from pylab import *
import time

deg2rad = pi/180.0
ai = 0.12*deg2rad
x = linspace(13,40,676)
y = linspace(10,22,801)

x = x*deg2rad
y = y*deg2rad
[x,y] = meshgrid(x,y)
z = (x**2+y**2)

xi = linspace(x.min(),x.max(),x.shape[1])
yi = linspace(y.min(),y.max(),y.shape[0])
tic= time.time()
zi = griddata(x.flatten(),y.flatten(),z.flatten(),xi,yi)
toc = time.time()
print toc-tic

fac = 2*pi/1.2681
nx = fac * (cos(y)*cos(x) - cos(ai))
ny = fac * (cos(y)*sin(x))
nz = fac * (sin(y) + sin(ai))
np = sqrt(nx**2 + ny**2)

z = (np**2+nz**2)*exp(-0.001*nz)

xi = linspace(np.min(),np.max(),x.shape[1])
yi = linspace(nz.min(),nz.max(),y.shape[0])
tic = time.time()
zi = griddata(np.flatten(),nz.flatten(),z.flatten(),xi,yi)
toc = time.time()
print toc-tic

Armin_Moser · February 20, 2009, 9:22pm

Jeff Whitaker wrote:

Armin Moser wrote:

Jeff Whitaker wrote:

Armin Moser wrote:


Hi,

I would like to interpolate an array of shape (801,676) to regularily
spaced datapoints using griddata. This interpolation is quick if the
(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this
condition is not fullfilled the delaunay triangulation is extremely
slow, i.e. not useable. Is this a known property of the used
triangulation? The triangulation can be performed with matlab without
any problems.

Armin


Armin: You could try installing the natgrid toolkit and see if that
speeds up griddata at all. If not, please post a test script with data
and maybe we can figure out what is going on.


I have already tried natgrid and it didn't improve the situation. As
suggested I append a script demonstrating the problem.

Thanks
Armin

Armin: On my mac, your two benchmarks take 15 and 14 seconds. Do you
consider that too slow?

No definitely not but up to now I always killed the process (after more
than 10 minutes). I tried the script again and it worked... I have to
check the original on monday at work.

Perhaps this is just a toy example to test griddata, but I assume you
realize that you wouldn't normally use griddata to interpolate data on
one regular grid to another regular grid. griddata is strictly for
interpolating scatter data (not on a regular mesh) to a regular mesh.

Yes I do, but the supporting points of the second example are not on a
regular grid.

Thanks a lot
Armin

Armin_Moser · February 20, 2009, 9:25pm

Eric Firing wrote:

Jeff Whitaker wrote:

Armin Moser wrote:

Jeff Whitaker wrote:

Armin Moser wrote:


Hi,

I would like to interpolate an array of shape (801,676) to regularily
spaced datapoints using griddata. This interpolation is quick if the
(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this
condition is not fullfilled the delaunay triangulation is extremely
slow, i.e. not useable. Is this a known property of the used
triangulation? The triangulation can be performed with matlab without
any problems.

Armin


Armin: You could try installing the natgrid toolkit and see if that
speeds up griddata at all. If not, please post a test script with data
and maybe we can figure out what is going on.


I have already tried natgrid and it didn't improve the situation. As
suggested I append a script demonstrating the problem.

Thanks
Armin

Armin: On my mac, your two benchmarks take 15 and 14 seconds. Do you
consider that too slow?

I got 10 s and 8 s on my Thinkpad--but that was without natgrid. When I
installed natgrid, the benchmark was taking so long I killed the window.

Thats the behaviour I observed (on Windows and Linux) with and without
natgrid. I'm very puzzled at the moment.

Thanks a lot
Armin

Eric_Firing1 · February 20, 2009, 10:07pm

Jeff Whitaker wrote:

Armin Moser wrote:

Jeff Whitaker wrote:


Armin Moser wrote:


Hi,

I would like to interpolate an array of shape (801,676) to regularily
spaced datapoints using griddata. This interpolation is quick if the
(x,y) supporting points are computed as X,Y = meshgrid(x,y). If this
condition is not fullfilled the delaunay triangulation is extremely
slow, i.e. not useable. Is this a known property of the used
triangulation? The triangulation can be performed with matlab without
any problems.

Armin


Armin: You could try installing the natgrid toolkit and see if that
speeds up griddata at all. If not, please post a test script with data
and maybe we can figure out what is going on.


I have already tried natgrid and it didn't improve the situation. As
suggested I append a script demonstrating the problem.

Reducing the original grid from 676x801 to 100x120, I get benchmarks of about 6 seconds with natgrid and 0.15 s with Robert's delaunay. This seems quite repeatable.

I also tried randomizing x and y in the first benchmark with natgrid, and it made only a slight difference.

Eric

···

Thanks
Armin

Armin: On my mac, your two benchmarks take 15 and 14 seconds. Do you consider that too slow?

Perhaps this is just a toy example to test griddata, but I assume you realize that you wouldn't normally use griddata to interpolate data on one regular grid to another regular grid. griddata is strictly for interpolating scatter data (not on a regular mesh) to a regular mesh.

-Jeff

------8<-------------
from numpy import *
from pylab import *
import time

deg2rad = pi/180.0
ai = 0.12*deg2rad
x = linspace(13,40,676)
y = linspace(10,22,801)

x = x*deg2rad
y = y*deg2rad
[x,y] = meshgrid(x,y)
z = (x**2+y**2)

xi = linspace(x.min(),x.max(),x.shape[1])
yi = linspace(y.min(),y.max(),y.shape[0])
tic= time.time()
zi = griddata(x.flatten(),y.flatten(),z.flatten(),xi,yi)
toc = time.time()
print toc-tic

fac = 2*pi/1.2681
nx = fac * (cos(y)*cos(x) - cos(ai))
ny = fac * (cos(y)*sin(x))
nz = fac * (sin(y) + sin(ai))
np = sqrt(nx**2 + ny**2)

z = (np**2+nz**2)*exp(-0.001*nz)

xi = linspace(np.min(),np.max(),x.shape[1])
yi = linspace(nz.min(),nz.max(),y.shape[0])
tic = time.time()
zi = griddata(np.flatten(),nz.flatten(),z.flatten(),xi,yi)
toc = time.time()
print toc-tic