更新：在每个单元格上应用（矢量化）函数以插值网格

Question

我有一个问题。我使用这些SO线程this, this和that来到了现在。

我有一个DEM文件和来自气象站的坐标+数据。现在，我想使用我的DEM按照GIDS模型（本文中的12号模型）插值空气温度数据。为了选择站点，我想使用KDTree使用8个最近邻。

简而言之，（我认为）我想使用我的DEM的坐标和高程在每个单元格评估函数。

我已经开发了一个工作函数，它使用x、y作为输入来评估我的网格的每个值。详见我的IPython Notebook。

但现在需要对整个numpy数组进行操作。我有点理解我必须矢量化我的函数，以便我可以将其应用于Numpy数组，而不是使用双重循环。请查看我简化的代码，以评估我的函数在数组上使用for循环和试用一个使用numpy meshgrid的矢量化函数的方法。这是正确的方法吗？

>>> data = [[0.8,0.7,5,25],[2.1,0.71,6,35],[0.75,2.2,8,20],[2.2,2.1,4,18]]
>>> columns = ['Long', 'Lat', 'H', 'T']
>>> df = pd.DataFrame(data, columns=columns)
>>> tree = KDTree(zip(df.ix[:,0],df.ix[:,1]), leafsize=10)
>>> dem = np.array([[5,7,6],[7,9,7],[8,7,4]])
>>> print 'Ground points\n', df
Ground points
   Long   Lat  H   T
0  0.80  0.70  5  25
1  2.10  0.71  6  35
2  0.75  2.20  8  20
3  2.20  2.10  4  18
>>> print 'Grid to evaluate\n', dem
Grid to evaluate
[[5 7 6]
 [7 9 7]
 [8 7 4]]
>>> def f(x,y):
...     [see IPython Notebook for details]
...     return m( sum((p((d(1,di[:,0])),2)))**-1 ,
...            sum(m(tp+(m(b1,(s(pix.ix[0,0],longp))) + m(b2,(s(pix.ix[0,1],latp))) + m(b3,(s(pix.ix[0,2],hp)))), (p((d(1,di[:,0])),2)))) )
... 
>>> #Double for-loop
... 
>>> tp = np.zeros([dem.shape[0],dem.shape[1]])
>>> for x in range(dem.shape[0]):
...     for y in range(dem.shape[1]):
...         tp[x][y] = f(x,y)
... 
>>> print 'T predicted\n', tp
T predicted
[[ 24.0015287   18.54595636  19.60427132]
 [ 28.90354881  20.72871172  17.35098489]
 [ 54.69499782  43.79200925  15.33702417]]
>>> # Evaluation of vectorized function using meshgrid
... 
>>> x = np.arange(0,3,1)
>>> y = np.arange(0,3,1)
>>> xx, yy = np.meshgrid(x,y, sparse=True)
>>> f_vec = np.vectorize(f) # vectorization of function f
>>> tp_vec = f_vec(xx,yy).T
>>> print 'meshgrid\nx\n', xx,'\ny\n',yy
meshgrid
x
[[0 1 2]] 
y
[[0]
 [1]
 [2]]
>>> print 'T predicted using vectorized function\n', tp_vec
T predicted using vectorized function
[[ 24.0015287   18.54595636  19.60427132]
 [ 28.90354881  20.72871172  17.35098489]
 [ 54.69499782  43.79200925  15.33702417]]

编辑

我使用了%%timeit来检查实际数据，网格大小为100x100，结果如下：

#double loop
for x in range(100):
    for y in range(100):        
        tp[x][y] = f(x,y)
1 loops, best of 3: 29.6 s per loop

#vectorized
tp_vec = f_vec(xx,yy).T
1 loops, best of 3: 29.5 s per loop

两者都不太好。

- Mattijn

---

虽然我很欣赏你在呈现背景、参考资料等方面非常仔细，但是你的问题既非常具体，又需要耗费大量时间来回答（这会降低别人帮助你的积极性）。我可以建议你写一个包含嵌入式测试数据的最小化难点示例。 - Henry Gomersall

我已经更新了代码，制作了一个简化的通用版本。 - Mattijn

现在将其制作为最小工作示例 - 创建一个包含所有内容的单个文件，运行并演示问题，并发布该文件。删除所有与问题无关的内容。顺便说一句，尽量更加遵守pep8；目前您的代码很难解析（最重要的是，您的函数f很难阅读 - 将其拆分为多行并使用更清晰的变量名称）。 - Henry Gomersall

@Henry，感谢您的建议，我真的很感激。我的问题是如何在每个单元格上应用矢量化函数。所以我已经得到了答案。现在我正在努力优化我的函数，但这与这个问题不再相关。无论如何，再次感谢。 - Mattijn

哦，我一直以为那很难！在我的经验中，人们很少发现 np.vectorize 能完全满足他们的需求。 - Henry Gomersall

Answer 1

如果在网格中使用矢量化函数，建议使用与依赖数组形状相同的meshgrid。使用从meshgrid派生的组件来使用矢量化函数评估每个网格单元。像这样：

def f(x,y):
    '...some code...'
    single_value = array[x,y] # = dependent array (e.g. DEM)
    '...some code...'
    return z

x = np.arange(array.shape[0])
y = np.arange(array.shape[1])
xx, yy = np.meshgrid(x,y, sparse=True)

f_vec = np.vectorize(f) # vectorization of function f

tp_vec = f_vec(xx,yy).T

- Mattijn