Search code examples
pythonscipy

How can I do Python 2D Interpolation at a point using SciPy?


I'm trying to interpolate the Z value at specific point on a 2D dataset that was taken experimentally, by providing X and Y and measuring Z. I've seen the examples that generate random functions that are then interpolated across, but I can't seem to get my data in the same grid format that works properly.

Relatively new to interpolation beyond traditional 1D, so would appreciate know where I'm going wrong and theory behind it.

I've attempted to use SciPy Interpolate GridData function but am having issues with floats and integer values as I require the use of decimals in this case.

grid_x = np.array([0, 152, 304, 457, 609, 760])        # Meters from start point
grid_y = np.array(np.linspace(0,0.5,5))                # Normalized temperature    

grid = np.zeros([len(grid_x),len(grid_y)], dtype= 'i,i')    # Make an array of zeroes

for i in range(len(grid_x)):
    for j in range(len(grid_y)):
        grid[i,j] = (grid_x[i], grid_y[j])          # Fill array with premeasured coordinates for X and Y

sunIntensity = np.array(([200, 199, 198, 201, 220],     # Experimental data [Z coordinate]
                   [195, 170, 175, 180, 185],
                   [175, 150, 150, 160, 165],
                   [140, 135, 145, 150, 160],
                   [130, 125, 110, 115, 130],
                   [90, 95, 110, 100, 120]))

InterpPoint = [(100, 0.2)]                  # Point to be interpolated at

InterpSunItensity = griddata(grid, sunIntensity, InterpPoint, method='linear')  # Interpolation 

This provides the following error however:

 TypeError: Cannot cast array data from dtype([('f0', '<i4'), ('f1', '<i4')]) to dtype('float64') according to the rule 'unsafe'

What am I missing here to get this interpolation done?


Solution

  • Your data is completely not in the shape that griddata is asking for. It wants points to be of shape (n, D) which in your case means a 25x2 array, but you've given it a 5x5 array of tuples. And it wants values to be a 1D array of size n (i.e. 25), but you've given it a 2D, 5x5 array.

    Specifying the tuple to be of integers (dtype='i,i') when most of your Y-coordinates are non-integer doesn't really help either but that's beside the point a bit.

    How about this:

    grid_x = np.array([0, 152, 304, 457, 609, 760])
    grid_y = np.linspace(0, 0.5, 5)
    sun_intensity = np.array([[200, 199, 198, 201, 220],
                              [195, 170, 175, 180, 185],
                              [175, 150, 150, 160, 165],
                              [140, 135, 145, 150, 160],
                              [130, 125, 110, 115, 130],
                              [90, 95, 110, 100, 120]])
    interp_point = np.array([[100, 0.2]])
    
    grid_x, grid_y = np.meshgrid(grid_x, grid_y, indexing='ij')
    points = np.vstack((grid_x.flatten(), grid_y.flatten())).T
    val = griddata(points, sun_intensity.flatten(), interp_point, method='linear')
    

    meshgrid takes your two 1-D arrays and turns them into two 2-D arrays, one column-wise and one row-wise. It works like your for loop except the output format is closer to what we want.

    The next line with vstack and flatten takes those two arrays, flattens them, and pairs them up side-by-side, resulting in a 25x2 array with x being the fast-moving coordinate and y being the slow-moving one (I think that's what you wanted — if not, change indexing='ij' to indexing='xy').

    sun_intensity also gets flattened in the call to griddata to make it a 1-D array of 25 elements.