Scipy map_coordinates bilinear interpolation compared to interp and IDL interpolate

Question

I'm in the process of rewriting a coworkers IDL code into python and am coming up with some differences that I'm confused about. According to other SO questions and mailing list threads I've found if you use scipy.ndimage.interpolation.map_coordinates and specify order=1 it is supposed to do bilinear interpolation. When comparing results between the IDL code (run in GDL) and python (map_coordinates) I got different results. Then I tried using mpl_toolkits.basemap.interp and I got the same result as the IDL code. Below is a simple example that shows what is wrong. Could someone help me figure out what I am doing wrong with map_coordinates or is order=1 not bilinear?

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

I am perfectly fine using interp, but I was curious why map_coordinates didn't return the same result. I noticed that the map_coordinates documentation doesn't mention bilinear, is it actually bilinear? What am I missing?

Answer 1

When use map_coordinates, you need transpose the array or change you coordinates to (y, x) format, because the shape of the array is (height, width).

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)

This will output:

[ 25.89351463]
[ 25.89351439]