前两期我们介绍了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,
}));
结果如图:
本想拿五年的月平均数据给大家展示,奈何两年的都会计算超时,看来哨兵数据还是太多了。只好拿Landsat的产品再看看效果了(五年确实也没算出来,2018-2020还是可以的),看效果感觉还不错。
参考博客:
https://mygeoblog.com/2021/05/30/whittaker-smoothing-algorithm-for-gee/
原始代码平滑算法写的较多,大家可以认真研究研究,不想研究的同学可以直接把数据集部分修改好就行。公众号回复122901即可取走。
更多精彩内容请关注: