GEE（8）：使用MODIS填补由去云后的Landsat影像计算得到的NDVI数据

最近想要在GEE中使用Landsat影像计算一下广州的NDVI值，发现这片区域云覆盖较多，去云以后部分月份的数据很少，就造成NDVI计算结果缺失的问题。经过查阅相关资料，可以使用MODIS的NDVI产品来填补由Landsat计算NDVI数据缺失的部分。主要用到了GEE中的.blend()函数。

.

实现代码如下：

//设置研究区
var geometry = 
    /* color: #d63000 */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[112.38214918009086, 23.76360001525737],
          [112.38214918009086, 22.348421631032746],
          [114.09601636759086, 22.348421631032746],
          [114.09601636759086, 23.76360001525737]]], null, false);

//Landsat去云函数
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}


// 定义计算NDVI函数
var getNDVI = function(image){
var NDVI = image.normalizedDifference(['B5', 'B4']).rename('NDVI')
return image.addBands(NDVI).copyProperties(image, ["system:time_start"])
};

//由Landsat 8 计算NDVI
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.map(maskL8sr)
.filterBounds(geometry)
.select('B4','B5')
.map(getNDVI)
.select(['NDVI'])

//导入MODIS NDVI产品
var MOD13Q1 = ee.ImageCollection('MODIS/006/MOD13Q1').select('NDVI')
          .map(function(image){
              var img = image.multiply(0.0001);
              return img.set('system:time_start',image.get('system:time_start')) ;
           });


// 以2020年为例，计算2020年12个月NDVI最大值
var listMonths = ee.List.sequence(1,12);
var collectMonth = ee.ImageCollection(listMonths
  .map(function(month) {
  var start = ee.Date.fromYMD(2020, month, 1);
  var end = ee.Date.fromYMD(2020, month, 1).advance(1, 'month');
  var landsat = l8.filterDate(start, end).max()
  var modis = MOD13Q1.filterDate(start, end).max().reproject('EPSG:4326',null,250)
    return modis.blend(landsat)
                  .reduce(ee.Reducer.max()).float().clip(geometry)
                  .set('system:time_start',month)
                  .resample();
}));
//print (collectMonth);

var visParams = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};

var ndviVis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};


var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2020-01-01', '2020-02-01')
                  .map(maskL8sr);
var ndvil8 = dataset.map(getNDVI)
                    .select('NDVI')
                    .max()
                    .clip(geometry);
              

Map.addLayer(ndvil8,ndviVis,'landsat_NDVI'); //使用去云后Landsat数据计算得到的NDVI
Map.addLayer(ee.Image(collectMonth.toList(12).get(0)),ndviVis,'Tianbu'); //使用MODIS填补后
Map.addLayer(dataset,visParams,'OriginImage'); //原始影像

//批量展示2020年12个月的NDVI计算结果
for(var i = 0;i < 12;i++){
  var image = ee.Image(collectMonth.toList(12).get(i));
  Map.addLayer(image,ndviVis,'NDVI_'+String(i));
}

下面分别为填补前后的结果：

原始影像：

去云后计算的NDVI：

填补后的NDVI：

小结：

这只是一种比较简单的填补方法，由于文中所选取的月份云量较大，所以填补内容较多，效果还是可以的。在实际中可以结合sentinel 数据计算NDVI来替换Landsat中云量较大的时间段。在之前也见过一种方法，是通过计算哨兵数据两个月份（比如1月和2月）之间NDVI值的比例，如果现在需要使用Landsat计算2月的NDVI，就将这个比例乘以1月的Landsat计算的NDVI，从而得到2月的NDVI值，这也是一种替换方法。大家实际计算时可以自行选择。