GEEer成长日记七：Whittaker平滑应用于Sentinel-2时间序列分析

爱探索的GEEer

已于 2022-02-15 01:16:24 修改

阅读量1.9k

点赞数 4

分类专栏： Google earth engine 文章标签： google earth javascript 云计算时序模型

于 2022-02-08 00:26:35 首次发布

本文链接：https://blog.csdn.net/weixin_51691876/article/details/122816864

版权

Google earth engine 专栏收录该内容

23 篇文章 242 订阅

订阅专栏

　　前两期我们介绍了Sentinel-2 Level-2A数据在逐月和逐日时间序列方面的研究。因为数据集在过滤的时候，我们先筛选的云量，再进行去云处理，可能导致有些时间段内影像缺失严重，致使NDVI的数值偏差较大。可按照常理来说，NDVI的变化应该是平滑的，没有特殊情况的发生，很难大起大落。

　　于是有同学小窗问小编，时间序列可否加入滤波后再显示。大概意思是要平滑的意思吧。今天我们就介绍一个平滑器在时间序列的应用，希望可以帮到大家。

Whittaker Smoother

具体原理部分我就不赘述了，可以看官方简介：

https://eigenvector.com/wp-content/uploads/2020/01/WhittakerSmoother.pdf

参考文献也给大家附上：

https://www.mdpi.com/2072-4292/12/18/2888/htm

话不多说，我们继续搞代码：

//还是老样子哈，以广东省2020年为目标
var geometry = ee.FeatureCollection('users/ZhengkunWang/guangdongsheng')
Map.centerObject(geometry,6)
//此处省略几百行代码。。。。。。
//中间平滑的算法内容较多，就不复制粘贴了，大家在文末回复122901即可获取code
//使用QA波段去云
function maskS2clouds(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()));
}
var S2 = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2019-01-01', '2020-12-31')
.filterBounds(geometry)
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
.map(maskS2clouds);
var addNDVI = function(image) {
return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'));
};
// Add NDVI band to image collection
var S2 = S2.map(addNDVI);
var NDVI = S2.select(['NDVI']);

var years = ee.List.sequence(2020, 2020);
var months = ee.List.sequence(1, 12);
var S2_monthlymeanNDVI =  ee.ImageCollection.fromImages(
  years.map(function (y) {
    return months.map(function(m) {
    return NDVI.filter(ee.Filter.calendarRange(y,y, 'year')).filter(ee.Filter.calendarRange(m, m, 'month')).mean().set('year', y).set('month', m).set('system:time_start', ee.Date.fromYMD(y, m, 1));
    });
  }).flatten()
);
// getting rid of masked pixels
ndvi = S2_monthlymeanNDVI.map(function(img){return img.unmask(S2_monthlymeanNDVI.mean())});

var ndvi =  whittakerSmoothing(ndvi)[0];

// add chart
print(ui.Chart.image.series(
  ndvi.select(['NDVI', 'NDVI_fitted']), geometry, ee.Reducer.mean(), 1000)
    .setSeriesNames(['NDVI', 'NDVI_smoothed'])
    .setOptions({
      title: 'Sentinel-2 NDVI time series smoothed',
      lineWidth: 1,
      pointSize: 3,
}));

　　结果如图：