如果我想要随机非结构化数据的3D样条/平滑插值怎么办？

Question

我由此受到启发答案通过@詹姆斯,看看如何griddata和map_coordinates可能被使用.在下面的例子中,我正在显示2D数据,但我的兴趣在于3D.我注意到griddata只提供1D和2D的样条曲线,并且仅限于3D和更高的线性插值(可能有很好的理由).然而,map_coordinates似乎很好用3D使用更高阶(比分段线性更平滑)插值.

我的主要问题是: 如果我在3D中有随机的非结构化数据(我不能使用map_coordinates),那么在NumPy SciPy宇宙中,或者至少在附近,是否有某种方法可以比分段线性插值更平滑？

我的第二个问题是:3D的样条是不可用的,griddata因为实现起来很困难或乏味,或者是否存在根本性的困难？

下面的图片和可怕的蟒蛇显示了我目前对griddata和map_coordinates如何使用或不能使用的理解.沿着粗黑线进行插值.

结构化数据:

联合国结构化数据:

可怕的蟒蛇:

import numpy as np
import matplotlib.pyplot as plt

def g(x, y):
    return np.exp(-((x-1.0)**2 + (y-1.0)**2))

def findit(x, X):  # or could use some 1D interpolation
    fraction = (x - X[0]) / (X[-1]-X[0])
    return fraction * float(X.shape[0]-1)

nth, nr = 12, 11
theta_min, theta_max = 0.2, 1.3
r_min,     r_max     = 0.7, 2.0

theta = np.linspace(theta_min, theta_max, nth)
r     = np.linspace(r_min,     r_max,     nr)

R, TH = np.meshgrid(r, theta)
Xp, Yp  = R*np.cos(TH), R*np.sin(TH)
array = g(Xp, Yp)

x, y = np.linspace(0.0, 2.0, 200), np.linspace(0.0, 2.0, 200)
X, Y = np.meshgrid(x, y)
blob = g(X, Y)

xtest = np.linspace(0.25, 1.75, 40)
ytest = np.zeros_like(xtest) + 0.75

rtest = np.sqrt(xtest**2 + ytest**2)
thetatest = np.arctan2(xtest, ytest)

ir = findit(rtest, r)
it = findit(thetatest, theta)

plt.figure()

plt.subplot(2,1,1)

plt.scatter(100.0*Xp.flatten(), 100.0*Yp.flatten())
plt.plot(100.0*xtest, 100.0*ytest, '-k', linewidth=3)
plt.hold
plt.imshow(blob, origin='lower', cmap='gray')

plt.text(5, 5, "don't use jet!", color='white')


exact = g(xtest, ytest)

import scipy.ndimage.interpolation as spndint
ndint0 = spndint.map_coordinates(array, [it, ir], order=0)
ndint1 = spndint.map_coordinates(array, [it, ir], order=1)
ndint2 = spndint.map_coordinates(array, [it, ir], order=2)

import scipy.interpolate as spint
points = np.vstack((Xp.flatten(), Yp.flatten())).T   # could use np.array(zip(...))
grid_x = xtest
grid_y = np.array([0.75])

g0 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='nearest')
g1 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='linear')
g2 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='cubic')


plt.subplot(4,2,5)

plt.plot(exact, 'or')
#plt.plot(ndint0)
plt.plot(ndint1)
plt.plot(ndint2)
plt.title("map_coordinates")

plt.subplot(4,2,6)

plt.plot(exact, 'or')
#plt.plot(g0)
plt.plot(g1)
plt.plot(g2)
plt.title("griddata")

plt.subplot(4,2,7)

#plt.plot(ndint0 - exact)
plt.plot(ndint1 - exact)
plt.plot(ndint2 - exact)
plt.title("error map_coordinates")

plt.subplot(4,2,8)

#plt.plot(g0 - exact)
plt.plot(g1 - exact)
plt.plot(g2 - exact)
plt.title("error griddata")

plt.show()


seed_points_rand = 2.0 * np.random.random((400, 2))
rr = np.sqrt((seed_points_rand**2).sum(axis=-1))
thth = np.arctan2(seed_points_rand[...,1], seed_points_rand[...,0])
isinside = (rr>r_min) * (rr<r_max) * (thth>theta_min) * (thth<theta_max)
points_rand = seed_points_rand[isinside]

