GEE第一段程序:沙洲面积提取

GEE第一段程序:沙洲面积提取

打开GEE,先选定一块区域(使用ROI在地图上选定)可以看到此时代码框中出现了Imports (1 entry) ,我们可以对这个变量var geometry改名,例如这里我们将之改为roi,执行下边的代码:

Map.centerObject(roi, 13);
var Area_roi=roi.area().divide(1000000);
print(Area_roi);
var bounds= ee.Image().toByte()
 .paint({
     featureCollection:
    ee.FeatureCollection([ee.Feature(roi)]),
     color: null,
     width: 1
 });
Map.addLayer(bounds, {palette: "red"}, "bounds");
// 选择遥感影像的时间段(startDate 包含输入时间,endDate 不包含输入时
//间。下图中的遥感影像选择时段实际上是 2018 年 12 月 28 日-2019 年 5 月 26
//日)
var startDate = ee.Date("2018-12-28");
var endDate =ee.Date("2019-5-27");
// 引入所有遥感影像源(这些遥感影像已经经过了预处理,可以直接利用)
var l5col=ee.ImageCollection("LANDSAT/LT05/C01/T1_SR");
var l7col= ee.ImageCollection("LANDSAT/LE07/C01/T1_SR");
var l4col=ee.ImageCollection("LANDSAT/LT04/C01/T1_SR");
var l8col =ee.ImageCollection("LANDSAT/LC08/C01/T1_SR");
//选择某一遥感影像源,进行影像的区域,云量筛选,并进行去云处理,计算
// MNDWI(即影像二值化).
var Col=ee.ImageCollection(l8col)
    .filterBounds(roi)
    .filterDate(startDate, endDate)
    .filter(ee.Filter.lte("CLOUD_COVER", 20))
    .map(
        function(image) {
         var cloudShadowBitMask = 1 << 3;
         var cloudsBitMask = 1 << 5;
         var snowBitMask = 1 << 4;
         var qa = image.select('pixel_qa');
         var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
            .and(qa.bitwiseAnd(cloudsBitMask).eq(0))
            .and(qa.bitwiseAnd(snowBitMask).eq(0));
         return image.updateMask(mask);
 }
 )
 .map(function(image) {
 var time_start = image.get("system:time_start");
 image = image.multiply(0.0001);
 image = image.set("system:time_start", time_start);
 return image;
 })
 .map(function(image) {
 return image.addBands(image.normalizedDifference(["B3", "B6"])
 .rename("MNDWI"));
 })
 .select("MNDWI");
print("Col", Col);
//大津算法(otsu)选择最优阈值方法,并进行影像的分割,水体和陆地分
// 离,提取陆地面积。
function otsu(histogram) {
     var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
     var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
     var size = means.length().get([0]);
     var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
     var sum = means.multiply(counts).reduce(ee.Reducer.sum(),
    [0]).get([0]);
     var mean = sum.divide(total);
     var indices = ee.List.sequence(1, size);
     var bss = indices.map(function(i) {
     var aCounts = counts.slice(0, 0, i);
     var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
     var aMeans = means.slice(0, 0, i);
     var aMean = aMeans.multiply(aCounts)
     .reduce(ee.Reducer.sum(), [0]).get([0])
     .divide(aCount);
     var bCount = total.subtract(aCount);
     var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
     return aCount.multiply(aMean.subtract(mean).pow(2)).add(
     bCount.multiply(bMean.subtract(mean).pow(2)));
 });
 return means.sort(bss).get([-1]);
}
var sCol1 = Col.map(function(image) {
         var histogram = image.reduceRegion({
         reducer: ee.Reducer.histogram(),
         geometry: roi,
         scale: 2,
         maxPixels: 1e13,
         tileScale: 10
         });
         var threshold = otsu(histogram.get("MNDWI"));
         var mask = image.lte(threshold);
         var water = mask.updateMask(mask).rename("water");
         var newWater = water.addBands(image);
        //转矢量
         var vectors = newWater.reduceToVectors({
         geometry: roi,
         scale: 2,
         geometryType: 'polygon',
         reducer: ee.Reducer.mean(),
         maxPixels: 1e13,
         tileScale: 10
     });
     vectors = vectors.filterBounds(roi);
     water = water.clip(vectors);
     //计算陆地面积
     var areaImage = water.multiply(ee.Image.pixelArea());
     //统计指定区域中所有陆地像素和,也就是沙体的面积
     var dict = areaImage.reduceRegion({
         reducer: ee.Reducer.sum(),
         geometry: vectors.geometry(),
         scale: 2,
         maxPixels: 1e13,
         tileScale: 10
     });
     water = water.set("area",
    ee.Number(dict.get("water")).divide(1000000));
     water = water.set("threshold", threshold);
     return water.toByte();
});
 print("sCol1", sCol1);
 Map.addLayer(sCol1.first(), {palette: "red"}, "BAR");
// 整合所有影像的时间和计算出的沙洲面积并输出
Export.table.toDrive({collection: sCol1,
     description: "guanzhou",
     fileNamePrefix: "guanzhou",
     selectors:["system:index","area"]
 });
 

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注