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

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

19 голосов
/ 29 декабря 2011

Я хотел бы выполнить блинейную интерполяцию с использованием 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.
Мой сценарий правильный?
Можете ли вы предоставить другой подход?

Ответы [ 5 ]

38 голосов
/ 29 декабря 2011

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

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
6 голосов
/ 29 декабря 2011

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

Вдохновленный здесь , я придумал следующий фрагмент. 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

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

3 голосов
/ 19 января 2012

Вы также можете обратиться к функции interp в matplotlib .

2 голосов
/ 29 декабря 2011

Думаю, смысл выполнения функции 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

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

...