GEE随记（六）：Sentinel影像反演浊度

1. 选定探究区

这里随便框选了一个矩形

var styling ={color:'red',fillColor:'00000000',width:2};// display hollow roi
Map.addLayer(ee.FeatureCollection(geometry).style(styling), {}, "outline");
Map.centerObject(geometry,13)

2. 选取影像及去云处理&可视化

这里根据区域覆盖的影像选择日期为1月1日到5日

//去云处理 
function maskS2clouds(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);
}

//筛选数据 
var dataset = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
                  .filterDate('2024-01-01', '2024-01-05')
                  .filterBounds(geometry)
                  // Pre-filter to get less cloudy granules.
                  //.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
                  .map(maskS2clouds)
print(dataset)
var img = dataset.median().clip(geometry)
//可视化 
var visualization = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};
Map.addLayer(img, visualization, 'RGB');

 

3. 计算指数构建水体掩膜

（1）这里选择NIR+OSTU计算水体划分阈值，计算浊度并且添加两个波段

浊度计算公式根据自己反演结果，每个地方可能不一致，这里不展示公式

（2）OSTU算法创建水体掩膜提取水体

// otsu
function otsu(histogram) {
  var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));//frequency
  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]);
}
var band = 'B8'
var band_histo = 'B8_histogram'
var histogram = img2.select(band).reduceRegion({
  reducer: ee.Reducer.histogram()
      .combine('mean', null, true)
      .combine('variance', null, true),
  geometry: geometry, 
  scale: 10,
  bestEffort: true
});
//print(histogram);
print(Chart.image.histogram(img2.select(band),geometry,10));
var threshold = otsu(histogram.get(band_histo));
//get water layer
var classA = img2.select(band).lte(threshold);
var watermask1 = classA.updateMask(classA)
Map.addLayer(watermask1, {palette: 'blue'}, 'NIR_ostu');

注意：NIR+OSTU的方法提的水体是很充分的，往往把山体阴影和部分不透水面也提取了，可能不太准确，这里只是粗略示范 

4. 提取水体后的浊度空间分布可视化

var ntu = img2.updateMask(watermask1.mask())
Map.addLayer(ntu.select('NTU'),{min:0,max:5,
palette: '#c2523c,#efba10,#7b7fe1,#1aa87d,#0b2c7b'},'NTU');