在多维规则网格上比较插值并处理批次输出

本文译写自 SciPy 1.18.0 官方教程《Multivariate data interpolation on a regular grid》，原文由 SciPy 文档贡献者维护，页面没有单独个人署名或明确首次发表日期。以下示例需要 NumPy、SciPy 和 Matplotlib；代码与标注为“原文输出”的结果经过静态核对，本次没有执行计算或绘图。

规则网格不要求每条轴等间距

如果手中的 N 维数据落在一个规则网格上，可以使用 RegularGridInterpolator 插值。它支持最近邻、线性以及奇数阶张量积样条等方法。更严格地说，它处理的是直线网格（rectilinear grid）：由各坐标轴的一维节点组合而成，轴内相邻节点的距离可以不相等，各维度上的节点数也可以不同。

下面用已知的二维函数生成数据，然后在更密集的网格上比较插值和真值。这样的比较能够说明这个函数上的效果，但不能推出某一种方法对所有数据都最好。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import RegularGridInterpolator

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

两个输入轴都从 0 到 3，一个有 8 个节点，一个有 11 个节点。indexing='ij' 使数组前两维对应这两个轴。接着构造 80×80 的查询网格，并用原函数计算真值：

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("True values")
fig.tight_layout()
fig.show()

原文配图显示：在这个光滑函数上，高阶样条更接近真值，不过比线性或最近邻更费计算；slinear 与 linear 的插值结果相符。如果数据让样条产生振铃，可以考虑 method="pchip"：它使用各维度上的 PchipInterpolator，组合成 PCHIP 的张量积。

原文的六面板比较图可在SciPy 原图查看。本稿配图是后面的轴与批次示意图，未重新运行上述绘图代码，也未制作伪造的数值结果图。

函数接口与对象接口

如果不想显式创建类实例，可以用便利函数 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

这里比较的是两种接口默认的线性插值结果，绝对容差为 1e-15。它没有比较不同插值策略，也不是本次环境中的测试报告。

只有一个节点的轴如何处理区间外查询

数据有时实际上位于 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.])

在这个例子里，bounds_error=False 允许查询区间外的点。当 fill_value=None 时，沿单点轴外推，所以三个查询得到同样的值 2；改成 -101 后，只有轴上的查询仍为 2，其余返回填充值。这个例子不能被理解为任意维度的外推都等价于边缘常数延伸。

原文另外提醒：若输入维度使用难以直接比较的单位，数值量级又相差很多，插值可能出现数值伪影。应考虑先缩放坐标数据，再插值。

向量输出放在尾部批次维度

现在设函数的输入和输出都是向量，即 f(x)=y，且 x 与 y 的长度不必相同。我们希望在 x 的网格上采样并一起插值多个输出分量。RegularGridInterpolator 允许 values 带尾部维度，正好可以表示这些批次分量。

它将 values 的前 N 维解释成数据维度，对应输入网格；其后的维度才是批次轴。这与很多 NumPy 广播场景中前导批次维度的习惯不同，不能把分量维随意放在开头。

n = 5  # 批次分量数

# 建立三维网格
x1 = np.linspace(-np.pi, np.pi, 10)
x2 = np.linspace(0.0, np.pi, 15)
x3 = np.linspace(0.0, np.pi/2, 20)
points = (x1, x2, x3)

def f(x1, x2, x3, n):
    lst = [np.sin(np.pi*x1/2) * np.exp(x2/2) + x3 + i
           for i in range(n)]
    return np.asarray(lst)

X1, X2, X3 = np.meshgrid(x1, x2, x3, indexing="ij")
values = f(X1, X2, X3, n)
values.shape
# 原文输出：(5, 10, 15, 20)

# 把批次分量轴移动到最后
values = np.moveaxis(values, 0, -1)
values.shape
# 原文输出：(10, 15, 20, 5)
rgi = RegularGridInterpolator(points, values)

x = np.asarray([0.2, np.pi/2.1, np.pi/4.1])
rgi(x).shape
# 原文输出：(1, 5)

这里的五个函数在同一个三维网格上求值。原始函数返回的形状是 (5,10,15,20)，移动轴后才成为插值器所需的 (10,15,20,5)。查询一个三维坐标时，返回形状为 (1,5)，包含五个分量。

也可以有多个批次维度。结果形状来自查询点的形状再接上批次形状；在这个单点特例中，查询按一个点处理，因此前部为 (1,)，后面接 (5,)。对通常的查询数组来说，最后一维是坐标分量，不直接保留成输出的数据轴。

等间距笛卡尔网格的另一个选择

对于整数坐标的笛卡尔网格，例如图像重采样，规则网格插值器不一定是最合适的选择，可以考虑 scipy.ndimage.map_coordinates。如果数据坐标是浮点数，但各轴都等间距，可以把它包装成类似 RegularGridInterpolator 的接口。

下面是原文给出的简化包装，源自 Johannes Buchner 的 regulargrid 包。它把物理坐标映射成像素索引坐标：

from scipy.ndimage import map_coordinates

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 维空间中的坐标。
        """
        # 转成 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)

各轴上的换算公式是 (n−1)×(坐标−最小值)/(最大值−最小值)。因此这个简化实现要求输入轴升序且等间距；不能把前面允许的不等间距网格直接交给它。它也没有实现前述尾部批次维度的完整接口，不能把所有插值器行为都当作可以互换。

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')
cgi([[1.5, 1.5], [3.5, 2.6]])
# 原文输出：array([ 9. , 64.9])

这个线性例子得到相同结果，但包装使用的是 map_coordinates 的边界条件，因此 cubic 和 quintic 可能与 RegularGridInterpolator 不同。其他边界设置和参数应查阅 map_coordinates 文档。

编者静态审查：示例只构造本地数组并绘图，没有外部下载、秘密或破坏性命令。包装代码没有检验升序、等间距、输入维数或 hi == lo；单点轴会使坐标归一化出现零分母，不能套用前面单点轴的示例。本文保留原简化实现并明确边界，没有声称已经修复或通过运行验证。

来源：SciPy 官方教程。包装代码原始归属：Johannes Buchner / regulargrid。页面版权归 SciPy community；中文译写及原创配图由编者整理。

Copyright (c) 2001-2002 Enthought, Inc. 2003, SciPy Developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above
   copyright notice, this list of conditions and the following
   disclaimer in the documentation and/or other materials provided
   with the distribution.

3. Neither the name of the copyright holder nor the names of its
   contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Copyright (c) 2013, Johannes Buchner
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

    Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
    Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.