在结构化的二维数据上使用Scipy插值，但在非结构化点上进行评估？

Question

我有以下最小代码，使用scipy.interpolate.interp2d对2D网格数据进行插值。

import numpy as np
from scipy import interpolate
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**2+yy**2)
f = interpolate.interp2d(x, y, z, kind='cubic')

现在，f可以用来评估其他点。问题是我想要评估的点是完全随机的，而不是形成规则网格的点。

# Evaluate at point (x_new, y_new), in total 256*256 points
x_new = np.random.random(256*256)
y_new = np.random.random(256*256)
func(x_new, y_new)

这会导致我的电脑运行时出错，似乎它将x_new和y_new视为网格，生成了一个65536x65536的评估矩阵，这不是我的目的。

运行时错误: 无法生成大小为65536x65536的输出(大小太大)

一种解决方法是逐个评估点，使用以下代码:

z_new = np.array([f(i, j) for i, j in zip(x_new, y_new)])

然而，它很慢！

%timeit z_new = np.array([f(i, j) for i, j in zip(x_new, y_new)])

1.26秒±46.3毫秒每个循环（平均值±7次运行的标准差，1个循环）

是否有更快的方法来评估随机点？

这里的更快指的是与下面的时间相当：

x_new = np.random.random(256)
y_new = np.random.random(256)

%timeit f(x_new, y_new)

同样的256 * 256 = 65536次计算，在我的电脑上所需时间：

1.21 毫秒 ± 39.6 微秒每次循环（平均值±7次运行的标准偏差，每次1000个循环）

它并不需要在1.21毫秒的速度下运行，121毫秒完全可以接受。

- Jiadong

Answer 1

你要找的函数是 scipy.interpolate.RegularGridInterpolator 给定一组点 (x,y,z)，其中 x 和 y 定义在一个规则网格上，它允许你对中间的 (x,y) 点进行采样，得到相应的 z 值。在你的情况下，代码如下：

import numpy as np
from scipy import interpolate
x = np.arange(-5.01, 5.01, 0.25)
y = np.arange(-5.01, 5.01, 0.25)

def f(x,y):
    return np.sin(x**2+y**2)

z = f(*np.meshgrid(x, y, indexing='ij', sparse=True))
func = interpolate.RegularGridInterpolator((x,y), z)

x_new = np.random.random(256*256)
y_new = np.random.random(256*256)
xy_new = list(zip(x_new,y_new))
z_new = func(xy_new)func(xy_new)

更多细节请参见https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.interpolate.RegularGridInterpolator.html

- Paul Coppin