Scipy map_coordinates 双线性插值与 interp 和 IDL interpolate 比较

Question

我正在把同事的IDL代码重写成Python，但遇到了一些让我困惑的不同之处。根据我在其他StackOverflow问题和邮件列表中找到的信息，如果你使用scipy.ndimage.interpolation.map_coordinates并指定order=1，它应该是进行双线性插值。但我把IDL代码（在GDL中运行）和Python（使用map_coordinates）得到的结果进行了比较，发现结果不一样。然后我尝试使用mpl_toolkits.basemap.interp，结果和IDL代码一样。下面是一个简单的例子，展示了我哪里出错了。有人能帮我看看我在使用map_coordinates时做错了什么，或者order=1真的不是双线性插值吗？

from scipy.ndimage.interpolation import map_coordinates
from mpl_toolkits.basemap import interp
import numpy

in_data = numpy.array([[ 25.89125824,  25.88840675],[ 25.90930748,  25.90640068]], dtype=numpy.float32)

map_coordinates(in_data, [[0.0],[0.125]], order=1, mode='nearest')
# map_coordinates results in "25.89090157"
interp(in_data, numpy.array([0,1]), numpy.array([0,1]), numpy.array([0.0]), numpy.array([0.125]), order=1)
# interp results in "25.89351439", same as GDL's "25.8935" when printed

我使用interp完全没问题，但我很好奇为什么map_coordinates的结果不一样。我注意到map_coordinates的文档没有提到双线性插值，它真的算是双线性插值吗？我漏掉了什么吗？

Answer 1

使用 map_coordinates 时，你需要把数组转置，或者把你的坐标改成 (y, x) 的格式，因为数组的形状是 (高度, 宽度)。

from scipy.ndimage.interpolation import map_coordinates
from mpl_toolkits.basemap import interp
import numpy

in_data = numpy.array([[ 25.89125824,  25.88840675],[ 25.90930748,  25.90640068]], dtype=numpy.float32)

print map_coordinates(in_data.T, [[0.0],[0.125]], order=1, mode='nearest')
print interp(in_data, numpy.array([0,1]), numpy.array([0,1]), numpy.array([0.0]), numpy.array([0.125]), order=1)

这将会输出：

[ 25.89351463]
[ 25.89351439]