Xprand, Yprand = points_rand.T  # unpack
array_rand = g(Xprand, Yprand)

grid_x = xtest
grid_y = np.array([0.75])

plt.figure()

plt.subplot(2,1,1)

plt.scatter(100.0*Xprand.flatten(), 100.0*Yprand.flatten())
plt.plot(100.0*xtest, 100.0*ytest, '-k', linewidth=3)
plt.hold
plt.imshow(blob, origin='lower', cmap='gray')
plt.text(5, 5, "don't use jet!", color='white')


g0rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='nearest')
g1rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='linear')
g2rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='cubic')

plt.subplot(4,2,6)

plt.plot(exact, 'or')
#plt.plot(g0rand)
plt.plot(g1rand)
plt.plot(g2rand)
plt.title("griddata")


plt.subplot(4,2,8)

#plt.plot(g0rand - exact)
plt.plot(g1rand - exact)
plt.plot(g2rand - exact)
plt.title("error griddata")

plt.show()

Answer 1

好问题!(还有漂亮的地块!)

对于非结构化数据,您需要切换回用于非结构化数据的函数. griddata是一种选择,但使用三角测量和中间线性插值.这导致三角形边界处的"硬"边缘.

样条是径向基函数.用scipy术语,你想要的scipy.interpolate.Rbf.我建议使用function="linear"或function="thin_plate"超过三次样条曲线,但也可以使用立方体.(与线性或薄板样条相比,立方样条会加剧"过冲"的问题.)

需要注意的是,径向基函数的这种特定实现将始终使用数据集中的所有点.这是最准确和最平滑的方法,但随着输入观察点数量的增加,它的扩展性很差.有几种方法,但事情会变得更复杂.我会留下另一个问题.

无论如何,这是一个简化的例子.我们将生成随机数据,然后在常规网格上的点处进行插值.(请注意,输入不在常规网格上,插值点也不需要.)

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt
np.random.seed(1977)

x, y, z = np.random.random((3, 10))

interp = scipy.interpolate.Rbf(x, y, z, function='thin_plate')

yi, xi = np.mgrid[0:1:100j, 0:1:100j]

zi = interp(xi, yi)

plt.plot(x, y, 'ko')
plt.imshow(zi, extent=[0, 1, 1, 0], cmap='gist_earth')
plt.colorbar()

plt.show()

样条类型的选择

我选择"thin_plate"了样条曲线的类型.我们的输入观察点范围从0到1(它们是由它们创建的np.random.random).请注意,我们的插值略高于1且远低于零.这是"过度".

线性样条曲线将完全避免过冲,但你会看到"牛眼"模式(尽管接近IDW方法的严重程度).例如,这是使用线性径向基函数插值的完全相同的数据.请注意,我们的插值不会超过1或低于0:

高阶样条曲线将使数据中的趋势更加连续,但会超出更多.默认值"multiquadric"与薄板样条曲线非常相似,但会使事情更加连续并且过度使用更糟糕:

但是,当你进入更高阶的样条曲线时,如"cubic"(三阶):

和"quintic"(第五顺序)

一旦你的输入数据略微超出你的输入数据,你就可以得到不合理的结果.

无论如何,这是一个简单的例子来比较随机数据的不同径向基函数:

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt
np.random.seed(1977)

x, y, z = np.random.random((3, 10))
yi, xi = np.mgrid[0:1:100j, 0:1:100j]

interp_types = ['multiquadric', 'inverse', 'gaussian', 'linear', 'cubic',
                'quintic', 'thin_plate']

for kind in interp_types:
    interp = scipy.interpolate.Rbf(x, y, z, function=kind)
    zi = interp(xi, yi)

    fig, ax = plt.subplots()
    ax.plot(x, y, 'ko')
    im = ax.imshow(zi, extent=[0, 1, 1, 0], cmap='gist_earth')
    fig.colorbar(im)
    ax.set(title=kind)
    fig.savefig(kind + '.png', dpi=80)

plt.show()