指数法和Random Forest实现山东省丰水期地表水体信息

课程大作业而已，仅供参考

随机森林，功能比较全

var table = ee.FeatureCollection("users/gis418670826/shandong_province");
// 山东
var geometry = table
Map.centerObject(geometry, 7);
Map.addLayer(geometry, {}, "geometry", false)

// 根据该水体数据集描述，1没有数据，2不是水体，2季节性水体，3永久水体
// 但是实测数据中没有0和1，只能通过unmask把unmask值填充为1
// 并通过gt(1)，分为水体1和非水体0
var jrc_year = ee.ImageCollection("JRC/GSW1_3/YearlyHistory")
                 .filterDate('2020-01-01', '2021-01-01')
                 .first()
                 .clip(geometry)
                 .select('waterClass')

var jrc_year_Water = jrc_year.reduceRegion({
    crs: 'EPSG:4326',
    reducer: ee.Reducer.count(),
    geometry : geometry,
    scale : 30,
    maxPixels :1e13,
  })
jrc_year =  jrc_year.unmask(ee.Image(1).clip(geometry))
                    .gt(1)
jrc_year = jrc_year.unmask(ee.Image(1).clip(geometry))
Map.addLayer(jrc_year, {min:0,max:1,palette:['#ffffff','red','#0000ff']}, "jrc_year")   

// 水体和非水体各取500个点
var point = jrc_year.stratifiedSample({
  numPoints : 500,
  classBand : "waterClass",
  region : geometry,
  scale :30,
  geometries :true
});
Map.addLayer(point,{}, "point")
print(point)

// landsat8 数据处理
// 去云
function remove_cloud_L8sr(image) {
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;
  var qa = image.select('QA_PIXEL');
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask)
}

// 在已有7波段上再添加mndwi
function computermNdwi(image)
{
  var mndwi = image.normalizedDifference(['SR_B3','SR_B5']).rename('mndwi');
  mndwi.updateMask(mndwi.gt(-1).and(mndwi.lt(1)));
  return image.addBands(mndwi)
}
// LANDSAT/LC08/C02/T1_L2
var landsat8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
                .filterDate('2020-04-01', '2020-11-1')
                .filterBounds(geometry)
                .map(remove_cloud_L8sr)
                .map(computermNdwi)
                .mean()   // 无云数据合成
                .clip(geometry)
Map.addLayer(landsat8,{}, "landsat8")
print(landsat8)

var bands = ['SR_B1','SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7','mndwi']
// 把landsat8波段信息赋予样本点
var training = landsat8.select(bands).sampleRegions({
  collection: point,
  properties: ['waterClass'],
  scale: 30
});

print(training)


var withRandom = training.randomColumn('random');//样本点随机的排列
// 保留一些数据进行测试，以避免模型过度拟合。
var split = 0.7; 
var trainingPartition = withRandom.filter(ee.Filter.lt('random', split));//筛选70%的样本作为训练样本
var testingPartition = withRandom.filter(ee.Filter.gte('random', split));//筛选30%的样本作为测试样本

var classifier = ee.Classifier.smileRandomForest(10).train({
    features:trainingPartition,
    classProperty:"waterClass",
    inputProperties:bands
  })
// print(classifier)

//对Landsat-8进行分类
var class_img = landsat8.select(bands).classify(classifier);
//运用测试样本分类，确定要进行函数运算的数据集以及函数
var test = testingPartition.classify(classifier);
//计算混淆矩阵
var confusionMatrix = test.errorMatrix('waterClass', 'classification');

print('confusionMatrix',confusionMatrix);//面板上显示混淆矩阵
print('overall accuracy', confusionMatrix.accuracy());//面板上显示总体精度
print('kappa accuracy', confusionMatrix.kappa());//面板上显示kappa值


Map.addLayer(class_img, {min: 0, max: 1, palette: ['black','yellow']}, 'result');
print(class_img)
// 计算面积


function computeWaterPixel(image) {
  var Sum = image.reduceRegion({
    crs: 'EPSG:4326',
    reducer: ee.Reducer.count(),
    geometry : geometry,
    scale : 30,
    maxPixels :1e13,
  })
  // print(Sum)
  var mask = image.eq(1)
   var Water = image.updateMask(mask).reduceRegion({
    crs: 'EPSG:4326',
    reducer: ee.Reducer.count(),
    geometry : geometry,
    scale : 30,
    maxPixels :1e13,
  })
  // data.get('sur_refl_b01').getInfo()
  var water_Pixel = Water.get('classification').getInfo()
  var sum_Pixel =  Sum.get('classification').getInfo()
  print("水体所占像素数量",water_Pixel)
 
  print("山东省像素数量",sum_Pixel)
  var area = 155800
  print("已知山东省面积为155800平方千米,水体的面积为", water_Pixel/sum_Pixel*area)
   print("jrc", jrc_year_Water.get('waterClass').getInfo()/sum_Pixel*area)
  
}

computeWaterPixel(class_img)

mndwi + otsu算法

// 山东
var geometry = table
Map.centerObject(geometry, 7);
Map.addLayer(geometry, {}, "geometry")

var jrc_year = ee.ImageCollection("JRC/GSW1_3/YearlyHistory")
                 .filterDate('2020-01-01', '2021-01-01')
                 .first()
                 .clip(geometry)
                 .select('waterClass')
                 .eq(3)//3永久水体 2季节水体
var point = jrc_year.stratifiedSample({
  numPoints : 100,
  classBand : "waterClass",
  region : geometry,
  scale :30,
  geometries :true
});
Map.addLayer(point)
print(point)

// 去云
function remove_cloud_L8sr(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)
}
 

// 数据
var landsat8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
                 .filterDate('2020-05-01', '2020-10-1')
                 .select(['B3','B5','pixel_qa'])
                 .filterBounds(geometry)
                 .map(remove_cloud_L8sr)
                 .mean()
                 .clip(geometry)
print(landsat8)

// computer mndwi for landsat8
function computermNdwi(image)
{
  var mndwi = image.normalizedDifference(['B3','B5']).rename('mndwi');
  return mndwi.updateMask(mndwi.gt(-1).and(mndwi.lt(1)));
}
 
//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]);
}
 
// landsat8 
function computer_water(imageCol){
      var img_mean = ee.Image(computermNdwi(imageCol));
      Map.addLayer(img_mean,{}, "img_mean")
      var histogram = img_mean.reduceRegion({
        reducer: ee.Reducer.histogram(), 
        geometry: geometry, 
        scale: 300,
        maxPixels: 1e13,
        tileScale: 16 //提高速度
      });
      var threshold = otsu(histogram.get("mndwi"));
      // var threshold = -0.23836122208868166
      print(threshold)
      var mask = img_mean.gte(threshold);
      var water = mask.rename("water").setMulti({'time':"2020"});
      print(water)
      // Map.addLayer(ee.Image(water),{}, "water")
      return ee.Image(water)
}
var new_image = ee.List([]);
var new_imag = computer_water(landsat8,new_image);
print(new_imag)
Map.addLayer(new_imag,{min:0,max:1,palette:['black','red']},'water_landsat8_otsu_based')


 

参考： 

在GEE中使用时间序列哨兵1号提取永久水体 - 知乎