【GEE】利用哨兵2号计算NDBI并使用前后时段影像补全空值

//输入目标区域


var geometry = table;
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var filtered = s2
  .filter(ee.Filter.date('2019-07-01', '2019-08-31'))//修改数据时间范围
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
  .filter(ee.Filter.bounds(geometry))
// Write a function for Cloud masking
function maskCloudAndShadowsSR(image) {
   var qa = image.select('QA60');
  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000).set(image.toDictionary(image.propertyNames()))
      .select("B.*")
      .copyProperties(image, ["system:time_start"]);
}
var filtered = filtered.map(maskCloudAndShadowsSR)
//在有效像素下添加时间信息，便于后续利用
var filtered = filtered.map(function(image) {
  var timeImage = image.metadata('system:time_start').rename('timestamp')
  var timeImageMasked = timeImage.updateMask(image.mask().select(0))
  return image.addBands(timeImageMasked)
})

var days = 30 //前后时段长
var millis = ee.Number(days).multiply(1000*60*60*24)

var maxDiffFilter = ee.Filter.maxDifference({
  difference: millis,
  leftField: 'system:time_start',
  rightField: 'system:time_start'
})


var lessEqFilter = ee.Filter.lessThanOrEquals({
  leftField: 'system:time_start',
  rightField: 'system:time_start'
})
var greaterEqFilter = ee.Filter.greaterThanOrEquals({
  leftField: 'system:time_start',
  rightField: 'system:time_start'
})

var filter1 = ee.Filter.and(maxDiffFilter, lessEqFilter)
var join1 = ee.Join.saveAll({
  matchesKey: 'after',
  ordering: 'system:time_start',
  ascending: false})
var join1Result = join1.apply({
  primary: filtered,
  secondary: filtered,
  condition: filter1
})
var filter2 = ee.Filter.and(maxDiffFilter, greaterEqFilter)
var join2 = ee.Join.saveAll({
  matchesKey: 'before',
  ordering: 'system:time_start',
  ascending: true})
var join2Result = join2.apply({
  primary: join1Result,
  secondary: join1Result,
  condition: filter2
})
//插值
var interpolateImages = function(image) {
  var image = ee.Image(image)
  // We get the list of before and after images from the image property
  // Mosaic the images so we a before and after image with the closest unmasked pixel
  var beforeImages = ee.List(image.get('before'))
  var beforeMosaic = ee.ImageCollection.fromImages(beforeImages).mosaic()
  var afterImages = ee.List(image.get('after'))
  var afterMosaic = ee.ImageCollection.fromImages(afterImages).mosaic()
  // Get image with before and after times
  var t1 = beforeMosaic.select('timestamp').rename('t1')
  var t2 = afterMosaic.select('timestamp').rename('t2')
  var t = image.metadata('system:time_start').rename('t')
  var timeImage = ee.Image.cat([t1, t2, t])
  var timeRatio = timeImage.expression('(t - t1) / (t2 - t1)', {
    't': timeImage.select('t'),
    't1': timeImage.select('t1'),
    't2': timeImage.select('t2'),
  })
  // Compute an image with the interpolated image y
  var interpolated = beforeMosaic
    .add((afterMosaic.subtract(beforeMosaic).multiply(timeRatio)))
  // Replace the masked pixels in the current image with the average value
  var result = image.unmask(interpolated)
  return result.copyProperties(image, ['system:time_start'])
}


// function addNDVI(image){
//   var ndvi = image.normalizedDifference(['B8','B4']).rename('ndvi');
//   var doy = image.date().format("yyyyMMdd")
//   return image.addBands(ndvi).select('ndvi').clip(geometry)
//               .set({'Date':doy});
// }
function addNDBI(image){
// 方法 3 自己手动加减乘除
  var ndbi = image.select('B11').subtract(image.select('B8')).divide(image.select('B11').add(image.select('B8')))
                  .float().rename('NDBI');
  var doy = image.date().format("yyyyMMdd")
  return image.addBands(ndbi).select('NDBI').clip(geometry)
              .set({'Date':doy});
}
var interpolatedCol = ee.ImageCollection(
  join2Result.map(interpolateImages)).map(addNDBI)
print(interpolatedCol)
//数据输出函数
function exportImage(image,roi,fileName){
  Export.image.toDrive({
    image: image,
    description: fileName + '_NDBI',
    folder: 'NDBI',
    scale: 10, //分辨率
    region: roi,
    maxPixels:1e13,
    crs:'EPSG:4326',  //坐标系
    fileFormat: 'GeoTIFF',
  });
}
//生成列表，迭代下载
var indexList = interpolatedCol.reduceColumns(ee.Reducer.toList(),["Date"]).get("list");
print("indexList",indexList);
indexList.evaluate(function(indexs){
  var indexs1 = ee.List(indexs).distinct();
  var length = indexs1.length().getInfo();
  print(indexs1)
  for (var i=0;i<length;i++){
    var img =interpolatedCol.filter(ee.Filter.eq("Date",indexs1.get(i))).mean()
    //Map.addLayer(img,indexs1,get(i).getInfo(),false);
    exportImage(img,geometry,indexs1.get(i).getInfo())
  }
})