处理 DEM 数据 0 值插值的 Python 代码
以下代码基于rasterio读取 / 写入 TIFF、scipy实现插值，可对 DEM 中的 0 值（异常值）进行邻近插值或反距离加权插值（IDW）修复。
步骤说明：
- 读取 DEM 数据，提取有效数据（非 0 值）和待插值的 0 值位置；
- 可选两种插值方式：邻近插值（快速）、反距离加权插值（更平滑）；
- 插值修复 0 值区域后，保存为新的 TIFF 文件；
- 增加数据范围校验（确保插值结果在合理区间）。
import numpy as np import rasterio from scipy.interpolate import NearestNDInterpolator, RBFInterpolator from rasterio.transform import from_origin def repair_dem_zero_values(input_tif, output_tif, method="nearest", power=2): """ 修复DEM中0值（异常值）的函数 :param input_tif: 输入DEM文件路径（dem.tif） :param output_tif: 输出修复后的DEM文件路径 :param method: 插值方法，可选 "nearest"（邻近插值）或 "idw"（反距离加权） :param power: IDW插值的幂次（默认2，值越大近点权重越高） """ # 1. 读取DEM数据 with rasterio.open(input_tif) as src: dem_data = src.read(1) # 读取第一波段 profile = src.profile # 获取影像元数据（投影、分辨率、变换等） transform = src.transform # 地理变换参数 nodata = src.nodata # 原始nodata值（若有） # 2. 标记有效数据和待插值位置 # 有效数据：非0值（0为异常值），且在合理范围0-5500m valid_mask = (dem_data != 0) & (dem_data >= 0) & (dem_data <= 5500) # 待插值位置：0值且在影像范围内 interpolate_mask = dem_data == 0 # 检查是否有需要插值的点 if not np.any(interpolate_mask): print("无0值需要插值，直接保存原始数据") with rasterio.open(output_tif, "w", **profile) as dst: dst.write(dem_data, 1) return # 3. 提取有效数据的坐标和值 # 获取行列坐标（可转换为地理坐标，此处用行列坐标插值更高效） rows, cols = np.meshgrid(np.arange(dem_data.shape[0]), np.arange(dem_data.shape[1]), indexing="ij") # 有效数据的行列坐标 valid_rows = rows[valid_mask] valid_cols = cols[valid_mask] valid_values = dem_data[valid_mask] # 待插值的行列坐标 interp_rows = rows[interpolate_mask] interp_cols = cols[interpolate_mask] # 组合坐标为N×2的数组（适配scipy插值接口） valid_points = np.column_stack((valid_rows, valid_cols)) interp_points = np.column_stack((interp_rows, interp_cols)) # 4. 执行插值 if method == "nearest": # 邻近插值（快速，适合大区域） interpolator = NearestNDInterpolator(valid_points, valid_values) interp_values = interpolator(interp_points) elif method == "idw": # 反距离加权插值（更平滑，计算稍慢） # 用RBFInterpolator模拟IDW：函数类型选"inverse_multiquadric"等价IDW interpolator = RBFInterpolator( valid_points, valid_values, function="inverse_multiquadric", smoothing=0, # 无平滑，纯IDW kernel_size=None ) # 计算IDW权重（幂次power） interp_values = interpolator(interp_points, norm=power) else: raise ValueError("method仅支持 'nearest' 或 'idw'") # 5. 修正插值结果（确保在0-5500m范围内） interp_values = np.clip(interp_values, 0, 5500) # 6. 替换0值为插值结果 dem_repaired = dem_data.copy() dem_repaired[interpolate_mask] = interp_values # 7. 保存修复后的DEM with rasterio.open(output_tif, "w", **profile) as dst: dst.write(dem_repaired, 1) print(f"修复完成！输出文件：{output_tif}") print(f"插值点数：{len(interp_values)}") print(f"有效数据占比：{np.sum(valid_mask)/dem_data.size*100:.2f}%") # -------------------------- 调用示例 -------------------------- if __name__ == "__main__": # 输入输出路径（根据实际情况修改） input_dem = "dem.tif" # 原始DEM文件 output_dem = "dem_repaired.tif" # 修复后的DEM文件 # 方法1：邻近插值（快速，推荐先试用） repair_dem_zero_values(input_dem, output_dem, method="nearest") # 方法2：反距离插值（更平滑，计算稍慢） # repair_dem_zero_values(input_dem, output_dem, method="idw", power=2)依赖安装
执行代码前需安装以下库：
pip install rasterio scipy numpy关键参数说明
method：nearest：邻近插值，速度极快，结果为最近有效像素的值，适合快速修复；idw：反距离加权插值，结果更平滑，符合地形连续性，适合对精度要求高的场景。
power（仅 IDW）：反距离的幂次，默认 2，值越大，近点对插值结果的影响越大（通常取 1-3）。clip：强制插值结果在 0-5500m 范围内，避免插值出现异常值。
注意事项
- 若 DEM 中 0 值区域过大（无邻近有效数据），插值结果可能不准确，建议先检查数据有效性；
- 地理坐标插值：若需基于经纬度 / 投影坐标插值，可将
rows/cols转换为地理坐标（通过transform参数），代码中已预留地理变换接口； - 大文件优化：若 DEM 文件超大（如 GB 级），可分块处理（参考
rasterio的window功能），避免内存溢出。
验证结果
修复后可通过以下方式验证：
# 读取修复后的数据，检查0值是否被替换 with rasterio.open("dem_repaired.tif") as src: data = src.read(1) print(f"修复后0值数量：{np.sum(data == 0)}") print(f"数据范围：{np.min(data)} - {np.max(data)}")