Search code examples
pythonnumpyscipyinterpolation

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


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

enter image description here

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: enter image description here

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.


Solution

  • 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:

    Estimated Surface

    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?