我们需要对用户输入的多个 (x, y) 坐标在给定数据上进行插值。数据是一个 4x7 的数组，其中值根据 x 和 y 的范围进行组织。我们希望根据用户输入的 x 和 y 坐标从数据中提取相应的值。
| >=0 1 2 3 4 5 >=6
-------------------------------------------
>=09 <10 | 6.4 5.60 4.8 4.15 3.5 2.85 2.2
>=10 <11 | 5.3 4.50 3.7 3.05 2.4 1.75 1.1
>=11 <12 | 4.7 3.85 3.0 2.35 1.7 1.05 0.4
>=12 | 4.2 3.40 2.6 1.95 1.3 0.65 0.0
如果用户输入 x = 2.5 和 y = 9，模型应该返回 4.475。另一方面，如果用户输入 x = 2.5 和 y = 9.5，模型应该返回 3.925。
我们尝试使用 map_coordinates 函数来进行插值，因为它提供了将坐标映射到 x,y 范围的能力。
import numpy as np
from scipy.ndimage import map_coordinates
# 定义数组
z = np.array([[6.4, 5.60, 4.8, 4.15, 3.5, 2.85, 2.2],
[5.3, 4.50, 3.7, 3.05, 2.4, 1.75, 1.1],
[4.7, 3.85, 3.0, 2.35, 1.7, 1.05, 0.4],
[4.2, 3.40, 2.6, 1.95, 1.3, 0.65, 0.0]])
# 插值函数
def twoD_interpolate(arr, xmin, xmax, ymin, ymax, x1, y1):
"""
在二维空间中进行插值，具有“硬边界”
"""
nx, ny = arr.shape
x1 = np.array([x1], dtype=np.float)
y1 = np.array([y1], dtype=np.float)
# 如果 x1 超出边界，将其值设置为最接近的点，以便我们始终在范围内进行插值
x1[x1 > xmax] = xmax
x1[x1 < xmin] = xmin
# 如果 y1 超出边界，将其值设置为最接近的点，以便我们始终在范围内进行插值
y1[y1 > ymax] = ymax
y1[y1 < ymin] = ymin
# 将 x1 和 y1 转换为索引，以便我们可以对它们进行映射
x1 = (nx - 1) * (x1 - xmin) / (xmax - xmin)
y1 = (ny - 2) * (y1 - ymin) / (ymax - ymin)
y1[y1 > 1] = 2.0
return map_coordinates(arr, [y1, x1])
# 获取值函数
def test_val(x, y, arr):
return twoD_interpolate(arr, 0, 6, 9, 12, x, y)
# 测试代码
print test_val(4, 11, z) --> 3.00
print test_val(2, 10, z) --> 3.85
然而，测试结果是错误的。上面代码中给出的结果分别是 3.00 和 3.85，而正确的结果应该是 1.7 和 3.7。
2、解决方案
要解决这个问题，我们需要对 twoD_interpolate 函数进行修改。主要修改的地方包括：
-
将
np.atleast_1d(x1)和np.atleast_1d(y1)改为np.array([x1], dtype=np.float)和np.array([y1], dtype=np.float)。这可以确保x1和y1始终是 1 维数组，从而避免后续出现错误。 -
将
y1[y1 > 1] = 2.0改为y1[y1 > 1] = 1。这可以确保y1的值始终在 0 和 1 之间，从而避免插值时出现超出范围的情况。 -
将
order=0改为order=1。这可以提高插值的精度，使插值结果更加准确。 -
将
mode='nearest'改为mode='constant'。这可以避免在插值时出现超出范围的情况，从而使插值结果更加稳定。
修改后的代码如下：
def twoD_interpolate(arr, xmin, xmax, ymin, ymax, x1, y1):
"""
在二维空间中进行插值，具有“硬边界”
"""
ny, nx = arr.shape # 注意 ny 和 xy 的顺序
x1 = np.array([x1], dtype=np.float)
y1 = np.array([y1], dtype=np.float)
# 将坐标转换为与数组匹配。
x1 = (x1 - xmin) * (nx - 1) / float(xmax - xmin)
y1 = (y1 - ymin) * (ny - 1) / float(ymax - ymin)
# 需要 order=1 才能返回示例。
# mode='constant' 防止需要剪辑
return map_coordinates(arr, np.vstack((y1, x1)), order=1, mode='constant')
修改后的代码可以产生正确的结果：
print test_val(4, 11, z)
[ 1.7]
print test_val(2, 10, z)
[ 3.7]
print test_val(2.5, 9, z)
[ 4.475]
print test_val(2.5, 9.5, z)
[ 3.925]
# 甚至可以使用 1D numpy 数组
print test_val(np.arange(4),np.arange(4)+9,z)
[ 6.4 4.5 3. 1.95]
