从给定数据源和包含NoData值的网格点进行IDW插值

Question

我一直在寻找一个Python脚本来对我的数据进行插值处理。我找到了一段可以解决我问题的代码。现在我有几个问题想要修改下面的代码。1）我的数据是按照data.txt中的格式组织的，里面有经纬度和数值。我希望代码能够读取我的数据并处理这些文本数据。2）我还有经纬度的网格点数据。我想从这些网格点数据中进行插值，并将插值后的数值与提供的网格点数据合并。3）我的数据中有一些NoData/NULL值，因此我希望代码只从有数据的站点进行插值。

data.txt

XV   YV   v1  v2  v2  v4  v5  v6  v7
10  10  1   2   4   3   NA      4
10  15  4   4   3   NA  NA      NA
30  35  NA  NA  1   NA  5       18
5   20  4   NA  4   3   10      NA
15  15  NA  5   4   NA  NA      5
25  10  7   8   7   5   10      NA

我的网格点数据看起来是这样的：

grid.txt

 X  Y
10  10
10  15
10  20
10  25
10  30
10  35
10  40
15  10
15  15
15  20
15  25
15  30
15  35
15  40
20  10
20  15
20  20
20  25
20  30
20  35
20  40
25  10
25  15
25  20
25  25
25  30
25  35
25  40
30  10
30  15
30  20
30  25
30  30
30  35
30  40
35  10
35  15
35  20
35  25
35  30
35  35
35  40
40  10
40  15
40  20
40  25
40  30
40  35
40  40

我正在尝试的代码是：

#! /usr/bin/python

 from math import pow  
 from math import sqrt  
 import numpy as np  

 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)  
     print ZI

Answer 1

你真的确定要用IDW吗？IDW有个不太好的特点，就是在插值的时候，你必须考虑所有的点。这在很多情况下插值效果也不太好。所以，如果你只是想要对你的数据进行一个不错的插值，建议你看看scipy.interpolate.griddata，这可能正是你需要的。

另一方面，如果你想用IDW来插值一些未知点，它可以被简化成比较短的代码（但要记住，这仍然是O(m*n)，其中m和n分别是未知点和已知点的数量）。

import numpy as np

# inverse power to use
pow = 2

# some points:
known_pts = np.array([[0,0], [1,0], [3,0], [4,0],
                      [0,1], [2,1], [3,1], [4,1],
                      [1,3], [2,3], [3,3], [4,3]])
# and some data associated with them
known_values = known_pts[:,0] + known_pts[:,1]

# unknown points
unknown_pts = np.array([[2,0], [1,1],  [0,2], [1,2], [2,2], [3,2], [4,2], [0,3]])

# calculate all distances from unknown points to known points:
# (the array will have as many rows as there are unknown points and as many columns
#  as there are known points)
dist = np.sqrt((known_pts[:,0][None,:]-unknown_pts[:,0][:,None])**2 + (known_pts[:,1][None,:]-unknown_pts[:,1][:,None])**2)

# calculate the inverse distances, use a small epsilon to avoid infinities:
idist = 1. / (dist + 1e-12)**pow

# calculate the weighted average for each column with idist as the weight function
unknown_values = np.sum(known_values[None,:] * idist, axis=1) / np.sum(idist, axis=1)

当我们把这个例子画出来时（斑点区域反映了值，红色是已知点，蓝色是未知点），我们得到：