Scipy-interp2d Returned Function Sorts Input Argument Automatically And Undesirably

Following the documentation: import matplotlib.pyplot as plt from scipy import interpolate import numpy as np x = np.arange(-5.01, 5.01, 0.25) y = np.arange(-5.01, 5.01, 0.25) xx,

Solution 1:

Your f callable takes an assume_sorted parameter:

assume_sorted : bool, optional
    If False, values of `x` and `y` can be in any order and they are
    sorted first.
    If True, `x` and `y` have to be arrays of monotonically
    increasing values.

So, yes, the inputs are sorted internally if you haven't sorted them before hand. I don't see a way of getting the sorted coordinates back.

The x,y inputs to interp2d are also sorted before use. Apparently the interpolation calculations require sorted arrays.

You can recover a pre-sort order with a double argsort index

Make an array and shuffle it:

In [415]: xnew = np.arange(-10,11,2)
In [416]: xnew
Out[416]: array([-10,  -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10])
In [417]: np.random.shuffle(xnew)
In [418]: xnew
Out[418]: array([  0,   2,   6,  -2,  10,  -4,   8,  -8, -10,  -6,   4])

Get the recovery index:

In [419]: idx = np.argsort(np.argsort(xnew))
In [420]: idx
Out[420]: array([ 5,  6,  8,  4, 10,  3,  9,  1,  0,  2,  7], dtype=int32)

test it:

In [421]: np.sort(xnew)[idx]
Out[421]: array([  0,   2,   6,  -2,  10,  -4,   8,  -8, -10,  -6,   4])

Solution 2:

Based on hpaulj's answer, you could define a new class that returns the unsorted array, e.g.,

from scipy.interpolate import interp2d
import numpy as np

class unsorted_interp2d(interp2d):
    def __call__(self, x, y, dx=0, dy=0):
        unsorted_idxs = np.argsort(np.argsort(x))
        return interp2d.__call__(self, x, y, dx=dx, dy=dy)[unsorted_idxs]

And then, in your example, you would have:

x = np.arange(-5.01, 5.01, 0.25)
y = np.arange(-5.01, 5.01, 0.25)
xx, yy = np.meshgrid(x, y)
z = np.sin(xx+yy)
f = unsorted_interp2d(x, y, z, kind='cubic')

xnewrev=np.array(list(reversed(np.arange(-5.01, 5.01, 1e-2))))
f(xnewrev, 0)

# the expected output!
array([-0.96171273, -0.96171273, -0.96443103, ...,  0.96176018,
        0.9589498 ,  0.95603946])