当前位置： 文江博客 话题详情

如何在Python中执行双线性插值

发布于 2024-12-23 08:43:11 字数 1014 浏览 1 评论 0 原文

我想使用 python 执行双线性插值。
我想要插入高度的示例 GPS 点是：

B = 54.4786674627
L = 17.0470721369

使用具有已知坐标和高度值的四个相邻点：

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

and here's my primitive attempt: 
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

where z0 and z1
z01  z0  z11

     z
z00  z1   z10

I get 31.964 but from other software I get 31.961.
Is my script correct?
Can You provide another approach?



2022 年编辑：
我要感谢所有在这个问题发表十多年后仍然给出新答案的人。

原文

I would like to perform blinear interpolation using python.
Example gps point for which I want to interpolate height is:

B = 54.4786674627
L = 17.0470721369

using four adjacent points with known coordinates and height values:

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

and here's my primitive attempt:

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

where z0 and z1

z01  z0  z11

     z
z00  z1   z10

I get 31.964 but from other software I get 31.961.
Is my script correct?
Can You provide another approach?

2022 Edit:
I would like to thank everyone who, even more than a decade after publication of this question, gives new answers to it.

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（9）

分开我的手 2024-12-30 08:43:11

这是您可以使用的可重用函数。它包括文档测试和数据验证：

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)

您可以通过添加来运行测试代码：

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

在数据集上运行插值会产生：

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

Here's a reusable function you can use. It includes doctests and data validation:

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)

You can run test code by adding:

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

Running the interpolation on your dataset produces:

>>> 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
回复 原文
弱骨蛰伏 2024-12-30 08:43:11

不确定这是否有很大帮助，但在使用 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])

Not sure if this helps much, but I get a different value when doing linear interpolation using 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])
回复 原文
柳若烟 2024-12-30 08:43:11

受此处的启发，我想出了以下代码片段。该 API 针对多次重复使用同一个表进行了优化：

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

您可以像这样使用它：

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

此版本没有错误检查，如果您尝试在索引边界（或超出索引边界）使用它，您将遇到麻烦。有关代码的完整版本，包括错误检查和可选的外推，请查看此处。

Inspired from here, I came up with the following snippet. The API is optimized for reusing a lot of times the same table:

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

You can use it like this:

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

This version has no error checking and you will run into trouble if you try to use it at the boundaries of the indexes (or beyond). For the full version of the code, including error checking and optional extrapolation, look here.

回复 原文
夏の忆 2024-12-30 08:43:11

基于此公式的 numpy 实现：

在此处输入图像描述

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)

输入：
这里，x_是[x0,x1]的列表，y_是[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])

输出:

array([[31.95912739]])

A numpy implementation of based on this formula:

enter image description here

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)

Input:
Here,x_ is list of [x0,x1] and y_ is list of [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])

Output:

array([[31.95912739]])
回复 原文
拧巴小姐 2024-12-30 08:43:11

也可以参考matplotlib中的interp函数。

You can also refer to the interp function in matplotlib.

回复 原文
旧话新听 2024-12-30 08:43:11

我认为执行 floor 函数的要点是，通常您希望插入一个坐标位于两个离散坐标之间的值。然而，您似乎已经有了最近点的实际坐标值，这使得数学变得简单。

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

从您的示例翻译时可能存在一些错误的逻辑，但其要点是您可以根据每个点比其他邻居更接近插值目标点​​的程度来对每个点进行加权。

I think the point of doing a floor function is that usually you're looking to interpolate a value whose coordinate lies between two discrete coordinates. However you seem to have the actual real coordinate values of the closest points already, which makes it simple math.

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

There might be a bit of erroneous logic in translating from your example, but the gist of it is you can weight each point based on how much closer it is to the interpolation goal point than its other neighbors.

回复 原文
旧伤慢歌 2024-12-30 08:43:11

这与此处定义的解决方案相同，但应用于某些函数并与中可用的interp2d进行比较西皮。我们使用 numba 库使插值函数比 Scipy 实现更快。

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

