Scipy——3d网格——为什么需要将网格xi参数转换为tuple

Question

为什么下面对gridata的调用失败？

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')

出现以下异常：ValueError：xi中的维度数与x 不匹配

如果我将xi（grid_vals）参数强制转换为tuple:

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

错误消失了。为什么？

Answer 1

基本原因是griddata通过一个points = _ndim_coords_from_arrays(points)函数传递points和xi，该函数的文档为：

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

元组上的关键操作是：

p = np.broadcast_arrays(*points)

其他任何东西，包括列表，都会被转换成一个数组：

points = np.asanyarray(points)

实际插值需要最后具有"3d"维度的数组。

因此，您的3个(200,200,200)阵列的列表变成了(3,200,200,200)形状的阵列。但是你的points数组是(50,3)。number of dimensions in xi does not match x消息是由于200与3不匹配而产生的。

griddata文档对points很清楚，对xi则不那么清楚。但其示例使用CCD_ 15，使用来自CCD_ 16的数组。

所以这是可行的：

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

从meshgrid列表生成所需数组的另一种方法是使其成为数组，并滚动第一维

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

生成网格的另一种方法是np.ix_，它以元组的形式返回一个开放网格。像这样的开放式网格确实需要广播。

将使用以下任一项对单个点进行插值：

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

请参阅对John的4123拉取请求的反应，了解更多关于原因的讨论。

Answer 2

它失败了，因为底层的scipy插值模块有正确处理多维元组而不是多维列表的代码：

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

我已经创建了问题4123和拉取请求4124来修复这个问题，至少对于列表是这样。