Cómo realizar la interpolación bilineal en Python

Question

Me gustaría realizar la interpolación blinear utilizando Python.
Ejemplo gps punto para el cual quiero para interpolar altura es:

B = 54.4786674627 
L = 17.0470721369 

usando cuatro puntos adyacentes con coordenadas conocidas y valores de altura:

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 

y aquí está mi intento primitivo:

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 

donde Z0 y Z1

z01 z0 z11 

    z 
z00 z1 z10 

consigo 31.964 sino de otro software consigo 31.961.
¿Mi secuencia de comandos es correcta?
¿Puede proporcionar otro enfoque?

Answer 1

Aquí es una función reutilizable que puede utilizar. Incluye prueba unitaria y validación de datos:

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) 

Puede ejecutar código de prueba añadiendo:

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

Ejecución de la interpolación en el conjunto de datos produce:

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

Creo que el objetivo de hacer una función floor es que, por lo general, está buscando interpolar un valor cuya coordenada se encuentra entre dos coordenadas discretas. Sin embargo, parece que ya tiene los valores de coordenadas reales reales de los puntos más cercanos, lo que lo convierte en matemática simple.

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 

Puede haber un poco de lógica errónea en la traducción de su ejemplo, pero el quid de la cuestión es que puede ponderar cada punto en función de la cantidad de más cerca está el punto objetivo de interpolación de sus demás vecinos.

Answer 3

No estoy seguro si esto ayuda mucho, pero me da un valor diferente cuando se hace la interpolación lineal utilizando 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

También puede hacer referencia a la interp function in matplotlib.

Answer 5

Inspirada en here, se me ocurrió el siguiente fragmento.La API está optimizado para la reutilización de una gran cantidad de veces que el mismo cuadro:

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

Usted puede utilizar de esta manera:

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 

Esta versión no tiene ninguna comprobación de errores y que se ejecutará en problemas si se intenta para usarlo en los límites de los índices (o más allá). Para ver la versión completa del código, incluida la comprobación de errores y la extrapolación opcional, consulte here.