Interpolate Z values in a 3D surface, starting from an irregular set of points

Question

I have a surface that looks like Figure A, imagine that is top view. The surface has calculated Z value.

Now I need to find all Z values in new points like figure B. How to do this? I tried scipy.interpolate.interp2d but it gives some weird results like this:

I just want to find custom z for custom x and y inside the "figure".

Mimimal code example

func_int = scipy.interpolate.interp2d([point[0] for point in pointsbottom],[point[1] for point in pointsbottom],[point[2] for point in pointsbottom], kind = 'linear')
pointscaption = map(lambda point:(point[0],point[1],func_int(point[0],point[1])),pointscaption)

Where pointsbottom is list of (x,y,z) and pointscaption is list of (x,y,z) but I need to find new z.

Answer 1

Try using griddata instead:

    grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')

The difference is griddata expects regular data as input (hmm..., I think). It's not that you should have a different result but more that you'll be able to catch the problem faster. You can mask the "regular grid" data easily.

My first guess would be that those input coordinates are not what you are expecting them to be (perhaps with a different scale from the function you are calculating) but it's difficult to say without testing.

In any case it seems you want a surface, which by definition is a type of grid data, so it should be fairly easy to spot the problem using this different framework.

EDIT (further considerations about the doubts of the poster):

Let's say you want some object that you want to input some data inside. After doing this you want to be able to estimate any position using that data. For that purpose you can build a class like this:

    import numpy as np

    class Estimation():
        def __init__(self,datax,datay,dataz):
            self.x = datax
            self.y = datay
            self.v = dataz

        def estimate(self,x,y,using='ISD'):
            """
            Estimate point at coordinate x,y based on the input data for this
            class.
            """
            if using == 'ISD':
                return self._isd(x,y)
        
        def _isd(self,x,y):
            d = np.sqrt((x-self.x)**2+(y-self.y)**2)
            if d.min() > 0:
                v = np.sum(self.v*(1/d**2)/np.sum(1/d**2))
                return v
            else:
                return self.v[d.argmin()]

This example is using Inverse Squared Distance method which is very stable for estimation (if you avoid divisions by zero). It won't be fast but I hope it's understandable. From this point on you can estimate any point in 2D space by doing:

    e = Estimation(datax,datay,dataz)
    newZ = e.estimate(30,55) # the 30 and 55 are just example coordinates

If you were to do this to an entire grid:

    datax,datay = np.random.randint(0,100,10),np.random.randint(0,100,10)
    dataz       = datax/datay

    e = Estimation(datax,datay,dataz)

    surf = np.zeros((100,100))
    for i in range(100):
        for j in range(100):
            surf[i,j] = e.estimate(i,j)

You would obtain an image that you can see using, for instance, matplotlib (in which the color represents the height in your surface):

    import matplotlib.pyplot as plt
    plt.imshow(surf.T,origin='lower',interpolation='nearest')
    plt.scatter(datax,datay,c=dataz,s=90)
    plt.show()

The result of this experiment is this:

If you don't want to use ISD (Inverse Squared Distance) just implement a new method on Estimation class. Is this what you are looking for?