基于GEE平台使用大津法提取水体及Canny边缘检测

 

做水文研究一般第一步都需要进行水体提取，特将代码记录方便后续使用。

首先利用我先前编写好的代码获取遥感图像，对其进行下载，下载完成后再导入GEE中计算MNDWI

var jlj=image//先将我所需区域的TIF影像导入到GEE
var mndwi=jlj.normalizedDifference(['b2','b5'])//计算一下MNDWI
var visual={
  min:0,
  max:1,
  palette:['red','blue']
}
Map.addLayer(mndwi)

接下来使用otsu（大津法）进行一个水体阈值提取。

大津法：假定一张图片共有n个像素，其中灰度值小于阈值的像素为 n1 个，大于等于阈值的像素为 n2 个（ n1 + n2 = n ）。w1 和 w2 表示这两种像素各自的比重。而所有灰度值小于阈值的像素的平均值和方差分别为 μ1 和 σ1，所有灰度值大于等于阈值的像素的平均值和方差分别为 μ2 和 σ2。于是，可以得到：

 

类内差异 = w1(σ1的平方) + w2(σ2的平方)

类间差异 = w1w2(μ1-μ2)^2

 

要找到合适的阈值，让前者（指类内差异）最小或者后者（指类间差异）最大都可以，从计算角度出发，后者要好计算一些。

有关大津法的表述和代码取自下述作者文章：

https://mp.weixin.qq.com/s/vIVB7jfVOf1AuPDRcQrQzQ

接下来是在GEE 编写大津法的代码

// 计算mdwi影像的直方图.

var histogram = mndwi.reduceRegion({

reducer: ee.Reducer.histogram(50)

  .combine('mean', null, true)

  .combine('variance', null, true), 

// geometry: polygon, 

scale: 300,

bestEffort: true

});

//把 otsu 的计算步骤用 map 和 reduce 包裹起来，便于 GEE 调度资源：

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)));

});

// Return the mean value corresponding to the maximum BSS.

return means.sort(bss).get([-1]);

};


//从服务器取出阀值，应用到原图的 NDVI 波段上

var threshold = otsu(histogram.get('nd_histogram')); 

var water = mndwi.lt(threshold);
Map.addLayer(water,{},'WATER')

基于此，水体即可提取出来。 

在将二值化水体图像导出时，遇到的一些问题，做个总结

//导出的一个代码
Export.image.toDrive(
  {
    image:water,
    description:'water_export',
    scale:30,  //这个很关键，当时因为没把尺寸大小写入，导致导出的图像一直有问题
    region:table
  })

一般来说，用大津法提取水体精度还是比较可以的，但是偶尔需要用Canny边缘检测，然后建立缓冲区再进行水体二值化，下面是Canny边缘检测的代码：

  var canny = ee.Algorithms.CannyEdgeDetector(water, 1)//边缘检测，第一个是检测的图像，第二是设置阈值
  canny = canny.updateMask(canny)//用自身屏蔽图像以去除没有边缘的区域
  
Map.addLayer(canny,{},'Canny')  
Export.image.toDrive(
  {
    image:canny,
    description:'water_canny',
    scale:30,
    region:table
  })//边缘检测图像导出