5

我正在将同事IDL代码重写为python,并且提出了一些我感到困惑的差异。根据我发现的其他SO问题和邮件列表线程,如果您使用scipy.ndimage.interpolation.map_coordinates并指定order=1它应该执行双线性插值。当比较IDL代码(在GDL中运行)和python(map_coordinates)之间的结果时,我得到了不同的结果。然后我尝试使用mpl_toolkits.basemap.interp,并得到了与IDL代码相同的结果。下面是一个简单的例子,显示了什么是错的。有人可以帮我弄清楚我在做什么map_coordinatesorder=1不是双线性的吗?Scipy map_coordinates双线性插值与interp和IDL插值相比

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文档没有提到双线性,它实际上是双线性的吗?我错过了什么?

回答

7

使用map_coordinates时,由于数组的形状为(height, width),所以需要转置数组或将坐标更改为(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] 
+0

哇,我不能相信我没有想到这一点。我以为我试过了,非常感谢。 – daveydave400 2013-02-25 03:06:54