Google Earth Engine（GEE）深度学习入门教程-GEE预处理篇-哨兵2时间序列影像

GEE预处理篇-哨兵2时间序列影像

数据集的在GEE平台导出已经在上篇文章讲过了，本篇文章主要介绍哨兵2时间序列的影像的预处理方法，每个人的预处理方法不同，仅供参考。

本文利用Sentinel-2影像构建研究区的时间序列。研究区阴雨天气较多，云污染导致部分时相的影像覆盖率过低，无法构建完整的时间序列影像。综合考虑研究区的影像缺失情况，本文对去云后的时间序列影像进行线性插值，以补全云污染的区域，插值的时间单边窗口为30天。即缺失像素在分别向前后30天时间窗口中寻找最近高质量的观测值，采用线性插值填补当前的缺失值。

插值噪声、薄云等其他因素影响最终的分类结果，为了减轻这些影响，我们应用Savitzky-Golay滤波算法平滑和滤波插值后的时间序列。

本文将滤波的时间窗口设置25天，多项式阶数为3阶。采用时间窗口卷积的方式会导致时间序列两端出现不可避免的误差，为了消除这种影响，本文相应延长了插值时间序列和滤波时间序列的长度，并在滤波后截断多余的序列

以下是完整代码（包含测试代码）：

var getS2Image = function(BANDS,geometry,DEBUG){
     
     
    // 第一步：选择研究区，对影像数据进行去云、计算NDVI。
    var s2 = ee.ImageCollection("COPERNICUS/S2_SR"); 
    if(DEBUG === undefined){
        DEBUG = false
    }
    if(geometry===undefined){
        var geometry = ee.Geometry.Polygon(
          [[[115.81330192014848, 33.16593487478022],
            [115.81330192014848, 33.1643901294541],
            [115.81527602598344, 33.1643901294541],
            [115.81527602598344, 33.16593487478022]]], null, false)
    }
    //因为S-G滤波的开始会出现问题，所以增长时间序列 共NUM+ADD_SERIES*2个时间序列
    var ADD_SERIES = 3
    //因为需要保证插值数据的完整，所以需要对滤波的时间序列进行增长
    var interp_SERIES = 7
    var NUM = 24
    var startDate = ee.Date('2019-06-14').advance(-(ADD_SERIES)*5, 'day')
    var endDate = startDate.advance(5*(NUM+ADD_SERIES*2+interp_SERIES), 'day')
    var filtered = s2
      .filter(ee.Filter.date(startDate, endDate))
      .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))//30->70
      .filter(ee.Filter.bounds(geometry))
    // S2数据的去云函数
    function maskS2clouds(image) {
      var cloudProb = image.select('MSK_CLDPRB');
      var snowProb = image.select('MSK_SNWPRB');
      var cloud = cloudProb.lt(50);   //5-30
      var snow = snowProb.lt(5);
      var scl = image.select('SCL'); 
      var shadow = scl.eq(3); // 3 = cloud shadow
      var cirrus = scl.eq(10); // 10 = cirrus
     
      // 云的概率小于30%或云影分类
      var mask = (cloud.and(snow)).and(cirrus.neq(1)).and(shadow.neq(1));
      return image.updateMask(mask).divide(10000)
          .toFloat()       //.toFloat()  !!如果没有这一步可能会报类型错误， 计算指数也需要此强制类型转换
          .select("B.*")
          .copyProperties(image, ["system:time_start"]);
    }
    function resample(img) {
     return ee.Image(img
     .resample('bicubic')
//     .reproject({'crs':'EPSG:4326','scale': 10.0}) 
     );
    }
    // function addNDVI(image) {
    //   //!!! .toFloat()  添加其他指数时也要进行强制类型（float）转化，防止报错
    //   var ndvi = image.normalizedDifference(['B8', 'B4']).toFloat().rename('ndvi');
    //   return image.addBands(ndvi);
    // }
    var filtered = filtered.map(maskS2clouds)
                    .map(resample)
                  //.map(addNDVI)
                  .select(BANDS)
    // 第二步：设置好需要被插值的空影像，需要包含时间信息。这样方程计算出来以后我们就可以根据时间计算当时的像元值。 这里我们每隔5天就设置一个空影像。我们把空影像标记为“interpolated”，之后把真实影像和需要被插值的影像结合到一起。
    
    var n = 5;
    var totalDays = endDate.difference(startDate, 'day');
    var daysToInterpolate = ee.List.sequence(1, totalDays, n)
    //声明多波段的影像
    var emptyImage = BANDS.map(function(){return ee.Image()})
    var initImages = daysToInterpolate.map(function(day) {
      var image = ee.Image(emptyImage).rename(BANDS).toFloat().set({
        'system:index': ee.Number(day).format('%d'),
        'system:time_start': startDate.advance(day, 'day').millis(),
        // Set a property so we can identify interpolated images
        'type': 'interpolated'
      })
      return image
    })
    
    var initCol = ee.ImageCollection.fromImages(initImages)

    var mergedCol = filtered.merge(initCol)
    var mergedCol = mergedCol.map(function(image) {
      var timeImage = image.metadata('system:time_start').rename('timestamp')
      var timeImageMasked = timeImage.updateMask(image.mask().select(0))
      return image.addBands(timeImageMasked)
    })
    
    // 第三步：利用Join函数把每个影像对应的之前和之后的影像链接起来。
    // 线性插值的间隔
    var days = 45
    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: mergedCol,
      secondary: mergedCol,
      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 filtered = join2Result.filter(ee.Filter.eq('type', 'interpolated'))

    function interpolateImages(image) {
      var image = ee.Image(image)
    
      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()
    
      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'),
      })
    
      var interpolated = beforeMosaic
        .add((afterMosaic.subtract(beforeMosaic).multiply(timeRatio)))
      var result = image.unmask(interpolated)
      return result.copyProperties(image, ['system:time_start'])
    }
    var interpolatedCol = ee.ImageCollection(
      filtered.map(interpolateImages)).select(BANDS)
    var visualization = {
        min:[2955/10000,1223/10000,1461/10000],//%2线性拉伸
        max:[6008/10000,2424/10000,2416/10000],//%2线性拉伸
        bands: ['B8', 'B4', 'B3'],
      };
    // //调试代码
    // var smoothed_ = interpolatedCol.toList(100)
    // for (var index = 0;index< NUM+ADD_SERIES*2; index++){
    //     if(DEBUG)
    //         Map.addLayer(ee.Image(smoothed_.get(index)), visualization, index.toString());
    // }
   //...
    // print(filtered)
    // var debug =  ee.ImageCollection(filtered).filter(ee.Filter.eq("system:index",'2_106')).first()
    // print(debug)
    // print(ee.Image(interpolateImages(debug)).select(BANDS))
    // Map.addLayer(ee.Image(interpolateImages(debug)).select(BANDS), {}, 'interpolateImages')
    // print('Interpolated Collection', interpolatedCol)
    //...
    // 第五步：把插出来的数据进行SG滤波。
    
    var oeel=require('users/OEEL/lib:loadAll');
    // https://www.open-geocomputing.org/OpenEarthEngineLibrary/#.ImageCollection.SavatskyGolayFilter
    // S-G滤波的时间窗口 
    
    var window = ADD_SERIES*5    // 5*5= 25
    var millis = ee.Number(window).multiply(1000*60*60*24)
    
    // Use the same maxDiffFilter we used earlier
    var maxDiffFilter = ee.Filter.maxDifference({
      difference: millis,
      leftField: 'system:time_start',
      rightField: 'system:time_start'
    })
    
    // Use the default distanceFunction
    var distanceFunction = function(infromedImage, estimationImage) {
      return ee.Image.constant(
          ee.Number(infromedImage.get('system:time_start'))
          .subtract(
            ee.Number(estimationImage.get('system:time_start')))
            );
      }
    
    // 把SG滤波函数设置为三次函数
    var order = 3;
    var smoothed = oeel.ImageCollection.SavatskyGolayFilter(
      interpolatedCol, 
      maxDiffFilter,
      distanceFunction,
      order)
    //选择d_0_波段并重命名
    var d_0_BANDS = BANDS.map(function(band){return "d_0_"+band})
    var smoothed = smoothed.select(d_0_BANDS,BANDS)
    var allSmoothed = smoothed
    //截取掉延长后的序列
    
    smoothed = ee.ImageCollection(smoothed.toList(NUM+ADD_SERIES*2).slice(ADD_SERIES,ADD_SERIES+NUM))
    
    //调整编号
    smoothed = smoothed.map(function(img){return img.set("index",smoothed.toList(NUM).indexOf(img))})
    //合成缩减数据
    var visualization = {
        min:[2955/10000,1223/10000,1461/10000],//%2线性拉伸
        max:[6008/10000,2424/10000,2416/10000],//%2线性拉伸
        bands: ['B8', 'B4', 'B3'],
      };

    smoothed = smoothed.toList(NUM)
    var smoothedmean = ee.List([]);
    for (var index = 0;index< NUM/2; index++){
        var temp = ee.ImageCollection([
                  smoothed.get(2*index),
                  smoothed.get(2*index+1)])
                  .mean();
        smoothedmean=smoothedmean.add(temp);
        if(DEBUG)
            Map.addLayer(temp, visualization, index.toString());
    }
    smoothedmean = ee.ImageCollection(smoothedmean)
    //如果想返回滤波前的线性插值时间序列，DEBUG设置为true
    if(DEBUG)
        return [smoothedmean,smoothed,interpolatedCol,allSmoothed]
    else
        return smoothedmean
}
exports.getS2Image = getS2Image






  // 显示出来并绘图
