""" a small example of 2d interpolation on a map,
using the scipy function ndimage.map_coordinates
"""
import scipy as sp
from scipy.ndimage import map_coordinates

# Cities: an n x 2 array of [latitide,longitude] coordinates --
Paris = [48.9, 2.4]
Rome = [41.9, 12.5]
Greenwich = [51.5, 0]
Cities = sp.array([ Paris, Rome, Greenwich ])

# LatlongtoTemp: an array of temperatures at integer [lat,long] --
LatlongtoTemp = sp.array( range( 91*360 )) .reshape(( 91, 360 ))  # dummy
    # LatlongtoTemp[0,:] is temperature along the equator,
    # LatlongtoTemp[:,0] along the prime meridian through Greenwich,
    # LatlongtoTemp[90,0] at the north pole

order = 3  # 1: bilinear, 3 (the default): cubic spline

z = map_coordinates( LatlongtoTemp, Cities.T, order=order )
# --> temperatures at Paris, Rome and Greenwich -- approximately, depending on order

print "map_coordinates( LatlongtoTemp, Cities.T, order=%d )" % order
print z

If order is 0, map_coordinates rounds [lat,long] to the nearest integers: the temperature at Paris is approximated by LatlongtoTemp[49,2]. If 1, it does bilinear interpolation in the square with corners [48,2], [48,3], [49,2], [49,3] for Paris. If 2, it does quadratic spline interpolation over the 9 points LatlongtoTemp[48:50+1, 1:3+1]. And so on, up to order 5; the default is order=3, a cubic spline. (However, avoid order=2 unless you know what you're doing -- there are various definitions of quadratic local splines.) Order 1, bilinear, is much faster than 3, but may have visible seams between the square patches (is not C1 continuous).

The transpose Cities.T is used because map_coordinates takes columns, not rows. (If you forget the .T, map_coordinates will warn "RuntimeError: invalid shape for coordinate array".)

Going off the edge of an input array, e.g. LatlongtoTemp[51,-1] west of Greenwich, is a problem for any interpolation code; see the mode= option.

Map_coordinates can effortlessly interpolate vectors too, not just scalars: for example, LatLongtoColor, or LatLongtoWindgradient, or LatlongtoMyrecord.

Exercise: plot the 1d impulse response, i.e. interpolate [0,0,0,1,0,0,0], for orders 1, 2, 3.

See also:
scipy tutorial ndimage
scipy cookbook interpolation
Wikipedia Multivariate_interpolation
scipy.ndimage.zoom -- image zooming
scipy.interpolate -- f2py Fitpack
Chapter 3, Interpolation, in Press et al., Numerical Recipes (2007), isbn 9780521880688
For the reverse problem of turning scattered data to a regular grid, see griddata .

Last change: denis 6 oct