SciPy: Why is the input of 2D interpolation sorted internally?

Question

I hae a problem with scipy.interpolate.interp2d.

For sorted input, the interpolation is OK.

When I ask to get the interpolation values for an unsorted array, I get output as if it is sorted internally by SciPy. Why is that?

The way around is to get interpolation values in a loop.

Here is my demonstration code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
SciPy interp2d test.
Why is the input of 2D interpolation sorted internally?
'''
import matplotlib.pyplot as plt
import scipy.interpolate as itp
import numpy as np

def fMain():
    
 nx=11
 ny=21
 
 ax=np.linspace(0,1,nx)
 mx=np.empty((nx,ny))
 for i in range(ny):
  mx[:,i] = ax
 pass
  
 ay=np.linspace(0,1,ny)
 my=np.empty((nx,ny))
 for i in range(nx): # can I do this without loop?
  my[i,:] = ay
 pass   

 mz=np.empty((nx,ny))
 mz=mx**2 + my**3

 f2Di = itp.interp2d( mx, my, mz, kind='linear')
 #this provides identical results, ok
 #f2Di = itp.interp2d( ax, ay, mz.transpose(), kind='linear')

 if True :
  # just to check the interpolation   
  mzi = f2Di(ax,ay)   
  fig = plt.figure()
  axis = fig.add_subplot(projection='3d')
  axis.plot_wireframe( mx, my, mz )
  axis.scatter(mx, my, mzi.transpose(), marker="o",color="red")
  axis.set_xlabel("x")
  axis.set_ylabel("y")
  axis.set_zlabel("z")
  plt.tight_layout()
  plt.show()
  plt.close()
 pass 
 
 if True:
  y = 0.5
  az = f2Di( ax , y )
  axf=np.flip(ax)
  azf1 = f2Di( axf , y )
  azf2 = np.empty(nx)
  for i in range(nx):
   azf2[i] = f2Di( axf[i] , y )
  pass    
  plt.plot(ax,az,label="Normal",linewidth=3,linestyle="dashed")
  plt.plot(axf,azf1,label="Reversed")
  plt.plot(axf,azf2,label="Reversed loop")
  plt.legend()
  plt.xlabel("x")
  plt.ylabel("z")
  plt.tight_layout()
  plt.show()
  plt.close()  
 pass
pass    

if __name__ == "__main__":
 fMain()   
pass 

Answer 1

There's a builtin function for make grids like your mx,my:

In [68]: I,J = np.meshgrid(ay,ax)
In [69]: I.shape
Out[69]: (11, 21)

In [71]: np.allclose(I,my)
Out[71]: True

In [72]: np.allclose(J,mx)
Out[72]: True

alternatively you could have assigned the values with broadcasting

In [76]: my = np.empty((nx,ny)); my[:]=ay
In [78]: mx = np.empty((nx,ny)); mx[:]=ax[:,None]

The interp2d docs say that input arrays are flattened, even if input as 2d. And that the x,y can be the coordinates as in ax,ay; they don't have to be constructed from the full grid. So the 2 ways of setting up the f2Di are equivalent.

Full documentation for the use of f2Di(x,y) is

https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.__call__.html#scipy.interpolate.interp2d.__call__

It explicitly states that the inputs, x,y have to sorted, or it will do it for you.

One interpolation:

In [86]: mzi = f2Di(ax,ay)
In [87]: mzi.shape
Out[87]: (21, 11)

Another with the inputs reversed:

In [89]: azf1 = f2Di(ax[::-1], ay[::-1] )
In [90]: azf1.shape
Out[90]: (21, 11)
In [91]: np.allclose(mzi, azf1)
Out[91]: True

As you note, and attempt to show with a lot of plotting code, the results are the same - inputs have been sorted before they are used to interpolate.

If I falsely tell it that the coordinates are sorted:

In [94]: azf1 = f2Di(ax[::-1], ay[::-1] , assume_sorted=True)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
...

File ~\anaconda3\lib\site-packages\scipy\interpolate\_fitpack_impl.py:1054, in bisplev(x, y, tck, dx, dy)
   1052 z, ier = _fitpack._bispev(tx, ty, c, kx, ky, x, y, dx, dy)
   1053 if ier == 10:
-> 1054     raise ValueError("Invalid input data")
   1055 if ier:
   1056     raise TypeError("An error occurred")

ValueError: Invalid input data

Note that the error was raised by a function in _fitpack. The name implies that this using some sort of compiled interpolation code, a library that is probably written in C or Fortran. I'm not a developer, but I can imagine that it's easiest to write such code assuming that the inputs are sorted. Such shared libraries work best when they have clear, and relatively simple, expectations regarding the inputs.

Answer 2

To answer another question (from a comment in the code):

ax=np.linspace(0,1,nx)
mx=np.empty((nx,ny))
for i in range(ny):
    mx[:,i] = ax

(and similar for ay).

Can I do this with a loop?

Yes (technically no, since there'll be a C loop under the hood, but practially, yes). Use numpy.tile:

ax = np.linspace(0, 1, nx)
mx = np.tile(ax, (ny, 1)).T

---

And the np.empty doesn't make sense below: you allocate memory, but immediately (re)assign the variable to another value:

#mz = np.empty((nx, ny))   # This line is redundant
mz = mx**2 + my**3

This is also why np.empty has disappeared form the for-replacement code.