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: 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]
    

