How to perform bilinear interpolation in Python

Question

I would like to perform blinear interpolation using python.
Example gps point for which I want to interpolate height is:

B = 54.4786674627
L = 17.0470721369

using four adjacent points with known coordinates and height values:

n = [(54.5, 17.041667, 31.993), (54.5, 17.083333, 31.911), (54.458333, 17.041667, 31.945), (54.458333, 17.083333, 31.866)]

z01    z11

     z
z00    z10

and here's my primitive attempt:

import math
z00 = n[0][2]
z01 = n[1][2]
z10 = n[2][2]
z11 = n[3][2]
c = 0.016667 #grid spacing
x0 = 56 #latitude of origin of grid
y0 = 13 #longitude of origin of grid
i = math.floor((L-y0)/c)
j = math.floor((B-x0)/c)
t = (B - x0)/c - j
z0 = (1-t)*z00 + t*z10
z1 = (1-t)*z01 + t*z11
s = (L-y0)/c - i
z = (1-s)*z0 + s*z1

where z0 and z1

z01  z0  z11

     z
z00  z1   z10

I get 31.964 but from other software I get 31.961.
Is my script correct?
Can You provide another approach?

2022 Edit:
I would like to thank everyone who, even more than a decade after publication of this question, gives new answers to it.

Solution 1

Here's a reusable function you can use. It includes doctests and data validation:

