BUG: RegularGridInterpolator.__call__ fails when grid points arrays aren't C-contiguous in v1.10

[gmgunter]

Describe your issue.

scipy.interpolate.RegularGridInterpolator.__call__() raises an exception when one or more arrays of grid coordinates used to construct the interpolator are non-C-contiguous.

This appears to be new behavior in v1.10.0. In previous versions, there was no such restriction that these arrays must be C-contiguous -- the interpolator would produce valid results instead of raising an exception.

Reproducing Code Example

from scipy.interpolate import RegularGridInterpolator
import numpy as np

def f(x, y, z):
    return 2 * x**3 + 3 * y**2 - z

x = np.linspace(1, 4, 11)
y = np.linspace(4, 7, 22)
z = np.linspace(7, 9, 33)
xg, yg ,zg = np.meshgrid(x, y[::2], z[::3], indexing='ij', sparse=True)
data = f(xg, yg, zg)

# Take a strided slice of `y` and `z` so that these input arrays are not
# C-contiguous.
interp = RegularGridInterpolator((x, y[::2], z[::3]), data)

pts = np.array([[2.1, 6.2, 8.3],
                [3.3, 5.2, 7.1]])
interp(pts)

Error message

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/scratch/gunter/miniconda3/envs/tmp/lib/python3.9/site-packages/scipy/interpolate/_rgi.py", line 331, in __call__
    indices, norm_distances = self._find_indices(xi.T)
  File "/scratch/gunter/miniconda3/envs/tmp/lib/python3.9/site-packages/scipy/interpolate/_rgi.py", line 499, in _find_indices
    return find_indices(self.grid, xi)
  File "_rgi_cython.pyx", line 95, in scipy.interpolate._rgi_cython.find_indices
  File "stringsource", line 660, in View.MemoryView.memoryview_cwrapper
  File "stringsource", line 350, in View.MemoryView.memoryview.__cinit__
ValueError: ndarray is not C-contiguous

SciPy/NumPy/Python version information

1.10.0 1.24.1 sys.version_info(major=3, minor=9, micro=15, releaselevel='final', serial=0)

[ev-br]

Thanks for the report!
This is a regression of 1.10.0, and it has been fixed in main (by gh-17717):

In [1]: from scipy.interpolate import RegularGridInterpolator
   ...: import numpy as np
   ...:
   ...: def f(x, y, z):
   ...:     return 2 * x**3 + 3 * y**2 - z
   ...:
   ...: x = np.linspace(1, 4, 11)
   ...: y = np.linspace(4, 7, 22)
   ...: z = np.linspace(7, 9, 33)
   ...: xg, yg ,zg = np.meshgrid(x, y[::2], z[::3], indexing='ij', sparse=True)
   ...: data = f(xg, yg, zg)
   ...:
   ...: # Take a strided slice of `y` and `z` so that these input arrays are not
   ...: # C-contiguous.
   ...: interp = RegularGridInterpolator((x, y[::2], z[::3]), data)
   ...:
   ...: pts = np.array([[2.1, 6.2, 8.3],
   ...:                 [3.3, 5.2, 7.1]])
   ...: interp(pts)
Out[1]: array([125.84142857, 146.32518367])

In [2]: import scipy

In [3]: scipy.__version__
Out[3]: '1.11.0.dev0+1253.fdc32b8'

Cross-linking the duplicate reports in #17727 and #17716. The fix is tentatively labeled for 1.10.1.