2d非结构化网格数据的插值

Question

我有一些data（x，y，z）位于非结构化网格上，我想插入数据用于可视化目的。

我已经试过scipy.interpolate.griddata，插值假设处处都是相同的值。之后，我尝试了scipy.interpolate.Rbf，但是这给我一个内存错误（请参阅下面的代码）。

是否有其他方法或其他选项可以改善结果？

结果 - >

我的代码

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.interpolate import griddata, Rbf 

x, y, z = np.loadtxt('stackoverflow_example-data') 
# griddata 
points = np.reshape(np.array([x, y]),(z.size, 2)) 
grid_x, grid_y = np.mgrid[x.min():x.max():1000j,y.min():y.max():1000j] 
counts_I_grid_1 = griddata(points, z, (grid_x, grid_y), method='nearest') 
counts_I_grid_2 = griddata(points, z, (grid_x, grid_y), method='linear', fill_value=0) 
counts_I_grid_3 = griddata(points, z, (grid_x, grid_y), method='cubic', fill_value=0) 

# Rbf -- fails due to memory error 
#rbf = Rbf(x,y,z) 
#counts_I_Rbf = rbf(grid_x,grid_y) 

回溯（最近通话最后一个）： 文件 “/path/code.py” ，第14行 rbf = Rbf（x，y，z） 文件“/[...]/p （自我.xi，self.xi） 文件“/[...]/python3” .4/site-packages/scipy/interpolate/rbf.py“，第222行，在_call_norm中 return self.norm（x1，x2） 文件”/[...]/python3.4/site-packages/scipy /interpolate/rbf.py”，线114，在_euclidean_norm 返回SQRT（（（X1 - ×2）** 2）的.sum（轴= 0）） 的MemoryError

# plot the result 
fig = plt.figure() 

ax1 = plt.subplot(2,2,1) 
plt.title('Data') 
plt.gca().set_aspect((x.max() - x.min())/(y.max() - y.min())) 
plt.scatter(x, y, c=z, s=2, edgecolor='', marker=',') 
plt.colorbar(ax=ax1) 
plt.xlim(x.min(), x.max()) 
plt.ylim(y.min(), y.max()) 
plt.xticks([20.7,20.9,21.1,21.3]) 
plt.ticklabel_format(useOffset=False) 

ax2 = plt.subplot(2,2,2) 
plt.title('nearest') 
plt.imshow(counts_I_grid_1.T, origin='lower', 
      extent=(x.min(), x.max(), y.min(), y.max()), 
      aspect=(x.max() - x.min())/(y.max() - y.min()), 
      vmin=0,vmax=36) 
plt.colorbar(ax=ax2) 
plt.xticks([20.7,20.9,21.1,21.3]) 
plt.ticklabel_format(useOffset=False) 

ax2 = plt.subplot(2,2,3) 
plt.title('linear') 
plt.imshow(counts_I_grid_2.T, origin='lower', 
      extent=(x.min(), x.max(), y.min(), y.max()), 
      aspect=(x.max() - x.min())/(y.max() - y.min()), 
      vmin=0,vmax=36) 
plt.colorbar(ax=ax2) 
plt.xticks([20.7,20.9,21.1,21.3]) 
plt.ticklabel_format(useOffset=False) 

ax2 = plt.subplot(2,2,4) 
plt.title('cubic') 
plt.imshow(counts_I_grid_3.T, origin='lower', 
      extent=(x.min(), x.max(), y.min(), y.max()), 
      aspect=(x.max() - x.min())/(y.max() - y.min()), 
      vmin=0,vmax=36) 
plt.colorbar(ax=ax2) 
plt.xticks([20.7,20.9,21.1,21.3]) 
plt.ticklabel_format(useOffset=False) 
plt.tight_layout() 
plt.show() 

Answer 1

你的问题是由于到一个足够危险的微妙的错误，我相信这是值得一个完整的回答。

考虑这条线：

points = np.reshape(np.array([x, y]),(z.size, 2)) 

由于您的输入数组是一维，你想转换成[x,y]形 “（某事，2）”。请注意，它是有效的说

points = np.reshape(np.array([x, y]),(-1, 2)) 

让numpy推断你缺少的维度，这将仍然不是你想要的。当构造的2D阵列

np.array([x, y]) 

你被行定义矩阵行，从而产生形状的阵列“（2，东西）”。当你调用reshape时，它会默认逐行读取元素，因为numpy按照主要顺序存储数组（比如C/C++，而不像Fortran和MATLAB）。这意味着生成的两列数组将首先包含所有的值，然后是所有的值，而不是在每行中包含(x,y)对。

你实际上想要做的是交换你的数组的维度，而不接触它的结构。这意味着你必须转置你的矩阵。这意味着您必须使用

points = np.array([x, y]).T 
# or np.transpose([x,y]) 

改为。请注意，虽然这有相同shape你原来，它在正确的顺序中的元素：

In [320]: np.reshape(np.array([x,y]),(-1,2)) 
Out[320]: 
array([[ 20.7  , 20.702 ], 
     [ 20.704 , 20.706 ], 
     [ 20.708 , 20.71 ], 
     ..., 
     [ 852.356964, 852.356964], 
     [ 852.356964, 852.356964], 
     [ 852.356964, 852.356964]]) 

In [321]: np.array([x,y]).T 
Out[321]: 
array([[ 20.7  , 852.357235], 
     [ 20.702 , 852.357235], 
     [ 20.704 , 852.357235], 
     ..., 
     [ 21.296 , 852.356964], 
     [ 21.298 , 852.356964], 
     [ 21.3  , 852.356964]]) 

这将解决x/y点和z您的样品之间的矛盾，并产生预期的结果。根据我的经验，reshape-几乎从未被要求。通常你需要将一个ndarray变成一个1d之一，但是给定数组的方法是最好的。

结果证明：左边，你的原始三次插值;权，points版本修正：

注意为@MaxNoe suggested，我减少你的插值网格大小200×200。正如他暗示的那样，你的记忆错误与Rbf很可能源于这个大量的插值点。