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很可能源於這個大量的插值點。