Làm thế nào để thực hiện nội suy song tuyến trong Python

Question

Tôi muốn thực hiện nội suy blinear bằng python.
Ví dụ gps chỉ mà tôi muốn suy chiều cao là:

B = 54.4786674627 
L = 17.0470721369 

sử dụng bốn điểm liền kề với tọa độ biết và giá trị height:

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 

và đây là nỗ lực ban đầu của tôi:

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 

nơi z0 và z1

z01 z0 z11 

    z 
z00 z1 z10 

tôi nhận được 31,964 nhưng từ phần mềm khác tôi nhận được 31,961.
Tập lệnh của tôi có đúng không?
Bạn có thể cung cấp cách tiếp cận khác không?

Answer 1

Dưới đây là một chức năng tái sử dụng, bạn có thể sử dụng. Nó bao gồm doctests và xác nhận dữ liệu:

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) 

Bạn có thể chạy mã kiểm tra bằng cách thêm:

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

Chạy suy trên bộ dữ liệu của bạn sản xuất:

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

Tôi nghĩ rằng điểm làm một hàm floor là thường bạn đang tìm cách nội suy một giá trị có tọa độ nằm giữa hai tọa độ riêng biệt. Tuy nhiên, bạn dường như có các giá trị tọa độ thực tế thực sự của các điểm gần nhất, điều này làm cho nó trở nên đơn giản.

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 

Có thể có một chút logic có sai sót trong việc dịch từ ví dụ của bạn, nhưng các ý chính của nó là bạn có thể cân nặng mỗi điểm dựa trên cách gần gũi hơn nhiều đó là đến điểm mục tiêu suy hơn các nước láng giềng khác của nó.

Answer 3

Không chắc chắn nếu điều này giúp nhiều, nhưng tôi nhận được một giá trị khác nhau khi thực hiện nội suy tuyến tính sử dụng 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

Bạn cũng có thể tham khảo interp function in matplotlib.

Answer 5

Lấy cảm hứng từ here, tôi đã đưa ra đoạn mã sau.Các API được tối ưu hóa cho tái sử dụng rất nhiều lần cùng một bảng:

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

Bạn có thể sử dụng nó như thế này:

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 

Phiên bản này không có người kiểm tra lỗi và bạn sẽ gặp rắc rối nếu bạn cố gắng để sử dụng nó ở ranh giới của các chỉ mục (hoặc hơn thế nữa). Đối với phiên bản đầy đủ của mã, bao gồm kiểm tra lỗi và ngoại suy tùy chọn, hãy xem here.