基于GEE的Landsat 7ETM+条带填补

项目简介

该项目使用Google Earth Engine (GEE)平台，对Landsat 5和Landsat 7卫星影像进行预处理与影像填补操作。主要功能包括影像的选取、波段处理、缺失影像的填补以及最终影像的导出。代码中的核心功能是通过空间回归方法对Landsat 7和Landsat 5影像进行时序配准与插值填补。

主要步骤

1. 定义感兴趣区域 (ROI)：
    - 使用四个坐标点绘制矩形多边形，确定需要处理的地理区域。

2. 选择影像数据集：
    - 使用Landsat 5和Landsat 7数据集 (`ee.ImageCollection`)，分别选择时间范围为2010年5月1日至2010年5月31日之间的影像。

3. 影像预处理：
    - 定义影像预处理函数，将选择的波段进行裁剪、波段重命名，并转换像元值范围，以便与其他数据进行处理和分析。

4. 影像填补：
    - 使用空间回归方法，计算两个影像集之间的关系，通过滑动窗口对像元进行回归分析并填补缺失区域。主要步骤包括：
      - 选择匹配的填补影像；
      - 计算回归模型，生成比例和偏移量；
      - 处理填补失败的情况，使用默认值填补。

5. 显示与导出结果:
    - 显示原始Landsat 7、Landsat 5影像和填补后的Landsat 5影像，并将结果导出为影像文件，便于后续分析。

核心代码说明

- 预处理函数 `preprocessImage`:
  该函数处理Landsat影像集中的波段并应用裁剪区域。通过乘以系数与增加偏移量将影像值缩放至0-1范围。

- 影像填补函数 `GapFill`:
  这是实现影像插值填补的主要函数，输入Landsat 7影像集和Landsat 5影像集，并对它们进行空间回归分析，最终生成无缺失的填补影像。

- 导出函数 `Export.image.toDrive`:
  使用GEE的导出功能将填补后的影像导出为文件，以便本地存储和处理。

如何运行

1. 打开Google Earth Engine的代码编辑器，并将上述代码粘贴到编辑器中。
2. 修改`startDate`和`endDate`参数以更改需要处理的影像时间范围。
3. 运行代码，预处理并显示影像，同时将填补后的影像导出为本地文件。

var roi = 
    /* color: #b50000 */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[117.29699355756334, 40.12436832412577],
          [117.29699355756334, 39.890323465395774],
          [117.65610915814928, 39.890323465395774],
          [117.65610915814928, 40.12436832412577]]], null, false);

// 定义数据集
var landsat7 = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2');
var landsat5 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2');
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2');

// 定义波段
var bands = ['red', 'green', 'blue', 'nir', 'swir'];
var percentile = 30;
var imageParams = {min: 0.0, max: 0.3, bands: ['red', 'green', 'blue']};

// 输入时间范围
var startDate = '2010-05-01';
var endDate = '2010-05-31';

// 影像预处理函数
var preprocessImage = function(image) {
  var percentileImage = image
    .select(['SR_B3', 'SR_B2', 'SR_B1', 'SR_B4', 'SR_B5'])
    .multiply(2.75e-05).add(-0.2)
    .clip(roi)
    .rename(bands);  // 重命名波段，使其保持与原始影像一致
  return percentileImage
  //.updateMask(image.select(0).mask().focal_min(90, 'circle', 'meters'))
  .copyProperties(image, ['system:time_start']);  // 复制时间等属性;
};

// 选择 Landsat 7 影像
var landsat7Collection = landsat7
  //.select(['B3', 'B2', 'B1', 'B4', 'B5'], bands)
  .filterDate(startDate, endDate)
  .map(preprocessImage);
print('landsat7Collection:',landsat7Collection.limit(10))


// 选择 Landsat 5 影像
var landsat5Collection = landsat5
  //.select(['SR_B3', 'SR_B2', 'SR_B1', 'SR_B4', 'SR_B5'], bands)
  .filterDate(startDate, endDate)
  //.reduce(ee.Reducer.percentile([percentile]))
  .map(preprocessImage);
print('landsat5Collection:',landsat5Collection.limit(10))

