使用Python和Gdal对点数据进行IDW插值

Question

我有一个包含纬度、经度和降雨量信息的 CSV 文件。我想插入这些点并创建 .tiff 文件。任何人都可以建议我最简单的方法吗？

我正在尝试使用

gdal_grid

。我对在 Python 中使用 GDAL 非常陌生。

Answer 1

这其实是几个问题。假设您有一些经纬度的分散数据，您将构建您想要进行估计的所有位置（Tiff 图像像素的所有纬度和经度）。

一旦您拥有了，您就可以使用任何解决方案对您的数据进行 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 公式 计算距离（该公式根据球体上两点的经度和纬度给出了它们之间的大圆距离）。再次注意，您可以找到半正矢距离的各种解决方案，例如这个：

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

最后，一旦准备好数组，您应该使用 GDAL 构建 Tiff。为此，请检查以下问题，我引用了其解决方案的一部分：

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)