var geometry = ee.Geometry.Polygon(
          [[[115.81330192014848, 33.16593487478022],
            [115.81330192014848, 33.1643901294541],
            [115.81527602598344, 33.1643901294541],
            [115.81527602598344, 33.16593487478022]]], null, false)
// var geometry =geometry2
var S2 = getS2Image(["B4","B3","B8",],geometry,true)
var S2_SG = ee.ImageCollection(S2[1]).select(["B8"], ["SG"])
var S2_interp = ee.ImageCollection(S2[2]).select(["B8"], ["interp"])
print(S2_SG)
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)
// 把原始的和经过SG滤波的数据NDVI曲线画出来
var chart = ui.Chart.image.series({
  imageCollection: S2_interp.merge(S2_SG),
  region: geometry,
  reducer: ee.Reducer.mean(),
  scale: 10,
}).setOptions({
      lineWidth: 1,
      title: "滤波后的时序曲线",
      interpolateNulls: true,
      vAxis: {title: 'B8'},
      hAxis: {title: '', format: 'YYYY-MMM'},
      series: {
        0: {color: 'red', lineWidth: 2 }, // Smoothed NDVI
        1: {color: 'blue', lineWidth: 1, 
          lineDashStyle: [1, 1], pointSize: 1,
          }, // interpolatedCol NDVI
      },

    })
print(chart)