understanding inputs and outputs on scipy.ndimage.map_coordinates

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.

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


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.

Tags:

Python

Scipy