Is there any decent Pythonic way to interpolate in a triangular mesh, or would I need to implement that myself? That is to say, given a (X,Y) point we'll call P, and a mesh (vertices at (X,Y) with value Z, forming triangular facets), estimate the value at P. So that means first find the facet that contains the point, then interpolate - ideally a higher order interpolation than just "linearly between the facet's vertices" (i.e., taking into account the neighboring facets)?

I could implement it myself, but if there's already something available in Python....

(I checked scipy.interpolate, but its "meshes" seem to just be regular point grids. This isn't a grid, it's a true 2D mesh; the vertices can be located anywhere.)

I often use matplotlib.tri for this purpose. Here Xv,Yv are the vertices (or nodes) of the triangles, and Zv the values at those nodes:

from matplotlib.tri import Triangulation, LinearTriInterpolator, CubicTriInterpolator#you can add keyword triangles here if you have the triangle array, size [Ntri,3]
triObj = Triangulation(Xv,Yv) #linear interpolation
fz = LinearTriInterpolator(triObj,Zv)
Z = fz(X,Y)
#cubic interpolation
fzc = CubicTriInterpolator(triObj,Zv)
Zc = fz(X,Y)