interpolacion - Cómo realizar una interpolación bilineal en Python

scipy python 3

Question

Me gustaría realizar la interpolación blinear usando python.
Ejemplo de punto gps para el que quiero interpolar altura es:

B = 54.4786674627
L = 17.0470721369

utilizando 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

Obtengo 31.964 pero de otro software obtengo 31.961.
¿Mi guión es correcto?
¿Puede usted proporcionar otro enfoque?

Answer 1

Aquí hay una función reutilizable que puedes usar. Incluye doctests 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 el código de prueba agregando:

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

Ejecutar la interpolación en su 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 punto de hacer una función de floor es que generalmente se busca interpolar un valor cuya coordenada se encuentra entre dos coordenadas discretas. Sin embargo, parece que ya tiene los valores reales de coordenadas reales de los puntos más cercanos, lo que hace que sea 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 la esencia de esto es que puede ponderar cada punto en función de cuánto más cerca esté del punto de la interpolación que sus otros vecinos.

Answer 3

Inspirado desde here , se me ocurrió el siguiente fragmento. La API está optimizada para reutilizar muchas veces la misma tabla:

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

Puedes usarlo así:

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 comprobación de errores y se encontrará con problemas si intenta usarla en los límites de los índices (o más allá). Para obtener la versión completa del código, incluida la verificación de errores y la extrapolación opcional, consulte here .

Answer 4

No estoy seguro de si esto ayuda mucho, pero obtengo un valor diferente cuando hago interpolación lineal usando 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 5

También puede consultar la función de interp en matplotlib .