IDW插值 反距离权重插值 python代码

代码如图：

import numpy as np
from osgeo import gdal

def writeTIF_1_int(arr, dataset, name):
    gtif_driver = gdal.GetDriverByName("GTiff")

    # 创建切出来的要存的文件（4代表4个不都按，最后一个参数为数据类型，跟原文件一致）
    out_ds = gtif_driver.Create(name, arr.shape[1], arr.shape[0], 1,
                                gdal.GDT_UInt16)
    # 获取原图的原点坐标信息
    ori_transform = dataset.GetGeoTransform()
    # 设置裁剪出来图的原点坐标
    out_ds.SetGeoTransform(ori_transform)
    # 设置SRS属性（投影信息）
    out_ds.SetProjection(dataset.GetProjection())
    # 写入目标文件
    out_ds.GetRasterBand(1).WriteArray(arr)

    # 将缓存写入磁盘
    out_ds.FlushCache()
    del out_ds
    
def func(x, y, arr):
    return arr[x, y]

def idw_interpolation_90(points, values, unknown_point, p=2):
    # 只考虑权重大于90%的点，这样就可以排除距离很远的点的干扰！
    distances = np.sqrt(np.sum((points - unknown_point) ** 2, axis=1))
    weights = 1 / distances ** p
    weights_90 = np.percentile(weights, 90)  # 90%可自行调整
    weights[weights < weights_90] = 0
    weights = weights / np.sum(weights)
    interpolated_value = np.sum(weights * values)
    return interpolated_value

def idw_interpolation(points, values, unknown_point, p=2):
    distances = np.sqrt(np.sum((points - unknown_point) ** 2, axis=1))
    weights = 1 / distances ** p
    weights = weights / np.sum(weights)
    interpolated_value = np.sum(weights * values)
    return interpolated_value

spacing = np.linspace(250, 749, 500)
X1 = np.meshgrid(spacing, spacing)
X1 = np.reshape(X1, (2, -1)).T.astype(int)

spacing = np.linspace(0, 999, 1000)
X2 = np.meshgrid(spacing, spacing)
X2 = np.reshape(X2, (2, -1)).T.astype(int)

in_ds = gdal.Open("./dem_clip.tif")
in_band1 = in_ds.GetRasterBand(1).ReadAsArray().astype("float32")
out = np.zeros((1000, 1000)).astype("float32")


idx = 0
for unknown_point in X2:
    result = idw_interpolation_90(X1, in_band1.reshape(-1), unknown_point)
    out[unknown_point[0], unknown_point[1]] = result
    idx += 1
    print(idx)

out[250:750, 250:750] = in_band1
dataset1 = gdal.Open("./dem_clip_5delta.tif")
writeTIF_1_int(out, dataset1, "./idw_output.tif")

尺寸自行修改即可。

不仅提供了普通的反距离权重插值，并且有enhance的IDW方法，仅考虑权重大于90%的点，排除距离远的点的干扰。