遥感影像SG滤波（基于GEE）

var s2 = ee.ImageCollection("COPERNICUS/S2"); 
var geometry = ee.Geometry.Polygon([[
  [82.60642647743225, 27.16350437805251],
  [82.60984897613525, 27.1618529901377],
  [82.61088967323303, 27.163695288375266],
  [82.60757446289062, 27.16517483230927]
]]);
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)

var startDate = ee.Date('2017-01-01')
var endDate = startDate.advance(1, 'year')

var filtered = s2
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
  .filter(ee.Filter.bounds(geometry))

// Write a function for Cloud masking
function maskS2clouds(image) {
  var qa = image.select('QA60')
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0))
  return image.updateMask(mask)
      .select("B.*")
      .copyProperties(image, ["system:time_start"])
}

var filtered = filtered.map(maskS2clouds)

function addNDVI(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).toFloat().rename('ndvi');
  return image.addBands(ndvi);
}

var withNdvi = filtered.map(addNDVI);

var ndviCol = withNdvi.select('ndvi')
print('Original Collection', ndviCol)

var n = 5;
var totalDays = endDate.difference(startDate, 'day');
var daysToInterpolate = ee.List.sequence(1, totalDays, n)

var initImages = daysToInterpolate.map(function(day) {
  var image = ee.Image().rename('ndvi').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 = ndviCol.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)
})

var days = 60
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('ndvi')
print('Interpolated Collection', interpolatedCol)

var oeel=require('users/OEEL/lib:loadAll');
// https://www.open-geocomputing.org/OpenEarthEngineLibrary/#.ImageCollection.SavatskyGolayFilter

// 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_ndvi波段并重命名
var smoothed = smoothed.select(['d_0_ndvi'], ['smoothed'])

var title = 'Savitsky-Golay smoothing' +
  '(order = '+ order + ', window_size = ' + days + ')'

// 把原始的和经过SG滤波的数据NDVI曲线画出来
var chart = ui.Chart.image.series({
  imageCollection: ndviCol.merge(smoothed),
  region: geometry,
  reducer: ee.Reducer.mean(),
  scale: 20
}).setOptions({
      lineWidth: 1,
      title: title,
      interpolateNulls: true,
      vAxis: {title: 'NDVI'},
      hAxis: {title: '', format: 'YYYY-MMM'},
      series: {
        0: {color: 'blue', lineWidth: 1, 
          lineDashStyle: [1, 1], pointSize: 1,
          }, // Original NDVI
        1: {color: 'red', lineWidth: 2 }, // Smoothed NDVI
      },

    })
print(chart)


var ndviVis = {min:0, max: 0.7,  palette: [
    '#fdae61','#fee08b','#d9ef8b',
    '#a6d96a','#66bd63','#1a9850','#006837']
  }

var visualizeImage = function(image) {
  return image.visualize(ndviVis).clip(geometry)
}

var visCollectionOriginal = ndviCol
  .map(visualizeImage)
var visCollectionInterpolated = smoothed
  .map(visualizeImage)
  
var videoArgs = {
  dimensions: 400,
  region: geometry,
  framesPerSecond: 2,
  crs: 'EPSG:3857',
};

print('Original NDVI Images', visCollectionOriginal.getVideoThumbURL(videoArgs))

var videoArgs = {
  dimensions: 400,
  region: geometry,
  framesPerSecond: 5,
  crs: 'EPSG:3857',
};
print('SG Filtered NDVI Images', visCollectionInterpolated.getVideoThumbURL(videoArgs))