python - 使用python和gdal对点数据进行IDW插值

Question

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

我正在尝试使用 gdal_grid。我对在 python 中使用 gdal 很陌生。

Answer 1

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

一旦您掌握了这些解决方案，您就可以使用任何解决方案对您的数据进行 IWD(使用另一个 question 中的最新示例):

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 的距离(根据经度和纬度给出球体上两点之间的大圆距离)。再次注意，您可以找到半正弦距离的各种解决方案，例如 this one :

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。为此，请检查以下 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)