understanding inputs and outputs on scipy.ndimage.map_coordinates

Question

I am trying to map a straight line on a set of points in a grid. the data is in a list of x, y, z coordinates. I think map_coordinates is what i want, however i do not undestand the form of the inputs and outputs... any help would be much appreciated.

list = [[x1, y1, z1]
        [x2, y2, z2]
        ...
        [xn, yn, zn]]

the values I am trying to look up are x and y values.

look_up_values= [[x1, y1]
                 [x2, y2]
                 ...
                 [xn, yn]]

My question, why does map_coordinates expect a (2,2) array and what information is in the output ( a 2 value list).

here is a workable example:

from scipy.ndimage.interpolation import map_coordinates
import numpy as np
in_data = np.array([[0.,0.,0.]
                    ,[0.,1.,.2]
                    ,[0.,2.,.4]
                    ,[1.,0.,.2]
                    ,[1.,3.,.5]
                    ,[2.,2.,.7]])
z = map_coordinates(in_data, np.array([[1.,1.],[1.,2.]]), order=1)
print z #I do not understand this output...
#[1. .2]

if i had to guess i would say its interpolating between 2 points on the grid, buy what then would the output mean?

Answer 1

The output of map_coordinates is an interpolation of the value of the original array at the coordinates you've specified.

In your example you input [[1, 1], [1, 2]]. This means you want the interpolated value at two locations: the point x=1,y=1 and x=1,y=2. It needs two arrays because each array is the x- and y-coordinates. I.e. there are two coordinates you've asked for: x-coordinates at 1,1 and y-coordinates at 1,2.

The specific example you've chosen is a bit confusing because it looks like [1, 1] corresponds to [x0, y0] this interpretation is incorrect - it actually corresponds to [x0, x1]. This becomes clear with more than 2 points: The general format is [[x0, x1, x2, ...], [y0, y1, y2, ...]].

Your inputs can be as long or short as you like, but the arrays must be the same length since they are coupled.

Answer 2

Let me try to answer first examining a step by step 1d interpolation case:

import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import numpy as np

### 1d example of interpolation ###

in_data_x = np.array([1., 2., 3., 4., 5., 6.])
in_data_y = np.array([1.5, 2., 2.5, 3.,  3.5,  4.])  # y = .5 x - 1
f = interp1d(in_data_x, in_data_y, kind='linear')

print(f)
# f in all of the points of the grid (in_data_x): output coincides with in_data_y

    
print(f(1), f(1.), f(1.5), f(2.), f(2.5), f(3.))
# f in a point outside the grid:
print(f(1.8))
# this is equal to y = .5 x - 1 for x = 1.8, up to some point.
assert round(0.5 * 1.8 + 1, ndigits=10) == round(f(1.8), ndigits=10)

# plot up to this point
xnew = np.arange(1, 6, 0.1)
ynew = f(xnew)
plt.plot(in_data_x, in_data_y, 'o', xnew, ynew, '-')
# close the image to move forward.
plt.show()

### another 1d example of interpolation ###

in_data_x = np.array([1., 2., 3., 4., 5., 6.])
in_data_y = np.array([-1.8, -1.2, -0.2, 1.2, 3., 5.2])  # y = .2 x**2 - 2
f = interp1d(in_data_x, in_data_y, kind='cubic')

print(f)
# f in all of the points of the grid (in_data_x): output coincides with in_data_y
print(f(1), f(1.), f(1.5), f(2.), f(2.5), f(3.))
# f in a point outside the grid:
print(f(1.8))
# this is equal to y = .2 x**2 - 2 for x = 1.8, up to some precision.
assert round(0.2 * 1.8 ** 2 - 2, ndigits=10) == round(f(1.8), ndigits=10)

# plot up to this point
xnew = np.arange(1, 6, 0.1)
ynew = f(xnew)
plt.plot(in_data_x, in_data_y, 'o', xnew, ynew, '-')
plt.show()

The function interp1d provides you an interpolator that gives you the value that interpolates with some kind of algorithm (in this case linear) the function passing through x = [1., 2., 3., 4., 5., 6.] y = [-1.8, -1.2, -0.2, 1.2, 3., 5.2].

map_coordinates does the same. When your data have more than one dimension. First main difference is that the results is not an interpolator, but it is an array. Second main difference is that the x coordinates are given by the matrix coordinates of the dimension of your data. Third difference is that the input must be given as a column vector. See this example

from scipy.ndimage.interpolation import map_coordinates
import numpy as np


in_data = np.array([[0., -1., 2.],
                    [2., 1., 0.],
                    [4., 3., 2.]])  # z = 2.*x - 1.*y

# want the second argument as a column vector (or a transposed row)
# see on some points of the grid:
print('at the point 0, 0 of the grid the function z is: ')
print(map_coordinates(in_data, np.array([[0., 0.]]).T, order=1))
print('at the point 0, 1 of the grid the function z is: ')
print(map_coordinates(in_data, np.array([[0., 1.]]).T, order=1))
print('at the point 0, 2 of the grid the function z is: ')
print(map_coordinates(in_data, np.array([[0., 2.]]).T, order=1))

# see some points outside the grid
print()
print('at the point 0.2, 0.2 of the grid, with linear interpolation z is:')
print(map_coordinates(in_data, np.array([[.2, .2]]).T, order=1))
print('and it coincides with 2.*.2 - .2')
print()
print('at the point 0.2, 0.2 of the grid, with cubic interpolation z is:')
print(map_coordinates(in_data, np.array([[0.2, .2]]).T, order=3)

Finally answering your question, you gave as input

in_data = np.array([[0., 0., 0.],
                    [0., 1., .2],
                    [0., 2., .4],
                    [1., 0., .2],
                    [1., 3., .5],
                    [2., 2., .7]])

That is a function z(x,y) computed on the grid given by the matrix coordinates: z(0, 0) = 0. z(2, 2) = .7 Asking

z = map_coordinates(in_data, np.array([[1., 1.], [1., 2.]]), order=1)

means asking z(1,1) and z(1,2), where the second input array is read by column.

z = map_coordinates(in_data, np.array([[.5, .5]]).T, order=1)

means asking z(0.5, 0.5). Note the transpose .T in the input. Hope it does make sense and it is helpful.