/*// 选择 Landsat 8 影像
var landsat8Collection = landsat8
  .select(['SR_B4', 'SR_B3', 'SR_B2', 'SR_B5', 'SR_B6'], bands)
  .filterBounds(roi)
  .filterDate(startDate, endDate)
  .reduce(ee.Reducer.percentile([percentile]))
  .rename(bands)
  .map(preprocessImage);
*/
// 填补函数
var GapFill = function(srcCollection, fillCollection, kernelSize) {
  var MIN_SCALE = 1/3;
  var MAX_SCALE = 3;
  var MIN_NEIGHBORS = 144;
  var kernel = ee.Kernel.square(kernelSize * 30, "meters", false);

  var fillFunction = function(image) {
    var date = ee.Date(image.get('system:time_start'));

    // 筛选匹配的填补影像
    var fillImage = fillCollection
      .filterBounds(roi)
      .filterDate(date.advance(-10, 'day'), date.advance(10, 'day'))
      .median()
      .select(bands);

    var srcImage = srcCollection
      .filterBounds(roi)
      .filterDate(date.advance(-10, 'day'), date.advance(10, 'day'))
      .median()
      .select(bands);
      
    // 计算共同掩膜
    var commonMask = srcImage.mask().and(fillImage.mask());
    var fillMasked = fillImage.updateMask(commonMask);
    var srcMasked = srcImage.updateMask(commonMask);

    // 进行回归计算
    var regression = fillMasked.addBands(srcMasked);
    regression = regression.select(regression.bandNames().sort());
    var fit = regression.reduceNeighborhood(ee.Reducer.linearFit().forEach(srcImage.bandNames()), kernel, null, false);
    var offset = fit.select(".*_offset");
    var scale = fit.select(".*_scale");

    // 计算二级比例因子
    var reducer = ee.Reducer.mean().combine(ee.Reducer.stdDev(), null, true);
    var srcStats = srcImage.reduceNeighborhood(reducer, kernel, null, false);
    var fillStats = fillImage.reduceNeighborhood(reducer, kernel, null, false);
    var scaleFactor = srcStats.select(".*stdDev").divide(fillStats.select(".*stdDev"));
    var offsetFactor = srcStats.select(".*mean").subtract(fillStats.select(".*mean").multiply(scaleFactor));
    var invalidScale = scale.lt(MIN_SCALE).or(scale.gt(MAX_SCALE));
    scale = scale.where(invalidScale, scaleFactor);
    offset = offset.where(invalidScale, offsetFactor);

    // 处理所有其他方法都失败的情况
    var fallbackScale = scale.lt(MIN_SCALE).or(scale.gt(MAX_SCALE));
    scale = scale.where(fallbackScale, 1);
    offset = offset.where(fallbackScale, srcStats.select(".*mean").subtract(fillStats.select(".*mean")));
    var neighborCount = commonMask.reduceNeighborhood(ee.Reducer.count(), kernel, null, true, "boxcar");
    var scaledImage = fillImage.multiply(scale).add(offset).updateMask(neighborCount.gte(MIN_NEIGHBORS));

    return image.unmask(scaledImage, true).copyProperties(image, ['system:time_start']);
  };

  return srcCollection.map(fillFunction);
};

// 应用填补函数
var filledLandsat5 = GapFill(landsat7Collection, landsat5Collection, 10);
print('filledLandsat5:',filledLandsat5.limit(10))

//var filledLandsat8 = GapFill(landsat7Collection, landsat8Collection, 10);

// 组合填补结果
//var resultCollection = ee.ImageCollection.fromImages([filledLandsat5, filledLandsat8]);
//var finalResult = resultCollection.mosaic();

// 显示结果
Map.centerObject(roi, 12); // 调整地图中心和缩放级别
Map.addLayer(landsat7Collection.mean(), imageParams, "Landsat 7");
Map.addLayer(landsat5Collection.mean(), imageParams, "Landsat 5");
Map.addLayer(filledLandsat5.mean(), imageParams, "Landsat 5 Filled");
//Map.addLayer(filledLandsat8, imageParams, "Landsat 8 Filled");
//Map.addLayer(finalResult, imageParams, "Combined Filled");

// 导出处理后的影像
Export.image.toDrive({
  image: filledLandsat5.mean(), // 选择Landsat 5填补结果进行导出
  description: 'landsat7_filled_landsat5',
  folder: 'landsat7',
  scale: 30,
});