smooth, generalised 2D linear interpolation in python

Question

I'm trying to find a method of linear interpolation in 2D over a regular grid using python, but each proposed type in scipy seems to have it's disadvantages.

My aim is basically:

• Have smooth linearly interpolated data over a regular grid, or as close as possible
• The original data can be at arbitrary locations
• (optional) Linearly extrapolate to the edges

But all the functions seem to have problems with this:

• Functions like griddata, interp2d, LinearNDInterpolator appear to create triangles and interpolate within them, creating a bunch of hard lines/creases that I can't have.
• Rbf seems at first to do exactly what I want, but when presented with planes that are flat, it generates an interpolation based on some kind of sphere, creating a curved surface.

If Rbf would simply interpolate a flat plane as a flat plane using the linear setting, as would be expected, it'd be perfect.

Are there any ideas on how to achieve this, or if there's another function that does what I'm after? I've attached a basic example below.

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt

#create some objects to store data
x=np.empty((2,2))
y=np.empty((2,2))
f_shape=(100,100,100)

#generate coordinates
x[0, :] =  0
x[-1, :] = f_shape[0]
y[:, 0] = 0
y[:, -1] = f_shape[1]

#--------constant height----------
z=np.full((2,2),50)

#create interpolation function and interpolate across grid
interp=interpolate.Rbf(x,y,z,function='linear')
grid=np.mgrid[0:f_shape[0],0:f_shape[1]]
result=interp(grid[0,:,:],grid[1,:,:])

plt.imshow(result) #incorrect curved surface from constant height!!!

#--------random heights-----------
z=np.random.uniform(25,75,(2,2))

#create interpolation function and interpolate across grid
interp=interpolate.Rbf(x,y,z,function='linear')
grid=np.mgrid[0:f_shape[0],0:f_shape[1]]
result=interp(grid[0,:,:],grid[1,:,:])

plt.imshow(result) #Apparently nice smooth linear-ish interpolation

Incorrect curved surface from constant height:

Apparently nice smooth linear-ish interpolation:

Answer 1

I've managed to write a function that suits my purpose. It interpolates (fills in) a plane from a grid of coordinates by interpolating along the grid lines, then interpolating the plane in the x and y directions, and taking the average of the two.

It should be possible to speed this up a bit by reshaping the coordinates into a 1D vector, interpolating the plane in one go, then reshaping back into 2D. However this code is certainly fast enough already for reasonable plane sizes.

Seems to work ok if the coordinates are outside of the plane too. Extrapolation also works if the grid is approximately regular. It'll extrapolate regardless but you'll start to see some sharp creases away from the edge as the grid irregularity increases.

Here's the code. An example is provided in the docstring.

