Scipy map_coordinates bilinear interpolation compared to interp and IDL interpolate

idl-programming-language, interpolation, matplotlib-basemap, python, scipy

Solution

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]

Problem

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?

Original source