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'})
}

结果展示：