使用python和gdal的點數據的IDW插值

Question

我有一個包含Lat，Long和Rainfall信息的CSV文件。我想插入這些點並創建tiff文件。任何人都可以建議我最簡單的方法來做到這一點。

我想使用gdal_grid。我在python中使用gdal非常新。

Answer 1

這實際上是幾個問題。假設你有一些分散的數據用於lats和longs，那麼你需要建立所有的位置，然後進行估計（所有的位置和長度都是Tiff圖像的像素）。

一旦你有，你可以（在另一個question使用最近的一個例子）使用的任何解決方案，圍繞在你的數據做IWD：

class Estimation(): 
    # IWD. Check: https://stackoverflow.com/questions/36031338/interpolate-z-values-in-a-3d-surface-starting-from-an-irregular-set-of-points/36037288#36037288 
    def __init__(self,lon,lat,values): 
     self.x = lat 
     self.y = lon 
     self.v = values 

    def estimate(self,x,y,using='ISD'): 
     """ 
     Estimate point at coordinate x,y based on the input data for this 
     class. 
     """ 
     if using == 'ISD': 
      return self._isd(x,y) 

    def _isd(self,x,y): 
     #d = np.sqrt((x-self.x)**2+(y-self.y)**2) 
     d = x.copy() 
     for i in range(d.shape[0]): 
      d[i] = haversine(self.x[i],self.y[i],x,y) 
     if d.min() > 0: 
      v = np.sum(self.v*(1/d**2)/np.sum(1/d**2)) 
      return v 
     else: 
      return self.v[d.argmin()] 

上面的代碼實際上是用於計算與距離Haversine formula（它給出了經度和緯度上球體上兩點之間的大圓距離）。最後

def haversine(lon1, lat1, lon2, lat2): 
    """ 
    Check: https://stackoverflow.com/questions/15736995/how-can-i-quickly-estimate-the-distance-between-two-latitude-longitude-points 
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees) 
    """ 
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2]) 
    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 
    c = 2 * asin(sqrt(a)) 
    km = 6367 * c 
    return km 

一旦你有你的陣列準備好，你應該剛剛建立的Tiff使用GDAL：再次通知你可以找到各種各樣的半正矢遠程解決方案，像this one。對於此檢查以下內容question爲此，我引用它的解決方案的一部分：

driver = gdal.GetDriverByName('GTiff') 
ds = driver.Create('output.tif',xsize, ysize, 1, gdal.GDT_Float32,) 
# this assumes the projection is Geographic lat/lon WGS 84 
srs = osr.SpatialReference() 
srs.ImportFromEPSG(4326) 
ds.SetProjection(srs.ExportToWkt()) 
gt = [ulx, xres, 0, uly, 0, yres ] 
ds.SetGeoTransform(gt) 
outband=ds.GetRasterBand(1) 
outband.SetStatistics(np.min(mag_grid), np.max(mag_grid), np.average(mag_grid), np.std(mag_grid)) 
outband.WriteArray(mag_grid)