pythonscipysplinespatial-interpolation

Linear interpolation in two dimensions produces results outside of the data range


I am using the interp2d method to produce liner spline functions. With certain point sets the resulting function is producing results that I would not expect. E.g.:

from scipy import interpolate

x = [81, 81, 81, 83, 83, 83]
y = [ 9,  7,  5,  9,  7,  5]
z = [23.75374, 23.75416, 23.75376, 23.75621, 23.75581, 23.75686]
f = interpolate.interp2d(x, y, z, kind='linear')

print (str(f(82, 6)[0]))
print (str(f(82.5, 6.5)[0]))
print (str(f(81.5, 5.5)[0]))

Produces the following output:

8.07860599193
0.240930002164
15.9162159912

I would try to play around with the stiffness parameter of this interpolation method, but apparently it is not available in interp2d.

What is causing these results? How can they be avoided?


Solution

  • How to avoid this

    When interpolating at the points of a rectangular grid, it is better to tell interp2d that this is what's going on, by giving grid-structured data: m values of x coordinate, n values of y coordinate, and z with shape (n, m). Using same points as yours, we get correct results:

    x = [81, 83]
    y = [9, 7, 5]
    z = np.array([[23.75374, 23.75416, 23.75376], [23.75621, 23.75581, 23.75686]]).T
    f = interpolate.interp2d(x, y, z, kind='linear')
    print(f(82, 6)[0], f(82.5, 6.5)[0], f(81.5, 5.5)[0])
    

    outputs 23.7551475 23.755569375 23.754544375

    Why this happens

    The spline construction routine needs a triangular grid to build a piecewise linear function on. It does not know that the values you pass lie on rectangular grid; it just sees a bunch of points, many of them collinear and does not understand how to form triangles from them. So... it adds three more vertices with coordinates around (82.5, 5), (82.5, 7), (82.5, 9), in order to have such triangles. Problem is, since we don't have values at those points, they are just taken to be zero. Of course the interpolant is worthless then.

    One may ask, why didn't the algorithm show some warning that the result is flaky? It did. Passing your points directly to bisplrep with quiet mode disabled

    spl = interpolate.bisplrep(x, y, z, kx=1, ky=1, s=0, quiet=0)
    

    shows

    RuntimeWarning: Warning. The coefficients of the spline have been computed as the minimal norm least-squares solution of a rank deficient system. kx,ky=1,1 nx,ny=5,5 m=6 fp=0.000000 s=0.000000 warnings.warn(RuntimeWarning(_mess))

    Translation: "I didn't have enough data, so I made some up". The method interp2d calls bisplrep internally but suppresses its error messages...