GEE Savitzky-Golay （SG滤波）

2020-11-01

9045

原创

### 数据预览
[![](http://www.seibert.cn/static/static/images/uploadimages/1_1593098286305.jpg)](http://www.seibert.cn/static/static/images/uploadimages/1_1593098286305.jpg)

### 源码分析

/**
     *  Savitzky-Golay smoothing for the time-series images
     *
     * 
     * 
     * @Source
     * https://www.seibert.cn/
     * 
     * @Author Saibo Li, IGSNRR, CAS
     * 
     */ 
	var geometry = 
        /* color: #d63000 */
        /* displayProperties: [
          {
            "type": "rectangle"
          }
        ] */
        ee.Geometry.Polygon(
            [[[113.98422435778511, 42.00616321967621],
              [113.98422435778511, 38.24435054695926],
              [117.41195873278511, 38.24435054695926],
              [117.41195873278511, 42.00616321967621]]], null, false);
    var dataset = ee.ImageCollection("MODIS/006/MCD15A3H");
    /*
    Preprocess data
    **/
    // set time
    var start_date = '2010-01-01'
    var end_date = '2010-12-31'
    Map.centerObject(geometry,6);
    var lai = dataset.select('Lai')
                     .filterDate(start_date,end_date)
                     .map(function(img){
                       return img.clip(geometry)
                     });
    print(lai)
    
    /*
    Visualization premeter data
    **/
    var laiVisualizationVis = {
      min: 0.0,
      max: 100.0,
      palette: ['e1e4b4', '999d60', '2ec409', '0a4b06'],
    };
    Map.addLayer(
        lai.mean(), laiVisualizationVis, 'LAI');
        
    // rebuid images for gap (no data)
    function bestNeighbour(image) {
      var collection = ee.ImageCollection.fromImages(image.get('neighbour'));
      image = image.addBands(collection.select(["Lai"]).mean().rename("mean"));
      return image;
    }
    
    function rebuildImages(imgCol) {
      var windowSize =  1;
      var timeField = "system:time_start";
      
      var join = ee.Join.saveAll({
        matchesKey: 'neighbour'
      });
      var selFilter = ee.Filter.maxDifference({
        difference: windowSize,
        leftField: timeField, 
        rightField: timeField
      });
      var selNeighbour = join.apply({
        primary: imgCol, 
        secondary: imgCol, 
        condition: selFilter
      });
      
      var meanCol = ee.ImageCollection(selNeighbour)
                      .map(bestNeighbour)
                      .select(["mean"], ["Lai"]);
      return meanCol;
    }
    var collection = rebuildImages(lai);
    
    // Step 2: Set up Savitzky-Golay smoothing 设置Savitzky-Golay平滑
    var window_size = 9
    var half_window = (window_size - 1)/2
    
    // Define the axes of variation in the collection array.在集合数组中定义变化轴。
    var imageAxis = 0;
    var bandAxis = 1;
    
    // add first and last
    var y0 = collection.first();
    var firstvals = collection.toList(2000).slice(1, half_window+1).reverse().map(
      function(e) { 
        var dif = ee.Number(0).subtract(ee.Date(ee.Image(e).get('system:time_start')).difference(ee.Date(y0.get('system:time_start')), 'day'))
        var dstamp = ee.Date(y0.get('system:time_start')).advance(dif,'day')
        return ee.Image(e).subtract(y0).abs().multiply(-1).add(y0).set("system:time_start",dstamp) 
      }
    );
    var yend = ee.Image(collection.toList(2000).slice(-1).get(0));
    // print(yend)
    var lastvals = collection.toList(2000).slice(-half_window-1,-1).reverse().map(
      function(e) { 
        var dif = ee.Number(0).subtract(ee.Date(ee.Image(e).get('system:time_start')).difference(ee.Date(yend.get('system:time_start')), 'day'))
        var dstamp = ee.Date(yend.get('system:time_start')).advance(dif,'day')
        return ee.Image(e).subtract(yend).abs().add(yend).set("system:time_start",dstamp)  
      }
    );
    collection = ee.ImageCollection(firstvals).merge(collection).merge(ee.ImageCollection(lastvals))
    
    
    var col_res = collection.map(function(img) {
      var dstamp = ee.Date(img.get('system:time_start'))
      var ddiff = dstamp.difference(ee.Date(start_date), 'day')
      img = img.select('Lai').set('date', dstamp)
      return img.addBands(ee.Image(1).toFloat().rename('constant')).
        addBands(ee.Image(ddiff).toFloat().rename('t')).
        addBands(ee.Image(ddiff).pow(ee.Image(2)).toFloat().rename('t2')).
        addBands(ee.Image(ddiff).pow(ee.Image(3)).toFloat().rename('t3'))
    })
    
    // Set polynomial order 设置多项式阶数
    var order = 3
    var coeffFlattener = [['constant', 'x', 'x2', 'x3']]
    var indepSelectors = ['constant', 't', 't2', 't3']
    
    // Change to order = 2 as follows:
    //var order = 2
    //var coeffFlattener = [['constant', 'x', 'x2']]
    //var indepSelectors = ['constant', 't', 't2']
    
    // Convert the collection to an array.
    var array = col_res.toArray();
    var arrayLength = array.arrayLength(imageAxis);
    
    // Update the mask to ensure that the number of images is greater than or
    // equal to the number of predictors (the linear model is solveable).
    array = array.updateMask(arrayLength.gt(4));
    
    // Solve 
    function getLocalFit(i) {
      // Get a slice corresponding to the window_size of the SG smoother 
      //获取与SG平滑器的window_size对应的切片
      var subarray = array.arraySlice(imageAxis, ee.Number(i).int(), ee.Number(i).add(window_size).int())
      var predictors = subarray.arraySlice(bandAxis, 1, 1 + order + 1)
      var response = subarray.arraySlice(bandAxis, 0, 1); // NDVI
      var coeff = predictors.matrixSolve(response)
      coeff = coeff.arrayProject([0]).arrayFlatten(coeffFlattener)
      return coeff  
    }
    
    // For the remainder, use col_res as a list of images 
    // 对于其余部分，请使用modis_res作为图像列表
    col_res = col_res.toList(col_res.size()) // 返回集合中的元素数
    var runLength = ee.List.sequence(0, col_res.size().subtract(window_size))
    
    // Run the SG solver over the series, and return the smoothed image version 
    // 在系列上运行SG解算器，并返回平滑的图像版本
    var sg_series = runLength.map(function(i) {
      var ref = ee.Image(col_res.get(ee.Number(i).add(half_window)))
      return getLocalFit(i).multiply(ref.select(indepSelectors)).reduce(ee.Reducer.sum()).copyProperties(ref)
    })
    
    // Part 3: Generate some output
    // 3A. Get an example original image and its SG-ed version 获取示例原始图像及其SG版本
    var Lai = (ee.Image(col_res.get(5+half_window)).select('Lai')) // half-window difference in index
    var Laisg = (ee.Image(sg_series.get(5))) 
    var mean = ee.Image(sg_series.get(0)).reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: geometry,
      scale:60,
      maxPixels: 1e13,
      tileScale:8
    });
    // print(mean)
    
    // show all NDVIsg images
    // for(var i = 0;i<30;i++){
    //   Map.addLayer(ee.Image(sg_series.get(0)),{},i+"sg")
    // }
    
    // chart
    // 3B. Now get a profile for a buffered point. 
    // Make sure to place the point in the centre of a homogenous area (reduces noise) 
    
    // click point
    Map.style().set("cursor", "crosshair");
    var pointPoi = null;
    Map.onClick(function(coords) {
      if (pointPoi !== null) {
        Map.remove(pointPoi);
        pointPoi = null;
      }
      var point = ee.Geometry.Point(coords.lon, coords.lat);
      pointPoi = Map.addLayer(point, {color:"red"}, "point");
      // var pt = /* color: #d63000 */geometry2.buffer(500);
       // Build a stack for all images in the collection
      function stack(i1, i2)
      {
        return ee.Image(i1).addBands(ee.Image(i2))
      }
      
      var s1orig = col_res.slice(1).iterate(stack, col_res.get(0))
    
      var s1sged = sg_series.slice(1).iterate(stack, sg_series.get(0))
      
      // Get samples from both series
      var y = ee.Image(s1orig).select(['Lai(..)*']).reduceRegion(ee.Reducer.mean(), point,500).values().slice(4,col_res.size().subtract(4)).reverse()
      var xlabels = ee.Image(s1orig).select(['t(..)*']).reduceRegion(ee.Reducer.first(), point,500).values().slice(4,col_res.size().subtract(4))
      var smoothy = ee.Image(s1sged).select(['sum(..)*']).reduceRegion(ee.Reducer.mean(), point,500).values().reverse()
      //y = y.map(function(f) { return ee.Number(f).log10().multiply(10)})
      xlabels = xlabels.map(function(f) { 
        
        return ee.Date(y0.get('system:time_start')).advance(ee.Number(f),'day').format('yyyy-MM-dd')
        
      }).sort()
      print(xlabels)
      //smoothy = smoothy.map(function(f) { return ee.Number(f).log10().multiply(10)})
      
      // Chart
      var yValues = ee.Array.cat([y, smoothy], 1);
      
      var smoothChart = ui.Chart.array.values(yValues, 0, xlabels).setSeriesNames(['Raw', 'Smoothed']).setOptions(
        {
          title: 'Savitsky-Golay smoothing (order = ' + order + ', window_size = ' + window_size + ')', 
          hAxis: {title: 'Time (date after ' + start_date + ')'}, vAxis: {title: 'Lai'},
          legend: null,
          series: { 
            0: { lineWidth: 1,pointSize:2},
            1: { lineWidth: 1, pointSize: 0, color: 'red' }}
        })
      print("smoothChart", smoothChart);
      //LAI chart
      var LAIchart = ui.Chart.image.series({
                    imageCollection: lai,
                    region:point,
                    reducer: ee.Reducer.mean(),
                    scale: 500
                  })
                  .setOptions({
                    title: "LAI列表", 
                    hAxis: {title: "date"},
                    vAxis: {title: "LAI value"},
                    lineWidth:1,
                    pointSize:2
                  });
      print("LAIchart", LAIchart);
    });
### 在用上述代码的时候请引用以下文献:
Fan, Z.; Li, S.; Fang, H. Explicitly Identifying the Desertification Change in CMREC Area Based on Multisource Remote Data. Remote Sens. 2020, 12, 3170.