线性分段数据的Numpy scipy 2D插值

8e2ybdfx  于 12个月前  发布在  其他
关注(0)|答案(3)|浏览(143)

我有几点:

points = np.array([[0, 105],[5000, 105],[0, 135],[5000, 135],[0, 165],[5000, 165]])

字符串
和价值观

values = np.array([[300, 380, 300, 390, 300, 400]]).transpose()


我尝试插值的输入

xi = np.array([[2500, 105],[2500, 125],[2500, 135],[2500, 150],[2500, 165]])


具有双线性插值的预期结果(参考:https://en.wikipedia.org/wiki/Bilinear_interpolation

[340, 343.3333, 345, 347.5, 350]


我的工作为第二个例子使用双线性插值

x1=2500, y1=105 giving z1=340
x2=2500, y2=135 giving z2=345
Hence for x3=2500, y3=125 gives z3=343.3333


但随着

gd = griddata(points, values, xi, method='linear', rescale=True)


我正在得到结果

[340, 345, 345, 345, 350]


我一定是错过了一些简单的东西在这里,但没有得到任何尝试多种不同的方法。

g0czyy6m

g0czyy6m1#

如果您正确地提供了数据,则可以使用scipy.interpolate.interpn完成此操作。(对于2D情况)。然后以与网格对应的格式定义值,这是2D情况下np.meshgrid(x, y, indexing="ij")的结果。请确保严格按升序或降序提供x和y,否则interpn将抛出错误。

import numpy as np
from scipy.interpolate import interpn

x = np.array([0, 5000])
y = np.array([105, 135, 165])

# data corresponds to (x, y) given by np.meshgrid(x, y, indexing="ij")
values = np.array([[300, 300, 300],
                   [380, 390, 400]])

xi = np.array([[2500, 105],
               [2500, 125],
               [2500, 135],
               [2500, 150],
               [2500, 165]])

interpolated = interpn((x, y), values, xi) # array([340., 343.33333333, 345., 347.5, 350.])

字符串
这里它是作为一个函数编写的,尽管它没有一般使用所需的所有功能,即检查输入,确保正确的数据库等。

import numpy as np
from scipy.interpolate import interpn

# moved one of the points and values to show it works with both unsorted
# x and y values
points = np.array([[0, 105],[5000, 105],[5000, 135],[0, 165],[5000, 165],[0, 135]])
values = np.array([[300, 380, 390, 300, 400, 300]]).T

xi = np.array([[2500, 105],[2500, 125],[2500, 135],[2500, 150],[2500, 165]])

def bilinear_interpolation(points, values, xi):
    p = points.copy()
    v = values.copy()

    # sorting the points to ensure ascending order
    ysorted = np.argsort(p[:,1])
    p = p[ysorted]
    v = v[ysorted]
    xsorted = np.argsort(p[:,0])
    p = p[xsorted]
    v = v[xsorted]
    
    x = np.unique(p[:,0])
    y = np.unique(p[:,1])
    v = v.reshape(x.size, y.size)

    return interpn((x, y), v, xi)

interpolated = bilinear_interpolation(points, values, xi)

ukdjmx9f

ukdjmx9f2#

我运行代码,得到和你一样的结果。

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator
import numpy as np
from scipy.interpolate import griddata

points = np.array([[0, 105], [5000, 105], [0, 135], [5000, 135], [0, 165], [5000, 165]])
values = np.array([[300, 380, 300, 390, 300, 400]]).transpose()
xi = np.array([[2500, 105],[2500, 125],[2500, 135],[2500, 150],[2500, 165]])

gd = griddata(points, values, xi, method='linear', rescale=True)
print(gd)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the surface
ax.plot_trisurf(points[:,0], points[:,1], values.flatten(), cmap=cm.get_cmap('Blues'), linewidth=0.2)

# Customize the axes
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.xaxis.set_major_locator(LinearLocator(3))
ax.yaxis.set_major_locator(LinearLocator(3))
ax.zaxis.set_major_locator(LinearLocator(5))

plt.show()

字符串
结果是:

[[340.]
 [345.]
 [345.]
 [345.]
 [350.]]


如果我把linear改为cubic,我得到:

[[340.46073832]
 [343.60340822]
 [345.00000255]
 [347.06944598]
 [349.53926443]]


这是更接近你所期望的.所以它看起来像一个插值方法的问题。
它看起来就像你在这个曲面上插值:


的数据
如果你能展示你是如何达到预期结果的,那么这将是有帮助的。

up9lanfz

up9lanfz3#

这是一个很长的问题:你可以写一个函数来完成同样的任务:

import numpy as np
def interpolate(p, points, values):
    xpoint, ypoint = np.unique(points[:,0]), np.unique(points[:,1])
    z = values.reshape(ypoint.size, xpoint.size).T
    x, y = p[:,0], p[:,1]
    xind = np.searchsorted(xpoint, x)
    yind = np.searchsorted(ypoint, y)
    xind_min = np.maximum(xind - 1, 0)
    yind_min = np.maximum(yind - 1, 0)
    Q11 = z[xind_min, yind_min]
    Q12 = z[xind_min, yind]
    Q21 = z[xind, yind_min]
    Q22 = z[xind, yind]
    x1, x2 = xpoint[xind_min], xpoint[xind]
    y1, y2 = ypoint[yind_min], ypoint[yind]
    x_diff = np.where(x2 > x1, x2 - x1, 1)
    xy1 = ((x2 - x)*Q11 + (x - x1)*Q21) / x_diff
    xy2 = ((x2 - x)*Q12 + (x - x1)*Q22) / x_diff
    y_diff = np.where(y2 > y1, y2 - y1, 1)
    res = ((y2 - y) * xy1 + (y - y1) * xy2)/y_diff 
    equalX, equalY = x2 == x, y2 == y
    equalXY = equalX & equalY
    res[equalXY] =  z[xind[equalXY], yind[equalXY]]
    res[equalX & ~equalY] = z.mean(1)[xind[equalX & ~equalY]]
    res[~equalX & equalY] = z.mean(0)[yind[~equalX & equalY]]
    return res

points = np.array([[0, 105],[5000, 105],[0, 135],[5000, 135],[0, 165],[5000, 165]])
values = np.array([[300, 380, 300, 390, 300, 400]]).transpose()
xi = np.array([[2500, 105],[2500, 125],[2500, 135],[2500, 150],[2500, 165]])

interpolate(xi, points, values)
array([340.        , 343.33333333, 345.        , 347.5       ,
       350.        ])

字符串

相关问题