def bilinear_interpolation(x, y, points):
    '''Interpolate (x,y) from values associated with four points.

    The four points are a list of four triplets:  (x, y, value).
    The four points can be in any order.  They should form a rectangle.

        >>> bilinear_interpolation(12, 5.5,
        ...                        [(10, 4, 100),
        ...                         (20, 4, 200),
        ...                         (10, 6, 150),
        ...                         (20, 6, 300)])
        165.0

    '''
    # See formula at:  http://en.wikipedia.org/wiki/Bilinear_interpolation

    points = sorted(points)               # order points by x, then by y
    (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points

    if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:
        raise ValueError('points do not form a rectangle')
    if not x1 <= x <= x2 or not y1 <= y <= y2:
        raise ValueError('(x, y) not within the rectangle')

    return (q11 * (x2 - x) * (y2 - y) +
            q21 * (x - x1) * (y2 - y) +
            q12 * (x2 - x) * (y - y1) +
            q22 * (x - x1) * (y - y1)
           ) / ((x2 - x1) * (y2 - y1) + 0.0)

You can run test code by adding:

if __name__ == '__main__':
    import doctest
    doctest.testmod()

Running the interpolation on your dataset produces:

>>> n = [(54.5, 17.041667, 31.993),
         (54.5, 17.083333, 31.911),
         (54.458333, 17.041667, 31.945),
         (54.458333, 17.083333, 31.866),
    ]
>>> bilinear_interpolation(54.4786674627, 17.0470721369, n)
31.95798688313631

Solution 2

Not sure if this helps much, but I get a different value when doing linear interpolation using scipy:

>>> import numpy as np
>>> from scipy.interpolate import griddata
>>> n = np.array([(54.5, 17.041667, 31.993),
                  (54.5, 17.083333, 31.911),
                  (54.458333, 17.041667, 31.945),
                  (54.458333, 17.083333, 31.866)])
>>> griddata(n[:,0:2], n[:,2], [(54.4786674627, 17.0470721369)], method='linear')
array([ 31.95817681])

Solution 3

Inspired from here, I came up with the following snippet. The API is optimized for reusing a lot of times the same table:

from bisect import bisect_left

class BilinearInterpolation(object):
    """ Bilinear interpolation. """
    def __init__(self, x_index, y_index, values):
        self.x_index = x_index
        self.y_index = y_index
        self.values = values

    def __call__(self, x, y):
        # local lookups
        x_index, y_index, values = self.x_index, self.y_index, self.values

        i = bisect_left(x_index, x) - 1
        j = bisect_left(y_index, y) - 1

        x1, x2 = x_index[i:i + 2]
        y1, y2 = y_index[j:j + 2]
        z11, z12 = values[j][i:i + 2]
        z21, z22 = values[j + 1][i:i + 2]

        return (z11 * (x2 - x) * (y2 - y) +
                z21 * (x - x1) * (y2 - y) +
                z12 * (x2 - x) * (y - y1) +
                z22 * (x - x1) * (y - y1)) / ((x2 - x1) * (y2 - y1))

You can use it like this:

table = BilinearInterpolation(
    x_index=(54.458333, 54.5), 
    y_index=(17.041667, 17.083333), 
    values=((31.945, 31.866), (31.993, 31.911))
)

print(table(54.4786674627, 17.0470721369))
# 31.957986883136307

This version has no error checking and you will run into trouble if you try to use it at the boundaries of the indexes (or beyond). For the full version of the code, including error checking and optional extrapolation, look here.

Solution 4

You can also refer to the interp function in matplotlib.

Solution 5

A numpy implementation of based on this formula:

def bilinear_interpolation(x,y,x_,y_,val):

    a = 1 /((x_[1] - x_[0]) * (y_[1] - y_[0]))
    xx = np.array([[x_[1]-x],[x-x_[0]]],dtype='float32')
    f = np.array(val).reshape(2,2)
    yy = np.array([[y_[1]-y],[y-y_[0]]],dtype='float32')
    b = np.matmul(f,yy)

    return a * np.matmul(xx.T, b)

Input: Here,x_ is list of [x0,x1] and y_ is list of [y0,y1]

bilinear_interpolation(x=54.4786674627,
                       y=17.0470721369,
                       x_=[54.458333,54.5],
                       y_=[17.041667,17.083333],
                       val=[31.993,31.911,31.945,31.866])

Output:

array([[31.95912739]])

View more solutions

Comments

[Ben]

You've got rounding errors and you're rounding??? What happens if you remove floor?

[machine yearning]

What are L and B? The coordinates of the point at which you'd like to interpolate?

[daikini]

@machine yearning that's right

[chris]

One note - latitude and longitude are not planar coordinates, so this result won't get you what you want if you're dealing with large distances.

[ovgolovin]

Aren't you forgetting to divide left, right and z by dy1+dy2, dy1+dy2 and dx1+dx2 respectfully?

[machine yearning]

I'm not sure why you'd do that. dx1, dx2, dy1, and dy2 are all normalized to supplementary values between 0 and 1 (so dy1+dy2 always equals to 1) since dx is the total distance between the left neighbor and the right neighbor, and similarly for dy.

[daikini]

@machine yearning I'm not sure if it is clear that the goal is to interpolate height value for given point what is about 31 meters according to heights of adjacent points 31.993, 31.911, 31.945, 31.866.

[daikini]

@machine yearning Thanks for Your answer.

[machine yearning]

@daikini: Lol yeah that's what I was going for. What I was saying was that with bilinear interpolation you can just do linear interpolation along one axis for two pairs of points, and do linear interpolation along the other axis between the two resulting points. I think it makes more sense to normalize everything to [0, 1] than to try to requantize your discrete intervals.

[sastanin]

griddata interpolates linearly in a simplex (triangle) rather than bilinearly in a rectangle; that means it is doing triangulation (Delaunay?) first.

[DavidC.]

@Raymond Hettinger Thank you very much for this answer. Why would not scipy.interpolate.interp2d work in this case? Isn't the interp2d also a bilinear interpolation since it "Interpolates over a 2-D grid" (source: docs.scipy.org/doc/scipy-0.14.0/reference/generated/…) ?

[Sibbs Gambling]

@DavidC. AFAIK, it is bilinear interpolation when you use kind=linear. Empirically, I've also compared the results between this answer and interp2d with kind=linear -- they are exactly the same.

[Nirmal]

Can this be modified to handle missing (NaN) values?

[Khalil Al Hooti]

Yes, it can @Nirmal, but it needs more efforts

[Nirmal]

scipy.interpolate.griddata does the job perfectly, but Numba doesn't support it.