Как выполнить билинейную интерполяцию в Python

Question

Я хотел бы выполнить блиновую интерполяцию с помощью python.
Пример GPS-точка, для которой я хочу интерполировать высоту:

B = 54.4786674627 
L = 17.0470721369 

с использованием четырех смежных точек с известными координатами и значениями высот:

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 

и вот моя примитивная попытка:

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 

где z0 и z1

z01 z0 z11 

    z 
z00 z1 z10 

я получаю 31,964, но от другого программного обеспечения, я получаю 31.961.
Правильно ли мой сценарий?
Можете ли вы предоставить другой подход?

Answer 1

Вот многоразовой функцией, которую вы можете использовать. Она включает в себя doctests и проверки данных:

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) 

Вы можете запустить тестовый код, добавив:

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

Запуск интерполяцию набора данных производит:

>>> 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 

Answer 2

Я думаю, что точка выполнения функции floor состоит в том, что обычно вы хотите интерполировать значение, координата которого лежит между двумя дискретными координатами. Однако у вас, похоже, есть реальные реальные значения координат ближайших точек, что делает его простой математикой.

z00 = n[0][2] 
z01 = n[1][2] 
z10 = n[2][2] 
z11 = n[3][2] 

# Let's assume L is your x-coordinate and B is the Y-coordinate 

dx = n[2][0] - n[0][0] # The x-gap between your sample points 
dy = n[1][1] - n[0][1] # The Y-gap between your sample points 

dx1 = (L - n[0][0])/dx # How close is your point to the left? 
dx2 = 1 - dx1    # How close is your point to the right? 
dy1 = (B - n[0][1])/dy # How close is your point to the bottom? 
dy2 = 1 - dy1    # How close is your point to the top? 

left = (z00 * dy1) + (z01 * dy2) # First interpolate along the y-axis 
right = (z10 * dy1) + (z11 * dy2) 

z = (left * dx1) + (right * dx2) # Then along the x-axis 

Там может быть немного ошибочной логики в переводе с вашего примера, но суть его, что вы можете взвешивать каждую точку на основе того, насколько он ближе к точке интерполяции цели по сравнению с другими своими соседями.

Answer 3

Не уверен, что это очень помогает, но я получаю другое значение при выполнении линейной интерполяции с использованием 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]) 

Answer 4

Вы также можете обратиться к interp function in matplotlib.

Answer 5

Вдохновленный от here, я придумал следующий фрагмент.API оптимизирован для повторного использования много раз ту же таблицу:

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)) 

Вы можете использовать его как это:

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 

Эта версия не имеет проверки ошибок, и вы столкнетесь с проблемами, если вы пытаетесь использовать его на границах индексов (или за его пределами). Для полной версии кода, включая проверку ошибок и дополнительную экстраполяцию, посмотрите here.