def interlin2d(x,y,z,fsize):
    """
    Linear 2D interpolation of a plane from arbitrary gridded points.

    :param x: 2D array of x coordinates
    :param y: 2D array of y coordinates
    :param z: 2D array of z coordinates
    :param fsize: Tuple of x and y dimensions of plane to be interpolated.
    :return: 2D array with interpolated plane.

    This function works by interpolating lines along the grid point in both dimensions,
    then interpolating the plane area in both the x and y directions, and taking the 
    average of the two. Result looks like a series of approximately curvilinear quadrilaterals.

    Note, the structure of the x,y,z coordinate arrays are such that the index of the coordinates
    indicates the relative physical position of the point with respect to the plane to be interpoalted.

    Plane is allowed to be a subset of the range of grid coordinates provided. 
    Extrapolation is accounted for, however sharp creases will start to appear
    in the extrapolated region as the grid of coordinates becomes increasingly irregular.

    Scipy's interpolation function is used for the grid lines as it allows for proper linear extrapolation. 
    However Numpy's interpolation function is used for the plane itself as it is robust against gridlines
    that overlap (divide by zero distance).


    Example:
    #set up number of grid lines and size of field to interpolate
    nlines=[3,3]
    fsize=(100,100,100)

    #initialize the coordinate arrays
    x=np.empty((nlines[0],nlines[1]))
    y=np.empty((nlines[0],nlines[1]))
    z=np.random.uniform(0.25*fsize[2],0.75*fsize[2],(nlines[0],nlines[1]))

    #set random ordered locations for the interior points
    spacings=(fsize[0]/(nlines[0]-2),fsize[1]/(nlines[1]-2))
    for k in range(0, nlines[0]):
        for l in range(0, nlines[1]):
            x[k, l] = round(random.uniform(0, 1) * (spacings[0] - 1) + spacings[0] * (k - 1) + 1)
            y[k, l] = round(random.uniform(0, 1) * (spacings[1] - 1) + spacings[1] * (l - 1) + 1)

    #fix the edge points to the edge
    x[0, :] = 0
    x[-1, :] = fsize[1]-1
    y[:, 0] = 0
    y[:, -1] = fsize[0]-1

    field = interlin2d(x,y,z,fsize)
    """

    from scipy.interpolate import interp1d
    import numpy as np

    #number of lines in grid in x and y directions
    nsegx=x.shape[0]
    nsegy=x.shape[1]

    #lines along the grid points to be interpolated, x and y directions
    #0 indicates own axis, 1 is height (z axis)
    intlinesx=np.empty((2,nsegy,fsize[0]))
    intlinesy=np.empty((2,nsegx,fsize[1]))

    #account for the first and last points being fixed to the edges
    intlinesx[0,0,:]=0
    intlinesx[0,-1,:]=fsize[1]-1
    intlinesy[0,0,:]=0
    intlinesy[0,-1,:]=fsize[0]-1

    #temp fields for interpolation in x and y directions
    tempx=np.empty((fsize[0],fsize[1]))
    tempy=np.empty((fsize[0],fsize[1]))

    #interpolate grid lines in the x direction
    for k in range(nsegy):
        interp = interp1d(x[:,k], y[:,k], kind='linear', copy=False, fill_value='extrapolate')
        intlinesx[0,k,:] = np.round(interp(range(fsize[0])))
        interp = interp1d(x[:, k], z[:, k], kind='linear', copy=False, fill_value='extrapolate')
        intlinesx[1, k, :] = interp(range(fsize[0]))
    intlinesx[0,:,:].sort(0)

    # interpolate grid lines in the y direction
    for k in range(nsegx):
        interp = interp1d(y[k, :], x[k, :], kind='linear', copy=False, fill_value='extrapolate')
        intlinesy[0, k, :] = np.round(interp(range(fsize[1])))
        interp = interp1d(y[k, :], z[k, :], kind='linear', copy=False, fill_value='extrapolate')
        intlinesy[1, k, :] = interp(range(fsize[1]))
    intlinesy[0,:,:].sort(0)

    #interpolate plane in x direction
    for k in range(fsize[1]):
        tempx[k, :] = np.interp(range(fsize[1]),intlinesx[0,:,k], intlinesx[1,:,k])

    #interpolate plane in y direction
    for k in range(fsize[1]):
        tempy[:, k] = np.interp(range(fsize[0]), intlinesy[0, :, k], intlinesy[1, :, k])

    return (tempx+tempy)/2

Example of interpolation based on 9 points (shown as red dots)

Answer 2

Scipy's griddata works just fine:

import numpy as np
from scipy import interpolate as intp
import matplotlib.pyplot as plt
%matplotlib inline

grid_size = G = 100
height = H = 50

points = np.array([
    (0, 0),
    (G-1, 0),
    (0, G-1),
    (G-1, G-1)
], dtype=np.float32)

gy, gx = np.mgrid[:G, :G]
result = intp.griddata(points, np.full(points.shape[0], H), (gy, gx))

And the plot:

plt.imshow(result, interpolation='none')
plt.colorbar()

And just to be sure:

>>> np.allclose(result, 50)
True