 <p>1) Griddata cannot do this, since it only does volume interpolation. Surface has zero volume, so this use case is out of scope for <code>griddata</code>, hence the errors. Moreover, since floating point has only finite precision, usually <em>no point</em> lies exactly on the surface.</p> <p>Your initial attempt with <code>griddata</code> also cannot work correctly in any case: <code>griddata</code> is triangle-based, and cannot produce correct maps for bilinear surfaces.</p> <p>Do you really need the surface to be bilinear? If yes, there is in general no way to avoid solving a nonlinear minimization problem to find the <code>r</code>, <code>s</code> coordinates, as the coordinate map itself is nonlinear (you multiply <code>r</code> with <code>s</code> in it). Such problems can be solved e.g. with <code>scipy.optimize.leastsq</code>:</p> <pre><code># untested code, but gives the idea def func(c): r, s = c xyz = P0*(1-r)*(1-s) + P1*(r)*(1-s) + P2*(r)*(s) + P3*(1-r)*(s) return xyz - Point_on_Surface res = scipy.optimize.leastsq(func, [0, 0]) r, s = res[0] </code></pre> <p>2) If you are happy with a triangular surface where each triangle is flat, then you can solve the problem with linear projections onto the triangles.</p> <p>So let's assume your surface consists of flat triangles. The problem splits into two subproblems: (A) finding the projection of a 3D point onto a given triangle, and (B) finding the triangle inside which the point projects to.</p> <p>Given corners A, B, C of a triangle (in 3D), any point inside the triangle can be written in the form <code>P = c_1 A + c_2 B + (1 - c_1 - c_2) C</code>. <code>c_k</code> specify the barycentric coordinate system on the triangle. Given a point <code>P</code> in 3D, we can then find the barycentric coordinates of a point closest to it on the triangle via <code>c = np.linalg.lstsq([A-C, B-C], P-C)</code>. So, with <code>V[0], V[1], V[2] = A, B, C</code>,</p> <pre><code>def project_onto_triangle(P, V): c, distance, rank, s = np.linalg.lstsq(np.c_[V[0] - V[2], V[1] - V[2]], P - V[2]) return c, distance Point_on_Surface = 0.25 * P1[0] + 0.25 * P1[1] + 0.5 * P1[2] c, dist = project_onto_triangle(Point_on_Surface, P1[0:3]) # &gt;&gt;&gt; c # array([ 0.25, 0.25]) # barycentric coordinates # &gt;&gt;&gt; dist # array([ 4.07457566e-33]) # the distance of the point from the triangle plane # second triangle c2, dist2 = project_onto_triangle(Point_on_Surface, [P1[0], P1[2], P1[3]]) # &gt;&gt;&gt; c2 # array([ 0.45, 0.75]) # &gt;&gt;&gt; dist2 # array([ 0.05]) </code></pre> <p>If you need to do this repeatedly, note that <code>np.linalg.lstsq(V, P - C) == Q.dot(P - C)</code> if <code>Q = np.linalg.pinv(V)</code> --- so you can precompute the projector matrixes <code>Q</code>. The <code>c</code> coordinates obtained then then tell you the coordinates on the unit square.</p> <p>You will need to go through all the triangles to find the one in which the point lies. A point is on the triangle if both of the following are true:</p> <ul> <li><p><code>dist &lt; epsilon</code>, where <code>epsilon</code> is some small tolerance (point near the triangle plane)</p></li> <li><p><code>0 &lt;= c_1 &lt;= 1</code>, <code>0 &lt;= c_2 &lt;= 1</code> and <code>0 &lt;= 1 - c_1 - c_2 &lt;= 1</code> (point projection inside triangle)</p></li> </ul> <p>In the above example case, we find that the point is in the first triangle. For he second triangle, the distance is small (0.05), but the barycentric coordinates indicate that the projection lies outside the triangle.</p> <p>There are some optimizations one can make in the above code if many points need to be projected at once etc.</p> <p>If this is still too slow, the problem gets harder. <code>griddata</code> uses special trick for finding the correct simplex, but these work only for volume grids. I don't know a fast algorithm in this case, although there might be something suitable in the computational geometry literature.</p>