Why Is Scipy's Ndimage.map_coordinates Returning No Values Or Wrong Results For Some Arrays?

Code Returning Correct value but not always returning a value In the following code, python is returning the correct interpolated value for arr_b but not for arr_a. Event though, I

Solution 1:

This is actually my fault in the original question.

If we examine the position it is trying to interpolate twoD_interpolate(arr, 0, 1, 1400, 3200, 0.5, 1684) we get arr[ 170400, 0.1] as the value to find which will be clipped by mode='nearest' to arr[ -1 , 0.1]. Note I switched the x and y to get the positions as it would appear in an array.

This corresponds to a 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.

The issues comes in the stride. If we take an array that goes from 50 to 250:

>>>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 using map_coordinates:

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

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

Here is what the code should be:

deftwoD_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 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 for arrays that are either 1D or pseudo 1D. Will interpolate the x dimension only unless the input array is of the proper 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]