栅格篇 - 去除影像黑边

@ 20240723 & lth
目标：从数据层面去除栅格影像在arcmap中显示时出现的黑边问题
逻辑：黑边的出现其实就是栅格无效像元值的显示

前置

GEE-哨兵2栅格影像下载

// 自己绘制roi或传入坐标或传入文件都可以，我这里选择绘制，然后下载对应的哨兵2 影像
function exportS2ToDrive(roi, startDate, endDate, cloudThresh, folderName) {
  // 统一 roi 为 Geometry
  roi = roi instanceof ee.FeatureCollection ? roi.geometry() :
        roi instanceof ee.Feature           ? roi.geometry() :
        roi;

  // 查询影像
  var s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
             .filterBounds(roi)
             .filterDate(startDate, endDate)
             .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', cloudThresh))
             .select(['B2','B3','B4','B8']);

  // 预览
  Map.centerObject(roi, 10);
  Map.addLayer(s2.mosaic().clip(roi),
               {bands:['B4','B3','B2'], min:0, max:3000}, 'RGB');

  // 创建任务
  var imgList = s2.toList(s2.size());
  var n = imgList.size().getInfo();          // 取客户端数量
  for (var i = 0; i < n; i++) {
    var img  = ee.Image(imgList.get(i));
    var date = img.date().format('yyyyMMdd').getInfo();
    var id   = img.id().getInfo().split('/').pop();
    var desc = 'S2_' + date + '_' + id;

    Export.image.toDrive({
      image: img.clip(roi),
      description: desc,
      folder: folderName,
      region: roi,
      scale: 10,
      crs: 'EPSG:4490',
      maxPixels: 1e13
    });
  }
}

/************ 示例调用 ************/
// var roi = /* color: #d63000 */ee.Geometry.Polygon(
//   [[[116.3, 39.9], [116.4, 39.9], [116.4, 40.0], [116.3, 40.0]]]);

exportS2ToDrive(roi, '2025-05-01', '2025-05-31', 50, 'S2_May2025');

去除影像黑边

实现代码

from osgeo import gdal

def remove_black_borders(input_path:str, output_path:str, no_data_value:int=0):
    """

    使用GDAL设置NoData值来去除黑边。

    :param input_path: 输入影像路径

    :param output_path: 输出影像路径

    :param no_data_value: 黑边像素的值，默认为0

    """

    try:

        # 打开输入影像

        dataset = gdal.Open(input_path, gdal.GA_ReadOnly)

        if dataset is None:

            raise Exception(f"无法打开输入文件 {input_path}")

        # 设置Translate选项，指定NoData值 + 分块 + 压缩

        translate_options = gdal.TranslateOptions(

            noData=no_data_value,  # 注意这里是noData而非nodata

            format='GTiff',

            creationOptions=[

                'BLOCKXSIZE=512', 'BLOCKYSIZE=512',

                'COMPRESS=LZW', 'TILED=YES',

                'BIGTIFF=IF_SAFER'  # 支持大于4GB的TIFF文件

            ],

            callback=gdal.TermProgress_nocb  # 显示进度条（可选）

        )

        # 使用Translate函数转换并输出新的影像

        gdal.Translate(output_path, dataset, options=translate_options)

        print(f"✅ 成功设置NoData值，已保存至 {output_path}")

    except Exception as e:

        print(f"❌ 处理失败: {e}")

    finally:

        dataset = None  # 确保释放资源