开篇之作,笔者标码一枚,GEE开发工作四载有余,涉及项目百来余个。独乐乐不如众乐乐,遂分享之。
/*******************************************
* 愿诸君 :节之鼓而行之,以遂八风。
******************************************/
言归正传
对于GEE批量处理而言,JavaSript中for虽然是循环标配,但其API提及:云计算实现For Loops时,非不得已不能用For,因为使用for调用的是本地资源,而官方标配为map,那么存在以下几个矛盾:
- 使用map时,不能进行print或者Map.addLayer操作,对复杂逻辑而言,容易抛出异常
- 使用for拉低运行内存,稍有不慎,轻则弹窗警告之,重则遂现Limit
权衡利弊,还是使用map,但是需要稍加改动方可,话不多说,下面是一段可自定义时间区间的code。分为数据集定义、通用指数方程定义、批量下载实现和批量run task。
一、数据集定义
集结Landsat 4、5、7和8系列数据集,包括预处理 。此码无需改动
/****************************************************************************************************************************
* function of data collection
* **************************************************************************************************************************/
{{{ //preprocess dataset
var createLandsatCollection=function createLandsatCollection(params) {
var defaultParams = {
region: Map.getBounds(true),
start: '1982-01-01',
end: formatDate(new Date()),
mapImage: function (image) { return image }
}
params = mergeObjects([defaultParams, params])
var filter = ee.Filter.and(
ee.Filter.bounds(params.region),
ee.Filter.date(params.start, params.end )
)
var l4 = ee.ImageCollection('LANDSAT/LT04/C01/T1_SR')
.merge(ee.ImageCollection('LANDSAT/LT04/C01/T2_SR'))
.filter(filter)
.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
)
.map(cloudMaskL457)
var l5 = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
.merge(ee.ImageCollection('LANDSAT/LT05/C01/T2_SR'))
.filter(filter)
.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
)
.map(cloudMaskL457)
var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
.merge(ee.ImageCollection('LANDSAT/LE07/C01/T2_SR'))
.filter(filter)
.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
)
.map(cloudMaskL457)
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.merge(ee.ImageCollection('LANDSAT/LC08/C01/T2_SR'))
.filter(filter)
.select(
['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
)
.map(cloudMaskL8)
return l4.merge(l5).merge(l7).merge(l8).sort('system:time_start')
//Landsat4 5 7
function cloudMaskL457(image) {
var qa = image.select('pixel_qa');
// If the cloud bit (5) is set and the cloud confidence (7) is high
// or the cloud shadow bit is set (3), then it's a bad pixel.
var cloud = qa.bitwiseAnd(1 << 5)
.and(qa.bitwiseAnd(1 << 7))
.or(qa.bitwiseAnd(1 << 3));
// Remove edge pixels that don't occur in all bands
var mask2 = image.mask().reduce(ee.Reducer.min());
return image.updateMask(cloud.not()).updateMask(mask2);
};
//Landsat8
function cloudMaskL8(image) {
var cloudShadowBitMask = (1 << 3);
var cloudsBitMask = (1 << 5);
var qa = image.select("pixel_qa");
var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return image.updateMask(mask);
}
function excludeBand(bandName, image) {
var bandNames = image.bandNames()
var bandIndexes = ee.List.sequence(0, bandNames.size().subtract(1))
.filter(
ee.Filter.neq('item', bandNames.indexOf(bandName))
)
return image.select(bandIndexes)
}
function formatDate(date) {
var d = new Date(date),
month = '' + (d.getMonth() + 1),
day = '' + d.getDate(),
year = d.getFullYear()
if (month.length < 2)
month = '0' + month
if (day.length < 2)
day = '0' + day
return [year, month, day].join('-')
}
function mergeObjects(objects) {
return objects.reduce(function (acc, o) {
for (var a in o) { acc[a] = o[a] }
return acc
}, {})
}
}
}}}
二、通用指数方程定义
假设一个场景,如果在上述码块中,您没有使用.select[A,B]对 4,5,7,8的波段统一命名,那么首先任务一:计算NDVI,您是不是需要写2-3个fucntion对其NDVI的计算进行定义?如果还需要计算NDWI、NDBI、LSWI等等一系列的指数呢?实际上可以定义为带参的function,码块如下,使用方法见后附。此码无需改动
/****************************************************************************************************************************
* function of calculate index
* **************************************************************************************************************************/
{{{
//计算归一化指数 通用方程
var Norm_ID=function Norm_ID(B1,B2,name) {
var norm_id= function(img){
var temp_id= img.normalizedDifference([B1,B2]);
return img.addBands(temp_id.rename(name));
}
return norm_id;
}
//计算归比例指数 通用方程
var Two_Divide=function Two_Divide(B1,B2,name) {
var two_divide= function(img){
var temp_td=img.expression(
'b1/b2',
{b1:img.select(B1),
b2:img.select(B2)
});
return img.addBands(temp_td.rename(name));
}
return two_divide;
}
}}}
三、批量下载
上述function固定好以后,那么接下来即是数据整合了,在这里笔者使用了getInfo(),虽然API对getInfo()投了反对票,但是权衡引言中的矛盾一二,这里选择使用getInfo()是最佳方案。
/****************************************************************************************************************************
*
* defination Time range
*
* **************************************************************************************************************************/
var date_start='2000-01-01'
var date_end='2020-12-31'
var date_interval=12 //自定义时间片段,如 3 则为季度,6则为半年,12为1年
var roi=table //研究区自定义
//==========================================> step1 find date interval
var date_start_format=ee.Date(date_start)
var all_days=ee.Date(date_end).difference(date_start_format,'month')
var data_range=all_days.divide(date_interval).ceil();
var data_range_list=ee.List.sequence(1,data_range).getInfo()
// print(data_range_list)
//==========================================> step2 map the date
var ndvi_col=data_range_list.map(function(date_item){//******************for each data_range_list
var start_data=date_start_format.advance((date_item-1)*date_interval,'month')
var end_data=start_data.advance(date_interval,'month')
// print(start_data)
// print(end_data)
//========================================= landsat
var collection_landsat = createLandsatCollection({
region: roi,
start:start_data,
end:end_data,
})
// calculate index 这时使用定义的万能函数。
var ndvis_landsat=collection_landsat.map(Norm_ID('nir','red','NDVI'))
.map(Norm_ID("nir","swir1","LSWI"))
.map(Norm_ID("green","nir","NDWI"))
.map(Norm_ID("swir1","nir","NDBI"))
.map(Norm_ID("green","swir1","mNDWI"))
Export.image.toDrive({
image:ndvis_landsat.median().clip(roi),
folder:'NDVI_landsat',
fileNamePrefix :'landsat_'+date_start+'_'+date_item,
description:'landsat_'+date_start+'_'+date_item,
region:roi.geometry().bounds(),
scale:120,
crs:"EPSG:4326",
maxPixels:1e13
})
})
四、批量Run tasks
如果按季度下载,需要许run很多次。虽然之前有过处理code,但好像无效了。偶然看见杨志强老师的code,感谢杨老师无私分享
runTasks = function() {
const evt = new MouseEvent('click', {bubbles: true, cancelable: true, ctrlKey: true})
$$('.run-button' ,$$('ee-task-pane')[0].shadowRoot).forEach(function(e) {
e.dispatchEvent(evt)
})
}
runTasks()
尾记
如有不妥之处,请各位看官批评指正。