代码语言:javascript复制
//输入目标区域
var geometry = table;
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var filtered = s2
.filter(ee.Filter.date('2019-07-01', '2019-08-31'))//修改数据时间范围
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.bounds(geometry))
// Write a function for Cloud masking
function maskCloudAndShadowsSR(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()))
.select("B.*")
.copyProperties(image, ["system:time_start"]);
}
var filtered = filtered.map(maskCloudAndShadowsSR)
//在有效像素下添加时间信息,便于后续利用
var filtered = filtered.map(function(image) {
var timeImage = image.metadata('system:time_start').rename('timestamp')
var timeImageMasked = timeImage.updateMask(image.mask().select(0))
return image.addBands(timeImageMasked)
})
var days = 30 //前后时段长
var millis = ee.Number(days).multiply(1000*60*60*24)
var maxDiffFilter = ee.Filter.maxDifference({
difference: millis,
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var lessEqFilter = ee.Filter.lessThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var greaterEqFilter = ee.Filter.greaterThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var filter1 = ee.Filter.and(maxDiffFilter, lessEqFilter)
var join1 = ee.Join.saveAll({
matchesKey: 'after',
ordering: 'system:time_start',
ascending: false})
var join1Result = join1.apply({
primary: filtered,
secondary: filtered,
condition: filter1
})
var filter2 = ee.Filter.and(maxDiffFilter, greaterEqFilter)
var join2 = ee.Join.saveAll({
matchesKey: 'before',
ordering: 'system:time_start',
ascending: true})
var join2Result = join2.apply({
primary: join1Result,
secondary: join1Result,
condition: filter2
})
//插值
var interpolateImages = function(image) {
var image = ee.Image(image)
// We get the list of before and after images from the image property
// Mosaic the images so we a before and after image with the closest unmasked pixel
var beforeImages = ee.List(image.get('before'))
var beforeMosaic = ee.ImageCollection.fromImages(beforeImages).mosaic()
var afterImages = ee.List(image.get('after'))
var afterMosaic = ee.ImageCollection.fromImages(afterImages).mosaic()
// Get image with before and after times
var t1 = beforeMosaic.select('timestamp').rename('t1')
var t2 = afterMosaic.select('timestamp').rename('t2')
var t = image.metadata('system:time_start').rename('t')
var timeImage = ee.Image.cat([t1, t2, t])
var timeRatio = timeImage.expression('(t - t1) / (t2 - t1)', {
't': timeImage.select('t'),
't1': timeImage.select('t1'),
't2': timeImage.select('t2'),
})
// Compute an image with the interpolated image y
var interpolated = beforeMosaic
.add((afterMosaic.subtract(beforeMosaic).multiply(timeRatio)))
// Replace the masked pixels in the current image with the average value
var result = image.unmask(interpolated)
return result.copyProperties(image, ['system:time_start'])
}
// function addNDVI(image){
// var ndvi = image.normalizedDifference(['B8','B4']).rename('ndvi');
// var doy = image.date().format("yyyyMMdd")
// return image.addBands(ndvi).select('ndvi').clip(geometry)
// .set({'Date':doy});
// }
function addNDBI(image){
// 方法 3 自己手动加减乘除
var ndbi = image.select('B11').subtract(image.select('B8')).divide(image.select('B11').add(image.select('B8')))
.float().rename('NDBI');
var doy = image.date().format("yyyyMMdd")
return image.addBands(ndbi).select('NDBI').clip(geometry)
.set({'Date':doy});
}
var interpolatedCol = ee.ImageCollection(
join2Result.map(interpolateImages)).map(addNDBI)
print(interpolatedCol)
//数据输出函数
function exportImage(image,roi,fileName){
Export.image.toDrive({
image: image,
description: fileName '_NDBI',
folder: 'NDBI',
scale: 10, //分辨率
region: roi,
maxPixels:1e13,
crs:'EPSG:4326', //坐标系
fileFormat: 'GeoTIFF',
});
}
//生成列表,迭代下载
var indexList = interpolatedCol.reduceColumns(ee.Reducer.toList(),["Date"]).get("list");
print("indexList",indexList);
indexList.evaluate(function(indexs){
var indexs1 = ee.List(indexs).distinct();
var length = indexs1.length().getInfo();
print(indexs1)
for (var i=0;i<length;i ){
var img =interpolatedCol.filter(ee.Filter.eq("Date",indexs1.get(i))).mean()
//Map.addLayer(img,indexs1,get(i).getInfo(),false);
exportImage(img,geometry,indexs1.get(i).getInfo())
}
})