GEE计算时间序列植被指数-以哨兵2数据计算MTCI指数为例
// 首先是去云处理
function rmCloudByQA(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);}
// 计算MTCI指数(选择哨兵2数据的第4,5,6波段)
function s2_mtci(image) {
var mtci = image.expression('(RE2 - RE1)/ (RE1 - Red)', {
'RE2': image.select('B6'),
'RE1': image.select('B5'),
'Red': image.select('B4'),
}).float().rename('MTCI'); // 这里强调一点需要利用rename()将结果名字改成‘MTCI’,原因后面再说。
return image.addBands(mtci);
} // 这里强调,函数返回的如果直接是image的话,后面在进行时序分析就会报错,说没有波段匹配,因为缺少了相应的时间信息。所以这里是将计算的结果作为一个波段添加进原影像中。
// 导入研究区
var roi = ee.FeatureCollection("users/area");
// 导入哨兵2数据,根据时间、研究区和云量对影像进行筛选
var s2_col = ee.ImageCollection("COPERNICUS/S2")
.filterBounds(roi)
.filterDate('2022-05-14','2022-09-21')
.filter(ee.Filter.lte('CLOUD_COVERAGE_ASSESSMENT',20))
.map(rmCloudByQA)
.map(function(image){
return image.clip(roi)
});
print (s2_col)
// 计算研究区的MTCI指数
var mtci = s2_col.map(s2_mtci).select("MTCI"); // 上面只有rename()函数使用了以后,才能在这里根据名字进行select,不然的话,mtci指数的波段名字这里就是显示的'B6_1',不知道为什么是这样的。
print ('mtci is:', mtci)
Map.addLayer(mtci)
// 站点时间序列显示(以日期和日序数两种)
var chart1 = ui.Chart.image.series({
imageCollection:mtci,
region:ee.Geometry.Point([107.29424,40.73286]),
reducer:ee.Reducer.mean(), //计算均值
scale:10
}).setOptions({title:'MTCI IMAGE SERIES'});
print(chart1);
var chart2 = ui.Chart.image.doySeries({
imageCollection:mtci,
region:ee.Geometry.Point([107.29424,40.73286]),
regionReducer:ee.Reducer.mean(),
scale:10
}).setOptions({title:"ROI MTCI EACH DAY SERIES"})
print(chart2)
**结果:**可以看到选择以后只有mtci一个波段数据;时间序列数据也可以得到,通过右上角可以打开新的界面,导出图以及相应的数据表。
多点像元时间序列值提取
上述代码只是提供了单点的时间序列的表示,下面这个代码是在上述代码的基础上,提取多点像元时间序列,并导出。代码参考的博文GEE:提取多个点的时间序列数据
// 提取多点像元值
var pointValues = mtci.toBands().sampleRegions({ // 这里必须使用toBands才能进行计算,没搞懂这个是啥意思。
collection:syz,
/* properties:ee.List(['Number']),*/ //这里如果加上properties的话应该是可以根据矢量点的属性信息导出数据的,如果不加的话,属性信息应该都会导出,这里没有测试过。
scale:10
});
print(pointValues)
//导出
Export.table.toDrive({
collection: pointValues,
description:"MTCI",
folder: "MTCI",
fileFormat: "CSV"
});
导出的数据的编号和print(pointValues)中的feature是一致的,可以对照。我的矢量点文件中也有Number编号,所以也可以对照。如图:
申明: 自己也是查询了很多代码,加上自己的理解得到的结果,有问题大家可以指出来,互相交流学习,自己关注的是以下几个公众号来学习,推荐一下。