如果我想要随机非结构化数据的 3D spline/smooth 插值怎么办？
What to do if I want 3D spline/smooth interpolation of random unstructured data?
我受到@James 的启发，想看看如何使用 griddata 和 map_coordinates。在下面的示例中，我展示了 2D 数据，但我对 3D 感兴趣。我注意到 griddata 仅提供 1D 和 2D 的样条曲线，并且仅限于 3D 和更高的线性插值（可能有很好的理由）。但是，map_coordinates 使用更高阶（比分段线性更平滑）插值的 3D 似乎没问题。
我的主要问题： 如果我在 3D 中有随机的、非结构化的数据（我不能在其中使用 map_coordinates），有没有什么方法可以使 3D 更平滑-在 NumPy SciPy 宇宙中或至少在附近进行线性插值？
我的第二个问题：3D 样条曲线在 griddata 中不可用是因为实施起来困难或乏味，还是存在根本性困难？
下面的图像和可怕的 python 显示了我目前对如何使用 griddata 和 map_coordinates 的理解。沿粗黑线进行插值。
结构化数据：
UN结构化数据：
可怕python:
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()
好问题！ （还有漂亮的情节！）
对于非结构化数据，您需要切换回用于非结构化数据的函数。 griddata 是一种选择，但使用三角剖分并在两者之间进行线性插值。这导致三角形边界处有 "hard" 个边。
样条是径向基函数。在 scipy 术语中，您需要 scipy.interpolate.Rbf。我建议在三次样条曲线上使用 function="linear" 或 function="thin_plate"，但也可以使用三次样条曲线。 （与线性或薄板样条相比，三次样条会加剧 "overshooting" 的问题。）
一个警告是径向基函数的这种特殊实现将始终使用数据集中的所有点。这是最准确和最流畅的方法，但随着输入观察点数量的增加，它的扩展性很差。有几种方法可以解决这个问题，但事情会变得更加复杂。我会把它留给另一个问题。
无论如何，这是一个简化的例子。我们将生成随机数据，然后将其插值到规则网格上的点。 （请注意，输入不在规则网格上，插值点也不需要。）
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 且远低于零。这是"overshooting"。
线性样条将完全避免过冲，但您最终会得到 "bullseye" 模式（尽管远不及 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()
我受到@James griddata 和 map_coordinates。在下面的示例中，我展示了 2D 数据，但我对 3D 感兴趣。我注意到 griddata 仅提供 1D 和 2D 的样条曲线，并且仅限于 3D 和更高的线性插值（可能有很好的理由）。但是，map_coordinates 使用更高阶（比分段线性更平滑）插值的 3D 似乎没问题。
我的主要问题： 如果我在 3D 中有随机的、非结构化的数据（我不能在其中使用 map_coordinates），有没有什么方法可以使 3D 更平滑-在 NumPy SciPy 宇宙中或至少在附近进行线性插值？
我的第二个问题：3D 样条曲线在 griddata 中不可用是因为实施起来困难或乏味，还是存在根本性困难？
下面的图像和可怕的 python 显示了我目前对如何使用 griddata 和 map_coordinates 的理解。沿粗黑线进行插值。
结构化数据：
UN结构化数据：
可怕python:
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()
好问题！ （还有漂亮的情节！）
对于非结构化数据，您需要切换回用于非结构化数据的函数。 griddata 是一种选择，但使用三角剖分并在两者之间进行线性插值。这导致三角形边界处有 "hard" 个边。
样条是径向基函数。在 scipy 术语中，您需要 scipy.interpolate.Rbf。我建议在三次样条曲线上使用 function="linear" 或 function="thin_plate"，但也可以使用三次样条曲线。 （与线性或薄板样条相比，三次样条会加剧 "overshooting" 的问题。）
一个警告是径向基函数的这种特殊实现将始终使用数据集中的所有点。这是最准确和最流畅的方法，但随着输入观察点数量的增加，它的扩展性很差。有几种方法可以解决这个问题，但事情会变得更加复杂。我会把它留给另一个问题。
无论如何，这是一个简化的例子。我们将生成随机数据，然后将其插值到规则网格上的点。 （请注意，输入不在规则网格上，插值点也不需要。）
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 且远低于零。这是"overshooting"。
线性样条将完全避免过冲，但您最终会得到 "bullseye" 模式（尽管远不及 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()