插值過不規則柵格

Question

所以，我有存儲緯度，經度，和在網格上的一些屬性值3個numpy的陣列 - 也就是說，我有LAT（Y，X），LON（Y，X），和，說溫度T（y，x），對於x和y的一些限制。網格不一定是規則的 - 事實上，它是三極的。

然後，我想要將這些屬性（溫度）值插值到一堆不同的緯度/經度點（存儲爲lat1（t），lon1（t），約10,000 t ...）實際的網格點。我試過matplotlib.mlab.griddata，但這需要太長時間（畢竟，它並不是真正爲我所做的而設計的）。我也嘗試scipy.interpolate.interp2d，但我得到一個MemoryError（我的網格大約400x400）。

是否有任何形式的光滑，最好是快速的方式呢？我不禁想到答案是顯而易見的......謝謝！

Answer 1

嘗試在SO inverse-distance-weighted-idw-interpolation-with-python描述逆距離加權的組合和 scipy.spatial.KDTree 。 Kd-trees 在二維3d ...中很好地工作，反距離加權是平滑和局部的，並且k =最近鄰居的數量可以變化以權衡速度/準確度。

Answer 2

我建議你考慮看看草（一種開源GIS包）插補功能（http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html）。它不是在Python中，但你可以重新實現它或與C代碼接口。

Answer 3

我是正確的思維你的數據網格是這個樣子（紅色是舊數據，藍色是新的內插數據）？

alt text http://www.geekops.co.uk/photos/0000-00-02%20%28Forum%20images%29/DataSeparation.png

這可能是一個稍微強力十歲上下的方法，但有關呈現您的現有數據作爲位圖（OpenGL的會做的顏色簡單的插值爲您配置和正確的選項，你可以渲染什麼數據應該是相當快的三角形）。然後，您可以在新點的位置採樣像素。

或者，你可以在空間上的排序第一組點，然後找到最接近的舊點周圍的新點，並根據距離對這些點進行插值。

Answer 4

有一堆的選擇這裏，哪一個是最好的取決於你的數據... 但是我不知道出的現成的解決方案，你

你說你的輸入數據來自三極數據。關於如何構建這些數據，主要有三種情況。

1. 從三維空間中的三維網格採樣，投影回第2d個LAT，LON數據。
2. 從三維空間中的2d網格採樣，投影到2d LAT LON數據。投影到二維LAT LON數據

最簡單的，這些在三極空間

非結構化數據是2而不是插在LAT LON空間「只是」改變你點回源空間和內插在那裏。

，對於1和2的工作原理是搜索，從三極空間映射到覆蓋你的採樣點的細胞的另一種選擇。 （您可以使用BSP或網格類型結構來加速此搜索）選取其中一個單元格，然後在其中進行插值。

最後還有非結構化插選項的堆..但他們往往是緩慢的。 我個人最喜歡的是使用最近N個點的線性插值，發現這些N個點可以再次用網格或BSP來完成。另一個不錯的選擇是Delauney對非結構化點進行三角測量，並在生成的三角形網格上進行插值。

就個人而言，如果我的網格是案例1，我會使用非結構化策略，因爲我擔心必須處理通過具有重疊投影的單元格進行搜索。選擇「正確」的單元格會很困難。

Answer 5

有一個nice inverse distance example by Roger Veciana i Rovira以及一些使用GDAL編寫geotiff的代碼，如果你進入。

這對於規則網格來說很粗糙，但是假設您將數據首先投影到具有pyproj或其他東西的像素網格，一直注意什麼投影用於數據。

他的算法的副本測試：

from math import pow 
from math import sqrt 
import numpy as np 
import matplotlib.pyplot as plt 

def pointValue(x,y,power,smoothing,xv,yv,values): 
    nominator=0 
    denominator=0 
    for i in range(0,len(values)): 
     dist = sqrt((x-xv[i])*(x-xv[i])+(y-yv[i])*(y-yv[i])+smoothing*smoothing); 
     #If the point is really close to one of the data points, return the data point value to avoid singularities 
     if(dist<0.0000000001): 
      return values[i] 
     nominator=nominator+(values[i]/pow(dist,power)) 
     denominator=denominator+(1/pow(dist,power)) 
    #Return NODATA if the denominator is zero 
    if denominator > 0: 
     value = nominator/denominator 
    else: 
     value = -9999 
    return value 

def invDist(xv,yv,values,xsize=100,ysize=100,power=2,smoothing=0): 
    valuesGrid = np.zeros((ysize,xsize)) 
    for x in range(0,xsize): 
     for y in range(0,ysize): 
      valuesGrid[y][x] = pointValue(x,y,power,smoothing,xv,yv,values) 
    return valuesGrid 


if __name__ == "__main__": 
    power=1 
    smoothing=20 

    #Creating some data, with each coodinate and the values stored in separated lists 
    xv = [10,60,40,70,10,50,20,70,30,60] 
    yv = [10,20,30,30,40,50,60,70,80,90] 
    values = [1,2,2,3,4,6,7,7,8,10] 

    #Creating the output grid (100x100, in the example) 
    ti = np.linspace(0, 100, 100) 
    XI, YI = np.meshgrid(ti, ti) 

    #Creating the interpolation function and populating the output matrix value 
    ZI = invDist(xv,yv,values,100,100,power,smoothing) 


    # Plotting the result 
    n = plt.normalize(0.0, 100.0) 
    plt.subplot(1, 1, 1) 
    plt.pcolor(XI, YI, ZI) 
    plt.scatter(xv, yv, 100, values) 
    plt.title('Inv dist interpolation - power: ' + str(power) + ' smoothing: ' + str(smoothing)) 
    plt.xlim(0, 100) 
    plt.ylim(0, 100) 
    plt.colorbar() 

    plt.show() 

Answer 6

有一個FORTRAN庫調用BIVAR，這是非常適合這個問題。通過一些修改，你可以使用f2py在python中使用它。

從描述：

BIVAR是FORTRAN90庫，內插散射二元數據，由Hiroshi阿克瑪。 BIVAR接受一組散佈在2D中的（X，Y）數據點以及相關的Z數據值，並且能夠構造與給定數據一致的平滑插值函數Z（X，Y），並且可以在飛機上的其他點進行評估。