【数据优化】基于GEE填补遥感缺失数据

1.写在前面

  在遥感影像分析中，我们经常会遇到由于云层遮挡、传感器故障等多重因素导致的图像数据缺失问题。为了解决这一挑战，常用的技术包括利用一年内数据的均值或最小值进行填充，以及采用线性插值等方法。在本文中，我们将探索如何借助 Google Earth Engine (GEE) 这一强大工具，以简洁高效的方式实现这些数据填充技术。这里我先使用年内数据填充法对多年数据进行填充。
  天地图平台发布了带有审图号（审图号：GS（2024）0650号）的最新的中国省、市、县行政区划可视化。该数据已经上传，并开通共享公开方便大家正常使用，代码如下：

2.填充代码

研究区设置：

// GS（2024）0650号========================================================================================
var china_provinces = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_provinces");
var china_city = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_city");
var china_county = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_county");
var sichuan = china_provinces.filter(ee.Filter.eq('name','四川省'))
var chengdu = china_city.filter(ee.Filter.eq('name','成都市'))
var Jingtang = china_county.filter(ee.Filter.eq('name','金堂县'))

Map.addLayer(sichuan.style({fillColor:'00000000',color:'red'}),{},"四川省", 0) //ff0000
Map.addLayer(chengdu.style({fillColor:'00000000',color:'blue'}),{},"成都市", 1) //ffff00
Map.addLayer(Jingtang.style({color:"black"}),{},"金堂县")
Map.centerObject(Jingtang, 10);

2.1 年内中值数据填充MODIS NPP空值

function FillgapNPP(image){
  var median = image.reduceRegion({
    reducer:ee.Reducer.median(),
    geometry:Jingtang.geometry(),
    scale:500,
    maxPixels:1e13
  }).values();
  median = ee.Number(median);
  var FillImage = image.unmask(median).clip(Jingtang);
  return FillImage;
}

var NPP = ee.ImageCollection("MODIS/061/MOD17A3HGF")
var NPPdata = NPP.map(function(image){return image.clip(Jingtang)})
                 .toList(22)
                 .aside(print);

var NPP2011 = ee.Image(NPPdata.get(10)).select('Npp').aside(print);
var NPP2016 = ee.Image(NPPdata.get(15)).select('Npp').aside(print);
var NPP2022 = ee.Image(NPPdata.get(21)).select('Npp').aside(print);

var NPP2011_Fill = FillgapNPP(NPP2011);
var NPP2016_Fill = FillgapNPP(NPP2016);
var NPP2022_Fill = FillgapNPP(NPP2022);

var visualization = {
  bands: ['Npp'],
  min: 0,
  max: 19000,
  palette: ['ffffff', 'ce7e45', 'df923d', 'f1b555', 'fcd163', 
  '99b718', '74a901','66a000', '529400', '3e8601', '207401', 
  '056201', '004c00', '023b01','012e01', '011d01', '011301'
    ]
};

Map.addLayer(NPP2022,visualization,'NPP2022')
Map.addLayer(NPP2022_Fill,visualization,'NPP2022_Fill')

结果展示：

2.2 年内中值数据填充Landsat8 NDVI空值

  这里为了让数据出现空值，我使用了2014年-2023年6-8月的数据，若使用全年的数据，可能不会出现空值的情况。

//2.Landsat8 NDVI插值======================================================================
//函数定义=================================================================================
// Applies scaling factors.
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}
function rmCloudNew(image) {
  var cloudShadowBitMask = (1 << 4); 
  var cloudsBitMask = (1 << 3); 
  var qa = image.select('QA_PIXEL'); 
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) 
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0)); 
  return image.updateMask(mask)
              .copyProperties(image)
              .copyProperties(image, ["system:time_start"]);
}
var get_NDVI = function(image) {
  var NDVI=image.normalizedDifference(['nir','red']).rename(['NDVI']);
  image=image.addBands(NDVI)
  return image.select("NDVI")
};
// 定义一个函数来填充 Landsat 8 图像中的缺失值
function fillGapLandsat8(image) {
  // 计算整个 Jingtang 区域内所有像素值的中位数
  var median = image.reduceRegion({
    reducer: ee.Reducer.median(),
    geometry: Jingtang.geometry(),
    scale: 30,
    maxPixels: 1e13
  }).values(); 
  median = ee.Number(median);
  var fillImage = image.unmask(median).clip(Jingtang);
  return fillImage;
}


// Landsat8===================================================================================
var Landsat = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                .filterBounds(Jingtang)
                .filter(ee.Filter.calendarRange(2014,2023,'year'))
                .filter(ee.Filter.calendarRange(1,12,'month'))
                .map(applyScaleFactors)
                .select(['SR_B4','SR_B5','QA_PIXEL'],['red','nir','QA_PIXEL'])
                .map(rmCloudNew)
                .map(get_NDVI);
print("Landsat:", Landsat)

var startDate = 2014
var endDate = 2023
for(var i = startDate;i<=endDate;i++){
  var ndvi_year = Landsat.filterDate(i+'-06-01', i+'-08-31').select('NDVI')
  var ndvi_median = ndvi_year.median().clip(Jingtang)
  
  // 给每个月的NDVI图像指定一个波段名称
  var ndvi_band = ndvi_median.rename('NDVI_median_' + i);
    
  var Landsat8_fill = fillGapLandsat8(ndvi_band)
 
  Map.addLayer(ndvi_band, colorizedVis, 'NDVI_median_' + i, 0);
  Map.addLayer(Landsat8_fill, colorizedVis, 'NDVI_median_' + i + '_Fill', 0);

  Export.image.toDrive({
      image: Landsat8_fill,
      description: i+'year_median',
      region: Jingtang,
      scale: 30,
      maxPixels: 1e13,
      folder: 'NDVI_year'
    })
}

结果展示：