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

柔光的暖阳◎ 2024-05-28 10:30 47阅读 0赞

//输入目标区域
    
    
    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())
      }
    })