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 .