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.

Answer 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