Scipy -- 3d griddata -- Why is it necessary to cast griddata xi argument to tuple?

Question

Why does the below call to griddata fail?

import scipy.interpolate
import numpy as np
grid_vals = np.meshgrid(*([np.linspace(-1,1,200)] * 3))
interp_vals = scipy.interpolate.griddata(np.random.randn(50,3), np.random.randn(50), grid_vals, 'linear')

The following exception occurs: ValueError: number of dimensions in xi does not match x

If i cast the xi (grid_vals) argument to tuple:

interp_vals = scipy.interpolate.griddata(np.random.randn(50,3), np.random.randn(50), tuple(grid_vals), 'linear') 

The error goes away. Why?

Answer 1

The basic why is that griddata passes both points and xi through a points = _ndim_coords_from_arrays(points) function whose documentation reads:

Convert a tuple of coordinate arrays to a (..., ndim)-shaped array.

and key action on tuples is:

p = np.broadcast_arrays(*points)

anything else, including a list, is just converted to an array:

points = np.asanyarray(points)

The actual interpolation expects arrays with the '3d' dimension last.

So your list of 3 (200,200,200) arrays becomes an array of (3,200,200,200) shape. But your points array is (50,3). The number of dimensions in xi does not match x message results from 200 not matching 3.

The griddata documentation is clear about points, less so for xi. But its example uses (x, Y) using arrays from mgrid.

So this would work:

 X, Y, Z = np.meshgrid(*([np.linspace(-1,1,200)] * 3))
 interp_vals = scipy.interpolate.griddata(points, values, (X,Y,Z), 'linear')

Another way of generating the required array from your meshgrid list is to make it an array, and roll the 1st dimension

grid_vals = np.rollaxis(np.array(grid_vals),0,4)

Another way of generating a mesh is np.ix_, which returns an open mesh in the form of tuple. An open mesh like this does need the broadcasting.

A single point would be interpolated with either:

interpolate.griddata(points,values,[[[[0,0,0]]]],'linear')
interpolate.griddata(points,values,([0],[0],[0]),'linear')

See the reaction to John's 4123 pull request has more discussion about the whys.

Answer 2

It's failing because the underlying scipy interpolate module has code that correctly handles multi-dimensional tuples and not multi-dimensional lists:

https://github.com/scipy/scipy/blob/master/scipy/interpolate/interpnd.pyx#L167

I've created issue 4123 and pull request 4124 to fix the issue, at least for lists.