前言:
由于我要使用S2在2016-2020年时序数据进行作物分类,但GEE仅提供2018年后的S2 L2A数据('COPERNICUS/S2_SR'),因此需要进行大气校正,同时由于需要的数据量大,因此需要批量处理。
大气校正的方法论文:(英文)一种传感器不变的大气校正方法:适用于哨兵-2/MSI和陆地卫星8/OLI (researchgate.net)
我的代码:(主要包括去云、大气校正、矢量裁剪、均值融合、显示和数据下载)
//加载 shpfile
var shpfile1 = ee.FeatureCollection('projects/ee-wenqikou/assets/qixia_0_09').filter(ee.Filter.eq('grid_label','6')).geometry()
//聚焦shp范围
Map.centerObject(shpfile1)
//显示shp
Map.addLayer(shpfile1)
//引入大气纠正方法
var siac = require('users/marcyinfeng/utils:SIAC');
//去云函数
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);}
//数据函数,参数(开始日期,截止日期,云量筛选,融合形式,图像显示选择波段(这里仅支持真彩色))
// 其中,融合方式compositiontype: 0 min 1 mean 2 median 3 max
var compsition = function(startdate,enddate,cloudiness,compositiontype,bands){
var image = ee.ImageCollection('COPERNICUS/S2').filterDate(startdate,enddate).filterBounds(shpfile1).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', cloudiness))
//去云
image = image.map(rmCloudByQA)
//批量大气校正
image = image.map(siac.get_sur)
//融合(这样可以免本地拼接等
if(compositiontype === 0){image = image.min().multiply(10000)}
else if(compositiontype == 1){image = image.mean().multiply(10000)} //均值融合
else if(compositiontype == 2){image = image.median().multiply(10000)}
else image = image.max().multiply(10000)
//裁剪
image = image.clip(shpfile1);
//数据显示
if(bands == 432) {Map.addLayer(image,{min:0, max:3000, bands:["B4", "B3", "B2"]},'real'+startdate);}
//数据下载
Export.image.toDrive({
image: image,
description: startdate,
folder: 'qixia/grid/image/06',
crs:"EPSG:4326",
fileNamePrefix: startdate,
region: shpfile1,
scale: 10,
maxPixels: 1e13
});
return image
}
var ys = compsition('2021-02-15','2021-02-16',100,1,432)
关于Landsat7和Landsat8:
Landsat7:仅需要修改“siac.get_sur”为“siac.get_l7_sur”
Landsat8:仅需要修改“siac.get_sur”为“siac.get_l8_sur”
例子:
var ground_sensor = ee.Geometry.Point(-52.905090, -28.228550);
var polygon = ground_sensor.buffer(1000).bounds();//
var start_date = '2017-11-15';
var end_date = '2017-11-30';
var criteria = ee.Filter.and(
ee.Filter.bounds(ground_sensor), ee.Filter.date(start_date, end_date));
var cloud_perc = 60;//Max cloud percentile per scene.
// -- Collections of Landsat 7, 8 and Sentinel 2 TOA data
var L8_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filter(criteria)
.filter(ee.Filter.lt('CLOUD_COVER', cloud_perc));
var L7_col = ee.ImageCollection('LANDSAT/LE07/C01/T1_TOA')
.filter(criteria)
.filter(ee.Filter.lt('CLOUD_COVER', cloud_perc));
var S2_col = ee.ImageCollection("COPERNICUS/S2")
.filter(criteria)
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', cloud_perc));
// - Import the SIAC atmospheric correction module
var siac = require('users/marcyinfeng/utils:SIAC');
// - Apply SIAC and retrieve bottom of atmosphere (BOA) reflectance
var L7_boa = siac.get_l7_sur(L7_col.first());
var L8_boa = siac.get_l8_sur(L8_col.first());
var S2_boa = siac.get_sur(S2_col.first());
// - Check and visualization
var Color_comp_01 = {bands:"B4,B3,B2", min: 0.0, max: 0.2, gamma: 1};
var Color_comp = {bands:"B4,B3,B2", min:200, max:2000, gamma: 1};
Map.addLayer(S2_col.first(), Color_comp, 'TOA');
Map.addLayer(S2_boa, Color_comp_01, 'BOA');
Map.centerObject(S2_col.first())