Ik wil graag blineaire interpolatie uitvoeren met python.
Voorbeeld gps-punt waarvoor ik de hoogte wil interpoleren is:

B = 54.4786674627
L = 17.0470721369

met behulp van vier aangrenzende punten met bekende coördinaten en hoogtewaarden:

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

en hier is mijn primitieve poging:

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

waarbij z0 en z1

z01  z0  z11
     z
z00  z1   z10

Ik krijg 31.964 maar van andere software krijg ik 31.961.
Is mijn script correct?
Kunt u een andere benadering geven?

Antwoord 1, autoriteit 100%

Hier is een herbruikbare functie die u kunt gebruiken. Het omvat doctests en gegevensvalidatie:

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)

U kunt testcode uitvoeren door toe te voegen:

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

Het uitvoeren van de interpolatie op uw dataset levert:

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

Antwoord 2, autoriteit 20%

Ik weet niet zeker of dit veel helpt, maar ik krijg een andere waarde bij lineaire interpolatie met 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])

Antwoord 3, autoriteit 15%

Geïnspireerd door hier, bedacht ik het volgende fragment. De API is geoptimaliseerd voor het veelvuldig hergebruiken van dezelfde tabel:

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

Je kunt het als volgt gebruiken:

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

Deze versie heeft geen foutcontrole en u zult in de problemen komen als u het probeert te gebruiken op de grenzen van de indexen (of verder). Voor de volledige versie van de Code, inclusief foutcontrole en optionele extrapolatie, kijk dan hier .

Antwoord 4, Autoriteit 7%

U kunt ook verwijzen naar de interp-functie in matplotlib .

Antwoord 5, Autoriteit 4%

Ik denk dat het punt van een floor-functie doet, is dat gewoonlijk een waarde op zoek is om een ​​waarde te interpoleren waarvan de coördinaten tussen twee afzonderlijke coördinaten ligt. Je lijkt echter al de daadwerkelijke echte coördinaten van de dichtstbijzijnde punten te hebben, wat het eenvoudige wiskunde maakt.

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

Er is misschien een beetje onjuiste logica in het vertalen vanuit uw voorbeeld, maar de GIST ervan is dat u elk punt kunt afgeweren op basis van hoeveel dichter het is voor het interpolatie-doelpunt dan de andere buren.

Antwoord 6, Autoriteit 4%

Een numpy implementatie van op basis van deze formule:

def bilinear_interpolation(x,y,x_,y_,val):
    a = 1 /((x_[1] - x_[0]) * (y_[1] - y_[0]))
    xx = np.array([[x_[1]-x],[x-x_[0]]],dtype='float32')
    f = np.array(val).reshape(2,2)
    yy = np.array([[y_[1]-y],[y-y_[0]]],dtype='float32')
    b = np.matmul(f,yy)
    return a * np.matmul(xx.T, b)

Invoer:
Hier is x_een lijst van [x0,x1]en y_is een lijst van [y0,y1]

bilinear_interpolation(x=54.4786674627,
                       y=17.0470721369,
                       x_=[54.458333,54.5],
                       y_=[17.041667,17.083333],
                       val=[31.993,31.911,31.945,31.866])

Uitvoer:

array([[31.95912739]])

Antwoord 7, autoriteit 2%

Dit is dezelfde oplossing als hiergedefinieerd, maar toegepast op een bepaalde functie en vergeleken met interp2dverkrijgbaar in Scipy. We gebruiken de numba-bibliotheek om de interpolatiefunctie nog sneller te maken dan de Scipy-implementatie.

import numpy as np
from scipy.interpolate import interp2d
import matplotlib.pyplot as plt
from numba import jit, prange
@jit(nopython=True, fastmath=True, nogil=True, cache=True, parallel=True)
def bilinear_interpolation(x_in, y_in, f_in, x_out, y_out):
    f_out = np.zeros((y_out.size, x_out.size))
    for i in prange(f_out.shape[1]):
        idx = np.searchsorted(x_in, x_out[i])
        x1 = x_in[idx-1]
        x2 = x_in[idx]
        x = x_out[i]
        for j in prange(f_out.shape[0]):
            idy = np.searchsorted(y_in, y_out[j])
            y1 = y_in[idy-1]
            y2 = y_in[idy]
            y = y_out[j]
            f11 = f_in[idy-1, idx-1]
            f21 = f_in[idy-1, idx]
            f12 = f_in[idy, idx-1]
            f22 = f_in[idy, idx]
            f_out[j, i] = ((f11 * (x2 - x) * (y2 - y) +
                            f21 * (x - x1) * (y2 - y) +
                            f12 * (x2 - x) * (y - y1) +
                            f22 * (x - x1) * (y - y1)) /
                           ((x2 - x1) * (y2 - y1)))
    return f_out

We maken het een behoorlijk grote interpolatierey om de prestaties van elke methode te beoordelen.

De voorbeeldfunctie is,

x = np.linspace(0, 4, 13)
y = np.array([0, 2, 3, 3.5, 3.75, 3.875, 3.9375, 4])
X, Y = np.meshgrid(x, y)
Z = np.sin(np.pi*X/2) * np.exp(Y/2)
x2 = np.linspace(0, 4, 1000)
y2 = np.linspace(0, 4, 1000)
Z2 = bilinear_interpolation(x, y, Z, x2, y2)
fun = interp2d(x, y, Z, kind='linear')
Z3 = fun(x2, y2)
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(10, 6))
ax[0].pcolormesh(X, Y, Z, shading='auto')
ax[0].set_title("Original function")
X2, Y2 = np.meshgrid(x2, y2)
ax[1].pcolormesh(X2, Y2, Z2, shading='auto')
ax[1].set_title("bilinear interpolation")
ax[2].pcolormesh(X2, Y2, Z3, shading='auto')
ax[2].set_title("Scipy bilinear function")
plt.show()

Prestatietest

Python zonder numba-bibliotheek

bilinear_interpolation-functie is in dit geval hetzelfde als de numba-versie, behalve dat we prangewijzigen met python normaal rangein de for-lus, en verwijder functiedecorator jit

%timeit bilinear_interpolation(x, y, Z, x2, y2) 

Geeft 7,15 s ± 107 ms per lus (gemiddelde ± standaard ontwikkeling van 7 runs, elk 1 lus)

Python met numba numba

%timeit bilinear_interpolation(x, y, Z, x2, y2) 

Geeft 2,65 ms ± 70,5 µs per lus (gemiddelde ± standaard ontwikkeling van 7 runs, 100 lussen elk)

Scipy-implementatie

%%timeit
f = interp2d(x, y, Z, kind='linear')
Z2 = f(x2, y2)

Geeft 6,63 ms ± 145 µs per lus (gemiddelde ± standaard ontwikkeling van 7 runs, 100 lussen elk)

Prestatietests worden uitgevoerd op ‘Intel(R) Core(TM) i7-8700K CPU @ 3,70GHz’

Antwoord 8

Ik stel de volgende oplossing voor:

def bilinear_interpolation(x, y, z01, z11, z00, z10):
    def linear_interpolation(x, z0, z1):
        return z0 * x + z1 * (1 - x)
    return linear_interpolation(y, linear_interpolation(x, z01, z11),
                                   linear_interpolation(x, z00, z10))