基于GEE平台提取水体

本文主要介绍如何利用GEE平台与哨兵2号影像提取水体。水体的提取主要是基于NDWI指数进行，当然，通过改变波段运算，也可以根据MNDWI进行提取。主要代码如下：
下面展示一些 内联代码片。

//首先对数据进行去云处理   
var s2_rmcloud = function(image) {
 var quality = image.select("QA60").unmask();
 return image.updateMask(quality.eq(0));
};
//NDWI
var s2_ndwi = function(image) {
 return image.addBands(image.normalizedDifference(["B3", "B8"]).rename("NDWI"));
};
//筛选数据 

var s2_nocloud = s2_col.map(s2_rmcloud)
                       .filterBounds(roi)
                       .filterDate("2020-11-20", "2020-11-25");
var ndwi = s2_nocloud.map(s2_ndwi).select("NDWI").reduce(ee.Reducer.mean()).clip(roi);
print(ndwi)
var visParam = {
 min: -0.5,
 max: 0.5,
 palette: ['00FFFF', '0000FF']
};
Map.addLayer(ndwi, visParam, "ndwi");
Map.centerObject(roi,12);
//otsu算法 
var otsu = function(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);
  
  // Compute between sum of squares, where each mean partitions the data.
  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)));
  });
  
  print(ui.Chart.array.values(ee.Array(bss), 0, means));
  
  // Return the mean value corresponding to the maximum BSS.
  return means.sort(bss).get([-1]);
};
//归一化，去除异常值  
function normalization(image,region,scale){
var mean_std = image.reduceRegion({
  reducer: ee.Reducer.mean()
            .combine(ee.Reducer.stdDev(),null, true),
  geometry: region,
  scale: scale,
  maxPixels: 1e13,
  // tileScale: 16
}); 
var unitScale = ee.ImageCollection.fromImages(
  image.bandNames().map(function(name){
    name = ee.String(name);
    var band = image.select(name);
    var mean=ee.Number(mean_std.get(name.cat('_mean')));
    var std=ee.Number(mean_std.get(name.cat('_stdDev')));
    var max=mean.add(std.multiply(3))
    var min=mean.subtract(std.multiply(3))
    var band1=ee.Image(min).multiply(band.lt(min)).add(ee.Image(max).multiply(band.gt(max)))
                        .add(band.multiply(ee.Image(1).subtract(band.lt(min)).subtract(band.gt(max))))
    var result_band=band1.subtract(min).divide(max.subtract(min));
    return result_band;
})).toBands().rename(image.bandNames());
  return unitScale;
}
var ndwi2=ndwi.select('NDWI_mean')
var ndwi3=normalization(ndwi2,roi,10)
var ndwi1=ndwi3.multiply(10000)
print(ndwi1)
var histogram = ndwi1.select('NDWI_mean').reduceRegion({
  reducer: ee.Reducer.histogram(255, 2)
      .combine('mean', null, true)
      .combine('variance', null, true), 
  geometry: roi, 
  scale: 10,
  bestEffort: true
});
print(Chart.image.histogram(ndwi1.select('NDWI_mean'), roi, 12));
print(histogram);
var threshold = otsu(histogram.get('NDWI_mean_histogram'));
print('threshold', threshold);
var water_mask = ndwi1.select('NDWI_mean').gt(threshold);

Map.addLayer(water_mask.mask(water_mask), {palette: 'red'}, 'water');
Export.image.toDrive({
  image: water_mask.mask(water_mask),
  description: "202009",
  folder:"UltimateS2_R2",
  fileNamePrefix: "202009",
  scale: 10,
  region:roi,
  maxPixels: 1e13});
  Map.addLayer(water_mask.mask(water_mask), {palette: 'red'}, 'class A');
var countDictionary = water_mask.mask(water_mask).reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: roi,
  scale: 10,
  maxPixels: 1e13
});
//把统计结果打印出来
print(countDictionary);

虽然看着很长，但是基本思想时利用阈值法提取水体，借鉴了Otsu算法确定阈值，另外加了剔除异常值的归一化算法。
基于哨兵2号1-C的数据源，作者在写代码的时候发现哨兵2号2A和哨兵2号1C的数据在去云后数据数量上不相同了，暂时不知道是什么原因，欢迎各位小伙伴来交流哇。