多变量数据在规则网格上的插值(RegularGridInterpolator
)#
假设您在规则网格上有N维数据,并且您希望对其进行插值。在这种情况下,RegularGridInterpolator
可能会很有用。支持多种插值策略:最近邻、线性和奇数次张量积样条。
严格来说,此类高效处理在*矩形*网格上给出的数据:具有可能不等间距的点的高维超立方格。每个维度的点数可以不同。
以下示例演示了其用法,并比较了使用每种方法的插值结果。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import RegularGridInterpolator
假设我们想要插值这个2维函数。
>>> def F(u, v):
... return u * np.cos(u * v) + v * np.sin(u * v)
假设我们只知道规则网格上的一些数据。
>>> fit_points = [np.linspace(0, 3, 8), np.linspace(0, 3, 11)]
>>> values = F(*np.meshgrid(*fit_points, indexing='ij'))
创建测试点和真实值以进行评估。
>>> ut, vt = np.meshgrid(np.linspace(0, 3, 80), np.linspace(0, 3, 80), indexing='ij')
>>> true_values = F(ut, vt)
>>> test_points = np.array([ut.ravel(), vt.ravel()]).T
我们可以创建插值器并使用每种方法插值测试点。
>>> interp = RegularGridInterpolator(fit_points, values)
>>> fig, axes = plt.subplots(2, 3, figsize=(10, 6))
>>> axes = axes.ravel()
>>> fig_index = 0
>>> for method in ['linear', 'nearest', 'slinear', 'cubic', 'quintic']:
… im = interp(test_points, method=method).reshape(80, 80) … axes[fig_index].imshow(im) … axes[fig_index].set_title(method) … axes[fig_index].axis(“off”) … fig_index += 1 >>> axes[fig_index].imshow(true_values) >>> axes[fig_index].set_title(“真实值”) >>> fig.tight_layout() >>> fig.show()
如预期所示,更高阶的样条插值最接近真实值,尽管计算成本比 linear 或 nearest 更高。slinear 插值也与 linear 插值相匹配。
如果你的数据导致样条方法产生振铃现象,你可以考虑使用 method=”pchip”,它使用 PCHIP 插值器的张量积,每个维度使用一个 PchipInterpolator
。
如果你更喜欢函数接口而不是显式创建类实例,interpn
便捷函数提供了等效的功能。
具体来说,这两种形式给出相同的结果:
>>> from scipy.interpolate import interpn
>>> rgi = RegularGridInterpolator(fit_points, values)
>>> result_rgi = rgi(test_points)
和
>>> result_interpn = interpn(fit_points, values, test_points)
>>> np.allclose(result_rgi, result_interpn, atol=1e-15)
True
对于数据局限于 N 维空间中的 (N-1) 维子空间,即当其中一个网格轴的长度为 1 时,沿该轴的外推由 fill_value
关键字参数控制:
>>> x = np.array([0, 5, 10])
>>> y = np.array([0])
>>> data = np.array([[0], [5], [10]])
>>> rgi = RegularGridInterpolator((x, y), data,
... bounds_error=False, fill_value=None)
>>> rgi([(2, 0), (2, 1), (2, -1)]) # 沿轴外推值
array([2., 2., 2.])
>>> rgi.fill_value = -101
>>> rgi([(2, 0), (2, 1), (2, -1)])
array([2., -101., -101.])
如果输入数据的维度具有不可通约的单位且相差多个数量级,插值结果可能存在数值伪影。建议在插值前对数据进行缩放。
均匀分布的数据#
如果你处理的是整数坐标下的笛卡尔网格数据,例如重采样图像数据,这些例程可能不是最佳选择。建议使用 scipy.ndimage.map_coordinates
代替。
对于等间距网格上的浮点数据,可以将 map_coordinates
轻松封装成类似于 RegularGridInterpolator
的接口。以下是一个源自 Johanness Buchner 的 ‘regulargrid’ 包 的简单示例:
class CartesianGridInterpolator:
def __init__(self, points, values, method='linear'):
self.limits = np.array([[min(x), max(x)] for x in points])
self.values = np.asarray(values, dtype=float)
self.order = {'linear': 1, 'cubic': 3, 'quintic': 5}[method]
def __call__(self, xi):
"""
`xi` 是一个类数组对象(数组或列表),表示多个点。
每个“点”是一个 ndim 维的类数组对象,表示 ndim 维空间中的一个点的坐标。
"""
# 将 xi 数组转置为 ``map_coordinates`` 约定的形式,
# 即将点的坐标沿二维数组的列排列。
xi = np.asarray(xi).T
# 将数据坐标转换为像素坐标
ns = self.values.shape
coords = [(n-1)*(val - lo) / (hi - lo)
for val, n, (lo, hi) in zip(xi, ns, self.limits)]
# 插值
return map_coordinates(self.values, coords,
order=self.order,
cval=np.nan) # fill_value
这个包装器可以用作(几乎)`RegularGridInterpolator` 的直接替代品:
>>> x, y = np.arange(5), np.arange(6)
>>> xx, yy = np.meshgrid(x, y, indexing='ij')
>>> values = xx**3 + yy**3
>>> rgi = RegularGridInterpolator((x, y), values, method='linear')
>>> rgi([[1.5, 1.5], [3.5, 2.6]])
array([ 9. , 64.9])
>>> cgi = CartesianGridInterpolator((x, y), values, method='linear')
array([ 9. , 64.9])
请注意,上面的示例使用了 map_coordinates
的边界条件。
因此,cubic
和 quintic
插值的结果可能与 RegularGridInterpolator
的结果不同。
有关边界条件和其他附加参数的更多详细信息,请参阅 scipy.ndimage.map_coordinates
文档。
最后,我们注意到这个简化的示例假设输入数据按升序排列。