Scipy map_coordinates interpolazione bilineare rispetto a interp e IDL interpolate

Question

Sono in procinto di riscrivere un codice IDL di colleghi in python e sto arrivando con alcune differenze di cui sono confuso. Secondo altre domande SO e thread di mailing list ho trovato se si utilizza scipy.ndimage.interpolation.map_coordinates e si specifica order=1 che si suppone faccia un'interpolazione bilineare. Quando ho confrontato i risultati tra il codice IDL (eseguito in GDL) e python (coordinate_carta) ho ottenuto risultati diversi. Poi ho provato a utilizzare mpl_toolkits.basemap.interp e ho ottenuto lo stesso risultato del codice IDL. Di seguito è riportato un semplice esempio che mostra ciò che è sbagliato. Qualcuno potrebbe aiutarmi a capire cosa sto sbagliando con map_coordinates o è order=1 non bilineare?

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 

Sto perfettamente bene usando interp, ma ero curioso di sapere perchè map_coordinates non ha prodotto lo stesso risultato. Ho notato che la documentazione map_coordinates non menziona bilineare, in realtà è bilineare? Cosa mi manca?

Answer 1

Quando si utilizza map_coordinates, è necessario trasporre l'array o modificare le coordinate in formato (y, x), poiché la forma dell'array è (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) 

Questo stamperà:

[ 25.89351463] 
[ 25.89351439]