通过指定的标志掩码填充栅格的孔洞

from osgeo import gdal, gdal_array
import numpy as np

def fill_holes(raster_path, mask_path, output_path):
    # 打开栅格数据和掩码数据
    raster = gdal.Open(raster_path)
    mask = gdal.Open(mask_path)
    
    # 读取数据和掩码
    band = raster.GetRasterBand(1)
    data = band.ReadAsArray()
    mask_band = mask.GetRasterBand(1)
    mask_data = mask_band.ReadAsArray()
    
    # 应用掩码填充孔洞
    filled_data = np.where(mask_data == 1, data, np.nan)
    filled_data = np.nan_to_num(filled_data, nan=np.median(data))
    
    # 创建输出文件
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(output_path, raster.RasterXSize, raster.RasterYSize, 1, band.DataType)
    out_ds.SetGeoTransform(raster.GetGeoTransform())
    out_ds.SetProjection(raster.GetProjection())
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(filled_data)
    
    # 清理
    out_ds.FlushCache()
    raster = None
    mask = None
    out_ds = None

# 使用示例
fill_holes('input_raster.tif', 'mask.tif', 'output_filled_raster.tif')