Search code examples
python3dinterpolationspline

How to create a 3D surface from given coordinates with python?


I have some discrete coordinates with their heights and I need to create a smooth surface which I will continue to use. I need to have heights for all the coordinates in that surface. I was considering using a 3D spline but I am having trouble applying the methods suggested in comparable kind of questions. I am not an experienced programmer so I would probably find any suggestions helpful. Please forgive me if I may asked something very similar to what was already discussed before. I am using Python 3.6.

EDIT

I made a tiny list (similar to what I have).

Z=[]
Z.append([20.2, 20.1, 35])
Z.append([20.1, 24.5, 36])
Z.append([21.0, 23.2, 33])
Z.append([22.3, 20.0, 34])
Z.append([22.3, 19.5, 28])
Z.append([20.1, 19.5, 27])
Z.append([20.1, 24.6, 31])
Z.append([22.3, 24.6, 32])

The first number here represents longitude, the second - latitude and the third - altitude. I need to make a "surface" (not necessarily to plot it) which would contain the information about all the coordinates in between this rectangle.

I believe this could sound unclear but maybie the problem is that I do not know the exact structure or the type of data I am looking for.


Solution

  • One possibility is to use method griddata from scipy.

    Here is small example how to use neareast neighbour interpolation method with your data:

    import numpy as np
    from scipy.interpolate import griddata
    # --------------------
    Z=[]
    Z.append([20.2, 20.1, 35])
    Z.append([20.1, 24.5, 36])
    Z.append([21.0, 23.2, 33])
    Z.append([22.3, 20.0, 34])
    Z.append([22.3, 19.5, 28])
    Z.append([20.1, 19.5, 27])
    Z.append([20.1, 24.6, 31])
    Z.append([22.3, 24.6, 32])
    # ---------------------------
    xin=np.array(Z)[:,0];
    yin=np.array(Z)[:,1];
    zin=np.array(Z)[:,2];
    # ----------------------------
    xout=np.linspace(20.,23.,10);
    yout=np.linspace(19.,25.,10);
    xout,yout = np.meshgrid(xout,yout);
    # ----------------------------
    zout=griddata((xin,yin),zin,(xout,yout),'nearest');
    # -----------------------------
    from pylab import pcolormesh,show
    pcolormesh(xout,yout,zout);show();