Unexpected result when interpolating with Scippy ndimage.map_coordinates

Question

I would like to interpolate several user inputs (x, y) from the following data:

  | >=0 1 2 3 4 5 >=6 ------------------------------------------- >=09 <10 | 6.4 5.60 4.8 4.15 3.5 2.85 2.2 >=10 <11 | 5.3 4.50 3.7 3.05 2.4 1.75 1.1 >=11 <12 | 4.7 3.85 3.0 2.35 1.7 1.05 0.4 >=12 | 4.2 3.40 2.6 1.95 1.3 0.65 0.0 

If the user enters x = 2.5 and y = 9 , the model should return 4.475 . On the other hand, if the user enters x = 2.5 and y = 9.5 , the model should return 3.925 .

I used map_coordinates as it provides the ability to map coordinates to a range of x, y

Here is what I have done so far:

 import numpy as np from scipy.ndimage import map_coordinates # define array z = 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]]) # function to interpolate def twoD_interpolate(arr, xmin, xmax, ymin, ymax, x1, y1): """ interpolate in two dimensions with "hard edges" """ nx, ny = arr.shape x1 = np.array([x1], dtype=np.float) y1 = np.array([y1], dtype=np.float) # if x1 is out of bounds set its value to its closest point, so that we're always # interpolating within the range x1[x1 > xmax] = xmax x1[x1 < xmin] = xmin # if y1 is out of bounds set its value to its closest point, so that we're always # interpolating within the range y1[y1 > ymax] = ymax y1[y1 < ymin] = ymin # convert x1 and y1 to indices so we can map over them x1 = (nx - 1) * (x1 - xmin) / (xmax - xmin) y1 = (ny - 2) * (y1 - ymin) / (ymax - ymin) y1[y1 > 1] = 2.0 return map_coordinates(arr, [y1, x1]) # function to get the value def test_val(x, y, arr): return twoD_interpolate(arr, 0, 6, 9, 12, x, y) 

Code testing

 print test_val(4, 11, z) --> 3.00 print test_val(2, 10, z) --> 3.85 

These results are incorrect since they should return 1.7 and 3.7 respectively

Any ideas or thoughts on what I did wrong?

Answer 1

Using your function:

 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) # Change coordinates to match your array. x1 = (x1 - xmin) * (nx - 1) / float(xmax - xmin) y1 = (y1 - ymin) * (ny - 1) / float(ymax - ymin) # order=1 is required to return your examples. # mode='nearest' prevents the need for clip return map_coordinates(arr, np.vstack((y1, x1)), order=1, mode='nearest') # function to get the value def test_val(x, y, arr): return twoD_interpolate(arr, 0, 6, 9, 12, x, y) 

Now run some tests:

 print test_val(4, 11, z) [ 1.7] print test_val(2, 10, z) [ 3.7] print test_val(2.5, 9, z) [ 4.475] print test_val(2.5, 9.5, z) [ 3.925] #Can even use 1D numpy arrays now print test_val(np.arange(4),np.arange(4)+9,z) [ 6.4 4.5 3. 1.95] 

Explain:

np.atleast_1d is a function that ensures that your array has at least 1 dimension. np.array([x1]) will return a 2D array if x1 is a numpy array. This is undesirable.

 >>> np.atleast_1d(5) array([5]) >>> np.atleast_1d(np.arange(5)) array([0, 1, 2, 3, 4]) 

The setting order=1 refers to the order of spline interpolation. In the above example, you showed linear interpolation, therefore, 1 , if you want to take into account more values, you can increase it to achieve the desired effect.

np.vstack used to correctly position the x and y indices. In this terminology, map_coordinates needs data like:

 coords=[[y1,y2,y3,... [x1,y2,y3,...]] ycoords,xcoords=['y1','y2','y3'],['x1','x2','x3'] >>> np.vstack((xcoords,ycoords)) array([['y1', 'y2','y3'], ['x1', 'x2','x3']], dtype='|S2')