我们将其设置为一个相当大的插值数组来评估每种方法的性能。

示例函数为

z(x, y)=\sin \left(\frac{\pi x}{2}\right) e^{y / 2}

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

在此处输入图像描述

性能测试

没有 numba 库

bilinear_interpolation 函数的 Python，在本例中，与 numba 版本相同，只是我们进行了更改prange 在 for 循环中使用 python 正常 range，并删除函数装饰器 jit

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

每个循环给出 7.15 秒 ± 107 毫秒（平均值 ± 标准差） 7 次运行，每次 1 次循环）

Python 与 numba numba

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

每个循环给出 2.65 ms ± 70.5 µs（7 次运行的平均值 ± 标准差，每次 100 次循环）

Scipy 实现

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

每个循环给出 6.63 ms ± 145 µs（7 次运行的平均值 ± 标准偏差，每次 100 个循环）

性能测试在“Intel(R) Core(TM) i7-8700K CPU @ 3.70GHz”上执行

This is the same solution as defined here but applied to some function and compared with interp2d available in Scipy. We use numba library to make the interpolation function even faster than Scipy implementation.

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 make it quite a large interpolation array to assess the performance of each method.

The sample function is,

z(x, y)=\sin \left(\frac{\pi x}{2}\right) e^{y / 2}

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

enter image description here

Performance Test

Python without numba library

bilinear_interpolation function, in this case, is the same as numba version except that we change prange with python normal range in the for loop, and remove function decorator jit

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

Gives 7.15 s ± 107 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Python with numba numba

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

Gives 2.65 ms ± 70.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Scipy implementation

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

Gives 6.63 ms ± 145 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Performance tests are performed on 'Intel(R) Core(TM) i7-8700K CPU @ 3.70GHz'

回复 原文
音栖息无 2024-12-30 08:43:11

我建议采用以下解决方案：

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

I suggest the following solution:

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))
回复 原文
鹿童谣 2024-12-30 08:43:11

解决方案此处显示了双线性插值，我应用了他的方法此处。但我对这种方法的适应，简单地计算欧几里德距离到所有四个角的接近度，并将其用作简单的加权平均值，效果更好（我的适应在同一链接中）。

import numpy as np

def func(x, y):
    s1 = 0.2; x1 = 36.5; y1 = 32.5
    s2 = 0.4; x2 = 36.1; y2 = 32.8
    g1 = np.exp( -4 *np.log(2) * ((x-x1)**2+(y-y1)**2) / s1**2)
    g2 = np.exp( -2 *np.log(2) * ((x-x2)**2+(y-y2)**2) / s2**2)    
    return g1 + g2 

D = 20
x = np.linspace(36,37,D)
y = np.linspace(32,33,D)
xx,yy = np.meshgrid(x,y)
zz = func(xx,yy)

def find_corners(xi,yi):
    idx1 = np.searchsorted(x, xi, side="left")
    idx2 = np.searchsorted(y, yi, side="left")
    cs = [(idx2,idx1),(idx2-1,idx1),(idx2,idx1-1),(idx2-1,idx1-1)]
    return cs

def cdist(p1,p2):    
    distances = np.linalg.norm(p1 - p2, axis=1)
    return distances

def cell_interp(x, y, points):
    a = np.array([x,y]).reshape(-1,2)
    b = np.array(points)[:,:2]
    ds = cdist(a,b)
    ds = ds / np.sum(ds)
    ds = 1. - ds
    c = np.array(points)[:,2]
    iz = np.sum(c * ds) / np.sum(ds)
    return iz

def grid_interp(intx,inty):
    cs = find_corners(intx,inty)

    i,j = cs[0][0],cs[0][1]
    i,j = cs[1][0],cs[1][1]
    i,j = cs[2][0],cs[2][1]
    i,j = cs[3][0],cs[3][1]
    
    i0,j0 = cs[0][0],cs[0][1]
    i1,j1 = cs[1][0],cs[1][1]
    i2,j2 = cs[2][0],cs[2][1]
    i3,j3 = cs[3][0],cs[3][1]
    
    introw = [(xx[i0,j0],yy[i0,j0],zz[i0,j0]),
              (xx[i1,j1],yy[i1,j1],zz[i1,j1]),
              (xx[i2,j2],yy[i2,j2],zz[i2,j2]),
              (xx[i3,j3],yy[i3,j3],zz[i3,j3])]
    return cell_interp(intx,inty,introw)

