python - Interpolating a 3D surface known by its corner nodes and coloring it with a colormap

Question

I want to construct a 3D representation of experimental data to track the deformation of a membrane. Experimentally, only the corner nodes are known. However I want to plot the deformaiton of the overall structure and this why I want to interpolate the membrane to enable a nice colormap of it. By searching around, I came almost close to it with the following code:

import numpy
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.interpolate import griddata

x=numpy.array([0, 0, 1, 1])
y=numpy.array([0.5, 0.75, 1, 0.5])
z=numpy.array([0, 0.5, 1,0])

fig = plt.figure()
ax = Axes3D(fig)
verts = [zip(x, y, z)]
PC = Poly3DCollection(verts)
ax.add_collection3d(PC)

xi = numpy.linspace(x.min(),x.max(),20)
yi = numpy.linspace(y.min(),y.max(),20)
zi = griddata((x,y),z, (xi[None,:], yi[:,None]), method='linear')
xig, yig = numpy.meshgrid(xi, -yi)
ax.plot_surface(xig, yig, zi, rstride=1, cstride=1,  linewidth=0,cmap=plt.cm.jet,norm=plt.Normalize(vmax=abs(yi).max(), vmin=-abs(yi).max()))
plt.show()

and get the following plot:

The blue polygon is the surface known by its corner nodes and that I want to colormap. The colormapped surface is my best result so far. However, there are the black polygons near the top of the surface that are troubling me. I think it might be due to the fact that the surface doesn't fit the meshgrid and so the fourth corner is here a Nan.

Is there a workaround to avoid these black triangles or even better a better way of colormapping a surface known only by its corner nodes?

EDIT: Here is the figure with the triangulation solution given in my first comment by using the following command

triang = tri.Triangulation(x, y)
ax.plot_trisurf(x, y, z, triangles=triang.triangles, cmap=cm.jet,norm=plt.Normalize(vmax=abs(yi).max(), vmin=-abs(yi).max()))

Answer 1

Indeed it seems plot_trisurf should be perfect for this task! Additionally, you can make use of tri.UniformTriRefiner to get a Triangulation with smaller triangles:

import numpy
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import tri, cm

x = numpy.array([0, 0, 1, 1])
y = numpy.array([0.5, 0.75, 1, 0.5])
z = numpy.array([0, 0.5, 1, 0])

triang = tri.Triangulation(x, y)
refiner = tri.UniformTriRefiner(triang)
new, new_z = refiner.refine_field(z, subdiv=4)

norm = plt.Normalize(vmax=abs(y).max(), vmin=-abs(y).max())
kwargs = dict(triangles=new.triangles, cmap=cm.jet, norm=norm, linewidth=0.2)

fig = plt.figure()
ax = Axes3D(fig)
pt = ax.plot_trisurf(new.x, new.y, new_z, **kwargs)
plt.show()

Resulting in the following image:

Triangular grid refinement was only recently added to matplotlib so you will need version 1.3 to use it. Though if you would be stuck with version 1.2 you should also be able to use the source from Github directly, if you comment out the line import matplotlib.tri.triinterpolate and all of the refine_field method. Then you need to use the refine_triangulation method and use griddata to interpolate the new corresponding Z-values.

---

Edit: The above code uses cubic interpolation to determine the Z-values for the new triangles, but for linear interpolation you could substitute / add these lines:

interpolator = tri.LinearTriInterpolator(triang, z)
new, new_z = refiner.refine_field(z, interpolator, subdiv=4)

Alternatively, to do the interpolation with scipy.interpolate.griddata:

from scipy.interpolate import griddata

new = refiner.refine_triangulation(subdiv = 4)
new_z = griddata((x,y),z, (new.x, new.y), method='linear')