去除遥感影像黑边

ytkz2025-02-182025-02-18

为什么会出现黑边

某些遥感L1级别的影像本身自带黑边，所以在经过预处理后，L2级别的遥感影像存在黑边。

L1级别的影像的边缘如下，黑边的像素值是0。

L1级别的透明区域的像素值为Nodata，如下。

我们可以在arcgis上设置L1级别的像素为透明，但是这个设置是不改变影像，只是视觉效果上发生改变，操作如下：右击影像，选择属性，然后在Symbology勾选上Display Background Value按钮

可视化结果如下：

然而，经过6S大气校正、RPC几何校正后的影像边缘如下，黑边的像素值是负值。

怎么去除

逻辑简单，将L2级别的影像中的负值设置为0即可。

代码简单，如下：

from osgeo import gdal
import numpy as np
import os

def remove_black_edges_border_only(input_path, output_path, block_size=2000, edge_width=100):
    """
    以分块方式读取遥感影像，仅对影像边缘区域（图像边界内 edge_width 像素范围）进行检测，
    将小于0（负值）的像素值置为0；影像内部区域不作检测修改。

    参数：
      input_path  : 输入影像路径（例如 GeoTIFF 文件）
      output_path : 输出影像路径
      block_size  : 分块尺寸（默认2000，可以根据内存情况调整）
      edge_width  : 边缘区域宽度（单位：像素），仅对图像四周 edge_width 范围内的像素进行处理
    """
    # 打开输入影像（只读模式）
    ds = gdal.Open(input_path, gdal.GA_ReadOnly)
    if ds is None:
        raise FileNotFoundError(f"无法打开输入文件：{input_path}")

    cols = ds.RasterXSize
    rows = ds.RasterYSize
    bands = ds.RasterCount
    geotransform = ds.GetGeoTransform()
    projection = ds.GetProjection()

    # 创建输出影像，与输入影像在尺寸、波段和数据类型上保持一致
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(output_path, cols, rows, bands, ds.GetRasterBand(1).DataType)
    out_ds.SetGeoTransform(geotransform)
    out_ds.SetProjection(projection)

    # 针对每个波段分块处理
    for i in range(1, bands + 1):
        in_band = ds.GetRasterBand(i)
        out_band = out_ds.GetRasterBand(i)

        # 按块遍历整幅影像
        for y in range(0, rows, block_size):
            block_height = block_size if y + block_size < rows else rows - y
            for x in range(0, cols, block_size):
                block_width = block_size if x + block_size < cols else cols - x

                # 读取当前块数据
                data = in_band.ReadAsArray(x, y, block_width, block_height)
                # 保证数据为浮点型，便于比较（也可根据实际数据类型调整）
                data = data.astype(np.float32)

                # 生成当前块的全局像素坐标网格
                xs = np.arange(x, x + block_width)
                ys = np.arange(y, y + block_height)
                # X: 每行的列坐标； Y: 每列的行坐标，形状均为 (block_height, block_width)
                X, Y = np.meshgrid(xs, ys)
                # 定义边缘区域：全局坐标位于图像左边 edge_width 内、右边 edge_width 内、
                # 上边 edge_width 内或下边 edge_width 内的像素
                border_mask = ((X < edge_width) | (X >= cols - edge_width) |
                               (Y < edge_width) | (Y >= rows - edge_width))
                # 仅对边缘区域中小于0的像素值设置为0
                data[border_mask & (data < 0)] = 0

                # 将当前块处理后的数据写回输出影像
                out_band.WriteArray(data, xoff=x, yoff=y)
            # 刷新缓存
            out_band.FlushCache()

    # 关闭数据集
    ds = None
    out_ds = None
    print(f"分块处理完成，输出文件保存在：{output_path}")


# 示例调用
if __name__ == "__main__":
    input_file = r"D:\orth.tif"  # 替换为实际路径
    output_file = os.path.splitext(input_file)[0] + '_remove_blacks_border.tif'
    remove_black_edges_border_only(input_file, output_file, block_size=2000, edge_width=100)