Re: [Matplotlib-users] pcolor in action...

[Andrew S. <str...@as...>]

On Apr 14, 2004, at 11:37 AM, Peter Groszkowski wrote:

> Hello Everyone:
>
> Before I reinvent the wheel, I'd thought I'd ask:
>
> I have a set of x,y points in a plane, and to each a corresponding z.  
> They are spread out in circular concentric orbits.
> I need to do a flat surface plot (via pcolor) of these z values. A  
> simple algorithm would be:
>
> #get data for x, y, z.. this would be some fixed points x,y inside
> #a box -4<=x<=4 and -4<=y<=4 with a z value for each.
>
> #Create a grid
> xi=linspace(-4,4,200)
> yi=linspace(-4,4,200)
> [Xi,Yi]=meshgrid(xi,yi)
>
> # This is the key step; find a z value for each grid point.
> Zi=griddata(x,y,z,Xi,Yi)
>
> pcolor(Xi, Yi, Zi, shading='flat')
>
> show()
>
> The key here is this griddata function. It interpolates the surface at  
> the grid points based on the few data points I provide it with  
> (x,y,z). It uses Delaunay triangulation, and is part of the standard  
> MATLAB distribution (more about it here:  
> http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ 
> griddata.html#1007265). My question is whether anyone has come across  
> its implementation in python. I could not find any information in any  
> of the "standard" math packages. How do people do this?

VTK, including the Python API, does Delaunay triangulation, and I've  
used it to do what you're talking about.  IIRC, there are other  
Delaunay triangulation codes with Python bindings around. However, the  
Matlab function is a little more complicated;  they use interpolation  
in addition to a straight Delaunay triangulation, and it often produces  
nice-looking results.  They cite a reference in either the help file or  
the .m file; one can presumably read the reference and re-implement it  
for use in Python.

FYI, here's some relevant Python VTK code:

         profile = vtkPolyData()
         profile.SetPoints(points)

         delny = vtk.vtkDelaunay2D()
         delny.SetInput(profile)
         delny.SetTolerance(0.01)
         delny.BoundingTriangulationOff()

         normals = vtk.vtkPolyDataNormals()
         normals.SetInput(delny.GetOutput())

         lo,hi = zrange
         elevation = vtk.vtkElevationFilter()
         elevation.SetInput(delny.GetOutput())
         elevation.SetLowPoint(0, 0, lo)
         elevation.SetHighPoint(0, 0, hi)
         elevation.SetScalarRange(lo, hi)
         elevation.ReleaseDataFlagOn()

         delMesh = vtk.vtkPolyDataMapper()
         delMesh.SetInput(elevation.GetOutput())
         delMesh.SetLookupTable(lutLand) # doesn't work
         delActor = vtk.vtkActor()
         delActor.SetMapper(delMesh)
         ren1.AddActor( delActor )

Cheers!
Andrew