Why does Scipy ndimage.map_coordinates not return any values ​​or incorrect results for some arrays?

Question

Code Return The correct value, but not always returning a value

In the following code, python returns the correct interpolated value for arr_b , but not for arr_a .

The event, though, I looked at this problem for about a day, I'm really not sure what is happening.

For some reason, for arr_a, twoD_interpolate keeps returning [0], even if I play around or interfere with data and input.

How can I fix my code so that it actually interpolates by arr_a and returns the correct results?

 import numpy as np from scipy.ndimage import map_coordinates def twoD_interpolate(arr, xmin, xmax, ymin, ymax, x1, y1): """ interpolate in two dimensions with "hard edges" """ ny, nx = arr.shape # Note the order of ny and xy x1 = np.atleast_1d(x1) y1 = np.atleast_1d(y1) # Mask upper and lower boundaries using @Jamies suggestion np.clip(x1, xmin, xmax, out=x1) np.clip(y1, ymin, ymax, out=y1) # Change coordinates to match your array. x1 = (x1 - xmin) * (xmax - xmin) / float(nx - 1) y1 = (y1 - ymin) * (ymax - ymin) / float(ny - 1) # order=1 is required to return your examples. return map_coordinates(arr, np.vstack((y1, x1)), order=1) # test data arr_a = np.array([[0.7, 1.7, 2.5, 2.8, 2.9], [1.9, 2.9, 3.7, 4.0, 4.2], [1.4, 2.0, 2.5, 2.7, 3.9], [1.1, 1.3, 1.6, 1.9, 2.0], [0.6, 0.9, 1.1, 1.3, 1.4], [0.6, 0.7, 0.9, 1.1, 1.2], [0.5, 0.7, 0.9, 0.9, 1.1], [0.5, 0.6, 0.7, 0.7, 0.9], [0.5, 0.6, 0.6, 0.6, 0.7]]) arr_b = np.array([[6.4, 5.60, 4.8, 4.15, 3.5, 2.85, 2.2], [5.3, 4.50, 3.7, 3.05, 2.4, 1.75, 1.1], [4.7, 3.85, 3.0, 2.35, 1.7, 1.05, 0.4], [4.2, 3.40, 2.6, 1.95, 1.3, 0.65, 0.0]]) # Test the second array print twoD_interpolate(arr_b, 0, 6, 9, 12, 4, 11) # Test first area print twoD_interpolate( arr_a, 0, 500, 0, 2000, 0, 2000) print arr_a[0] print twoD_interpolate( arr_a_60, 0, 500, 0, 2000, 0, 2000)[0] print twoD_interpolate( arr_a, 20, 100, 100, 1600, 902, 50) print twoD_interpolate( arr_a, 100, 1600, 20, 100, 902, 50) print twoD_interpolate( arr_a, 100, 1600, 20, 100, 50, 902) ## Output [ 1.7] [ 0.] [ 0.7 1.7 2.5 2.8 2.9] 0.0 [ 0.] [ 0.] [ 0.] 

Code returning an invalid value:

 arr = np.array([[12.8, 20.0, 23.8, 26.2, 27.4, 28.6], [10.0, 13.6, 15.8, 17.4, 18.2, 18.8], [5.5, 7.7, 8.7, 9.5, 10.1, 10.3], [3.3, 4.7, 5.1, 5.5, 5.7, 6.1]]) twoD_interpolate(arr, 0, 1, 1400, 3200, 0.5, 1684) # above should return 21 but is returning 3.44 

Answer 1

Actually this is my mistake in the original question.

If we look at the position that he is trying to interpolate twoD_interpolate(arr, 0, 1, 1400, 3200, 0.5, 1684) , we get arr[ 170400, 0.1] as the value for the search, which will be truncated using mode='nearest' to arr[ -1 , 0.1] . Note. I switched x and y to get the positions to be displayed in the array.

This corresponds to interpolation from the values arr[-1,0] = 3.3 and arr[-1,1] = 4.7 , so the interpolation looks like 3.3 * .9 + 4.7 * .1 = 3.44 .

Problems go a step. If we take an array from 50 to 250 in size:

 >>> a=np.arange(50,300,50) >>> a array([ 50, 100, 150, 200, 250]) >>> stride=float(a.max()-a.min())/(a.shape[0]-1) >>> stride 50.0 >>> (75-a.min()) * stride 1250.0 #Not what we want! >>> (75-a.min()) / stride 0.5 #There we go >>> (175-a.min()) / stride 2.5 #Looks good 

We can view this with map_coordinates :

 #Input array from the above. print map_coordinates(arr, np.array([[.5,2.5,1250]]), order=1, mode='nearest') [ 75 175 250] #First two are correct, last is incorrect. 

So we really need (x-xmin) / stride , for the previous examples, the step was 1, so it did not matter.

Here is what the code should do:

 def twoD_interpolate(arr, xmin, xmax, ymin, ymax, x1, y1): """ interpolate in two dimensions with "hard edges" """ arr = np.atleast_2d(arr) ny, nx = arr.shape # Note the order of ny and xy x1 = np.atleast_1d(x1) y1 = np.atleast_1d(y1) # Change coordinates to match your array. if nx==1: x1 = np.zeros_like(x1.shape) else: x_stride = (xmax-xmin)/float(nx-1) x1 = (x1 - xmin) / x_stride if ny==1: y1 = np.zeros_like(y1.shape) else: y_stride = (ymax-ymin)/float(ny-1) y1 = (y1 - ymin) / y_stride # order=1 is required to return your examples and mode=nearest prevents the need of clip. return map_coordinates(arr, np.vstack((y1, x1)), order=1, mode='nearest') 

Note that a clip is not required with mode='nearest' .

 print twoD_interpolate(arr, 0, 1, 1400, 3200, 0.5, 1684) [ 21.024] print twoD_interpolate(arr, 0, 1, 1400, 3200, 0, 50000) [ 3.3] print twoD_interpolate(arr, 0, 1, 1400, 3200, .5, 50000) [ 5.3] 

Checking arrays that are 1D or pseudo 1D. Interpolates the measurement x only if the input array has the correct shape:

 arr = np.arange(50,300,50) print twoD_interpolate(arr, 50, 250, 0, 5, 75, 0) [75] arr = np.arange(50,300,50)[None,:] print twoD_interpolate(arr, 50, 250, 0, 5, 75, 0) [75] arr = np.arange(50,300,50) print twoD_interpolate(arr, 0, 5, 50, 250, 0, 75) [50] #Still interpolates the `x` dimension. arr = np.arange(50,300,50)[:,None] print twoD_interpolate(arr, 0, 5, 50, 250, 0, 75) [75]