x2 = np.linspace(36.0001,36.9999,D*2)
y2 = np.linspace(32.0001,32.9999,D*2)
xx2,yy2 = np.meshgrid(x2,y2)
zz2 = func(xx2,yy2)

grid_interp_vec = np.vectorize(grid_interp,otypes=[np.float])
zz2_grid = grid_interp_vec(xx2,yy2)
print (np.mean(np.square(zz2-zz2_grid)))

A solution here shows bilinear interpolation, I applied his method here. But my adaption of this method, simply computing a proximity to all four corners from Euclidian distance and using that as simple weighted average works much better (my adaption is in the same link).

import numpy as np

def func(x, y):
    s1 = 0.2; x1 = 36.5; y1 = 32.5
    s2 = 0.4; x2 = 36.1; y2 = 32.8
    g1 = np.exp( -4 *np.log(2) * ((x-x1)**2+(y-y1)**2) / s1**2)
    g2 = np.exp( -2 *np.log(2) * ((x-x2)**2+(y-y2)**2) / s2**2)    
    return g1 + g2 

D = 20
x = np.linspace(36,37,D)
y = np.linspace(32,33,D)
xx,yy = np.meshgrid(x,y)
zz = func(xx,yy)

def find_corners(xi,yi):
    idx1 = np.searchsorted(x, xi, side="left")
    idx2 = np.searchsorted(y, yi, side="left")
    cs = [(idx2,idx1),(idx2-1,idx1),(idx2,idx1-1),(idx2-1,idx1-1)]
    return cs

def cdist(p1,p2):    
    distances = np.linalg.norm(p1 - p2, axis=1)
    return distances

def cell_interp(x, y, points):
    a = np.array([x,y]).reshape(-1,2)
    b = np.array(points)[:,:2]
    ds = cdist(a,b)
    ds = ds / np.sum(ds)
    ds = 1. - ds
    c = np.array(points)[:,2]
    iz = np.sum(c * ds) / np.sum(ds)
    return iz

def grid_interp(intx,inty):
    cs = find_corners(intx,inty)

    i,j = cs[0][0],cs[0][1]
    i,j = cs[1][0],cs[1][1]
    i,j = cs[2][0],cs[2][1]
    i,j = cs[3][0],cs[3][1]
    
    i0,j0 = cs[0][0],cs[0][1]
    i1,j1 = cs[1][0],cs[1][1]
    i2,j2 = cs[2][0],cs[2][1]
    i3,j3 = cs[3][0],cs[3][1]
    
    introw = [(xx[i0,j0],yy[i0,j0],zz[i0,j0]),
              (xx[i1,j1],yy[i1,j1],zz[i1,j1]),
              (xx[i2,j2],yy[i2,j2],zz[i2,j2]),
              (xx[i3,j3],yy[i3,j3],zz[i3,j3])]
    return cell_interp(intx,inty,introw)

x2 = np.linspace(36.0001,36.9999,D*2)
y2 = np.linspace(32.0001,32.9999,D*2)
xx2,yy2 = np.meshgrid(x2,y2)
zz2 = func(xx2,yy2)

grid_interp_vec = np.vectorize(grid_interp,otypes=[np.float])
zz2_grid = grid_interp_vec(xx2,yy2)
print (np.mean(np.square(zz2-zz2_grid)))
回复 原文
~没有更多了~

关于作者

浊酒尽余欢

暂无简介

文章
评论
27 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

佚名

文章 0 评论 0

羁客

文章 0 评论 0

天天爱笑的徐老师

文章 0 评论 0

星

文章 0 评论 0

夏日落

文章 0 评论 0

隐诗

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文