bmax otsu for surface water mapping detection

Dynamic thresholding ftw


Find the scientific description of the algorithm in the paper of Markert et al. (2020). Comparing Sentinel-1 Surface Water Mapping Algorithms and Radiometric Terrain Correction Processing in Southeast Asia Utilizing Google Earth Engine

Figure 2 shows the workflow for the data processing applied.

Copy the code below or use this link


//Script for SAR water segmentation with Edge Otsu Algorithm in Cambodia

//Slope Correction Algorithms
// Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen)
function slopeCorrection(image) { 
  var imgGeom = image.geometry()
  var srtm = ee.Image('JAXA/ALOS/AW3D30/V2_2').select('AVE_DSM').clip(imgGeom) 
  var sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0))

  // Article ( numbers relate to chapters) 
  // 2.1.1 Radar geometry 
  var theta_i = image.select('angle')
  var phi_i = ee.Terrain.aspect(theta_i)
    .reduceRegion(ee.Reducer.mean(), theta_i.get('system:footprint'), 1000)
    .get('aspect')

  // 2.1.2 Terrain geometry
  var alpha_s = ee.Terrain.slope(srtm).select('slope')
  var phi_s = ee.Terrain.aspect(srtm).select('aspect')

  // 2.1.3 Model geometry
  // reduce to 3 angle
  var phi_r = ee.Image.constant(phi_i).subtract(phi_s)

  // convert all to radians
  var phi_rRad = phi_r.multiply(Math.PI / 180);
  var alpha_sRad = alpha_s.multiply(Math.PI / 180);
  var theta_iRad = theta_i.multiply(Math.PI / 180);
  var ninetyRad = ee.Image.constant(90).multiply(Math.PI / 180);

  // slope steepness in range (eq. 2)
  var alpha_r = (alpha_sRad.tan().multiply(phi_rRad.cos())).atan();

  // slope steepness in azimuth (eq 3)
  var alpha_az = (alpha_sRad.tan().multiply(phi_rRad.sin())).atan();

  // local incidence angle (eq. 4)
  var theta_lia = (alpha_az.cos().multiply((theta_iRad.subtract(alpha_r)).cos())).acos();
  var theta_liaDeg = theta_lia.multiply(180 / Math.PI);
  // 2.2 
  // Gamma_nought_flat
  var gamma0 = sigma0Pow.divide(theta_iRad.cos())
  var gamma0dB = ee.Image.constant(10).multiply(gamma0.log10());
  var ratio_1 = gamma0dB.select('VV').subtract(gamma0dB.select('VH'));

  // Volumetric Model
  var nominator = (ninetyRad.subtract(theta_iRad).add(alpha_r)).tan();
  var denominator = (ninetyRad.subtract(theta_iRad)).tan();
  var volModel = (nominator.divide(denominator)).abs();

  // apply model
  var gamma0_Volume = gamma0.divide(volModel);
  var gamma0_VolumeDB = ee.Image.constant(10).multiply(gamma0_Volume.log10());

  // we add a layover/shadow maskto the original implmentation
  // layover, where slope > radar viewing angle 
  var alpha_rDeg = alpha_r.multiply(180 / Math.PI);
  var layover = alpha_rDeg.lt(theta_i);

  // shadow where LIA > 90
  var shadow = theta_liaDeg.lt(85);

  // calculate the ratio for RGB vis
  var ratio = gamma0_VolumeDB.select('VV').subtract(gamma0_VolumeDB.select('VH'));

  var output = gamma0_VolumeDB.addBands(ratio).addBands(alpha_r).addBands(phi_s).addBands(theta_iRad)
    .addBands(layover).addBands(shadow).addBands(gamma0dB).addBands(ratio_1);

  return image.addBands(
    output.select(['VV', 'VH', 'slope_1', 'slope_2'], ['VV', 'VH', 'layover', 'shadow']),
    null,
    true
  ).addBands(image.select("angle"));
}


/***********End of Edge Otsu************/

/**
 * Perona-Malik (anisotropic diffusion) convolution
 * 
 * by Gennadii Donchyts see https://groups.google.com/forum/#!topic/google-earth-engine-developers/a9W0Nlrhoq0
 * I(n+1, i, j) = I(n, i, j) + lambda * (cN * dN(I) + cS * dS(I) + cE * dE(I), cW * dW(I))
 * 
 * I: ee.Image single band, natural units
 * iter: Number of interations to apply filter
 * K: kernal size
 * opt_method: choose method 1 (default) or 2, DETAILS
 * 
 * Returns: single band ee.Image in natural units
 * 
 * Example: image = PeronaMalik(image, 10, 3.5, 1)
 */

function PeronaMalik(I,iter, K, opt_method) {
  iter = iter || 10;
  K = K || 3;
  var method = opt_method || 1;
  // Define kernels
  var dxW = ee.Kernel.fixed(3, 3,
                           [[ 0,  0,  0],
                            [ 1, -1,  0],
                            [ 0,  0,  0]]);
  var dxE = ee.Kernel.fixed(3, 3,
                           [[ 0,  0,  0],
                            [ 0, -1,  1],
                            [ 0,  0,  0]]);
  var dyN = ee.Kernel.fixed(3, 3,
                           [[ 0,  1,  0],
                            [ 0, -1,  0],
                            [ 0,  0,  0]]);
  var dyS = ee.Kernel.fixed(3, 3,
                           [[ 0,  0,  0],
                            [ 0, -1,  0],
                            [ 0,  1,  0]]);
  var lambda = 0.2;
  var k1 = ee.Image(-1.0/K);
  var k2 = ee.Image(K).multiply(ee.Image(K));
  // Convolve
  for(var i = 0; i < iter; i++) {
    var dI_W = I.convolve(dxW);
    var dI_E = I.convolve(dxE);
    var dI_N = I.convolve(dyN);
    var dI_S = I.convolve(dyS);
  // Combine using choosen method
    switch(method) {
      case 1:
        var cW = dI_W.multiply(dI_W).multiply(k1).exp();
        var cE = dI_E.multiply(dI_E).multiply(k1).exp();
        var cN = dI_N.multiply(dI_N).multiply(k1).exp();
        var cS = dI_S.multiply(dI_S).multiply(k1).exp();
        I = I.add(ee.Image(lambda).multiply(cN.multiply(dI_N).add(cS.multiply(dI_S)).add(cE.multiply(dI_E)).add(cW.multiply(dI_W))));
        break;
      case 2:
        var cW = ee.Image(1.0).divide(ee.Image(1.0).add(dI_W.multiply(dI_W).divide(k2)));
        var cE = ee.Image(1.0).divide(ee.Image(1.0).add(dI_E.multiply(dI_E).divide(k2)));
        var cN = ee.Image(1.0).divide(ee.Image(1.0).add(dI_N.multiply(dI_N).divide(k2)));
        var cS = ee.Image(1.0).divide(ee.Image(1.0).add(dI_S.multiply(dI_S).divide(k2)));
        I = I.add(ee.Image(lambda).multiply(cN.multiply(dI_N).add(cS.multiply(dI_S)).add(cE.multiply(dI_E)).add(cW.multiply(dI_W))));
        break;
    }
  }
  return I;
};


// Bmax Otsu Thresholding Algorithm
function bmaxOtsu(img,kwargs){
  

  // define default parameterization for keywords
  var kwargKeys = [];
  for(var key in kwargDefaults) kwargKeys.push( key );
  var params;
  var i,k,v;

  // loop through the keywords and construct ee.Dictionary from them,
  // if the key is defined in the input then pass else use default
  params = ee.Dictionary(kwargs);
  for (i=0;i<kwargKeys.length;i++) {
    k = kwargKeys[i];
    v = kwargDefaults[k];
    params = ee.Dictionary(
      ee.Algorithms.If(params.contains(k),params,params.set(k,v))
    );
  }
  // parse out the keyword parameters as ee objects
  var initialThreshold = ee.Number( params.get('initialThreshold') ),
      reductionScale   = ee.Number( params.get('reductionScale') ),
      gridSize         = ee.Number( params.get('gridSize') ),
      bmaxThresh       = ee.Number( params.get('bmaxThreshold') ),
      maxN             = ee.Number( params.get('maxN') ),
      bandName         = ee.String( params.get('bandName') ),
      maxBuckets       = ee.Number( params.get('maxBuckets') ),
      minBucketWidth   = ee.Number( params.get('minBucketWidth') ),
      maxRaw           = ee.Number( params.get('maxRaw') ), 
      invert           = params.get('invert'),
      verbose          = params.get('verbose').getInfo();
  
  
  // get a list of geometry coordinates from images and find union
  var geom = img.geometry().transform("EPSG:4326",0.1)
  var bounds = geom.bounds()
  var coords = ee.List(bounds.coordinates().get(0));

  // get bounds of the image geometry
  var west = ee.Number(ee.List(coords.get(0)).get(0)),
      south = ee.Number(ee.List(coords.get(0)).get(1)),
      east = ee.Number(ee.List(coords.get(2)).get(0)),
      north = ee.Number(ee.List(coords.get(2)).get(1));


  // round the coordinates out to nearest grid size
  west = west.subtract(west.mod(gridSize));
  south = south.subtract(south.mod(gridSize));
  east = east.add(gridSize.subtract(east.mod(gridSize)));
  north = north.add(gridSize.subtract(north.mod(gridSize)));

  // contruct grid and see if it overlaps with the swath/image geom
  var grid = ee.FeatureCollection(
    ee.List.sequence(south,north.subtract(gridSize),gridSize).map(function(i){
      i = ee.Number(i);
      var out = ee.List.sequence(west,east.subtract(gridSize),gridSize).map(function(j){
        j = ee.Number(j);
        var box = ee.Feature(ee.Geometry.Rectangle(j,i,j.add(gridSize),i.add(gridSize)));
        var out = ee.Algorithms.If(geom.contains(box.geometry()),box,null);
        return ee.Feature(out);
      });
    return out;
    }).flatten()
  );

  Map.addLayer(grid,{},"grid")
  
  var bmax = grid.map(function(feature){
    var segment = img.clip(feature);
    var water = segment.lt(initialThreshold);
    var p1 = water.reduceRegion({
      reducer:ee.Reducer.mean(),
      geometry: feature.geometry(),
      bestEffort: true,
      scale:reductionScale,
    }).get(bandName);
    p1 = ee.Number(ee.Algorithms.If(p1,p1,0.99));
    var p2 = ee.Number(1).subtract(p1);
    
    var m = segment.updateMask(water).rename('m1').addBands(
      segment.updateMask(water.not()).rename('m2'));
      
    var mReduced = m.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: feature.geometry(),
      bestEffort: true,
      scale:reductionScale,
    });
    
    var m1 = ee.Number(mReduced.get('m1'));
    var m2 = ee.Number(mReduced.get('m2'));
    
    m1 = ee.Number(ee.Algorithms.If(m1,m1,-20));
    m2 = ee.Number(ee.Algorithms.If(m2,m2,-5));
    
    var sigmab = p1.multiply(p2.multiply(m1.subtract(m2).pow(2)));
    var sigmat = ee.Number(segment.reduceRegion({
      reducer: ee.Reducer.variance(),
      geometry: feature.geometry(),
      bestEffort: true,
      scale:reductionScale,
    }).get(bandName));
    sigmat = ee.Number(ee.Algorithms.If(sigmat,sigmat,2));
    var bmax = sigmab.divide(sigmat);
    return feature.set({'bmax':bmax});
  }).randomColumn('random');
  
  // filter for nodata
  var bmaxes = bmax.filter(ee.Filter.lt('bmax',1.0))
  var bmaxes = bmaxes.filter(ee.Filter.gt('bmax',bmaxThresh));
  var nBoxes = ee.Number(bmaxes.size());
  

  var randomSelection = maxN.divide(nBoxes);
  var selection = bmaxes.filter(ee.Filter.lt('random',randomSelection));
  Map.addLayer(selection,{},"selection")
  var histogram = ee.Dictionary(img.reduceRegion({
    reducer:ee.Reducer.histogram(maxBuckets, minBucketWidth,maxRaw)
      .combine('mean', null, true).combine('variance', null,true),
    geometry: selection,
    scale: reductionScale,
    maxPixels: 1e13,
    tileScale:16
  }).get(bandName.cat('_histogram')));
  print(histogram)
  var threshold = otsu(histogram);
  
   // segment image and mask 0 values (not water)
  var waterImg = ee.Image(ee.Algorithms.If(invert,img.gt(threshold),img.lt(threshold)));
  
  // if user specifies verbosity then print information
  // useful for inspecting what is going on within the algorithm
  var chart = constructHistChart(histogram,threshold)
      .setOptions({
        title: 'Bmax Search Histogram',
        hAxis: {
          title: 'Values',
        },
        vAxis:{
          title:'Count'
        } 
      });
    print('Algorithm parameters:',params);
    print("Calculated threshold:",threshold);
    print('Thresholding histogram:',chart);
 
 Map.addLayer(waterImg,{palette:"white,blue"},"water")  
 return waterImg;
}

function otsu(histogram) {
  // make sure histogram is an ee.Dictionary object
  histogram = ee.Dictionary(histogram);
  // extract relevant values into arrays
  var counts = ee.Array(histogram.get('histogram'));
  var means  = ee.Array(histogram.get('bucketMeans'));
  // calculate single statistics over arrays
  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);
  // compute between sum of squares, where each mean partitions the data
  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 the mean value corresponding to the maximum BSS
  return means.sort(bss).get([-1]);
}

function constructHistChart(histogram,threshold){
  var counts = ee.List(histogram.get('histogram'));
  var buckets = ee.List(histogram.get('bucketMeans'));
  // construct array for visualization of threshold in chart
  var segment = ee.List.repeat(0, counts.size());
  var maxFrequency   = ee.Number(counts.reduce(ee.Reducer.max()));
  var threshIndex    = buckets.indexOf(threshold);
  segment            = segment.set(threshIndex, maxFrequency);
  var histChart = ui.Chart.array.values(ee.Array.cat([counts, segment], 1), 0, buckets)
    .setSeriesNames(['Values', 'Threshold'])
    .setChartType('ColumnChart');
    
  return histChart;
} 
 
//Declare the period
var start = "2020-10-01";
var end = "2020-10-31";


// define default parameterization for keywords
var kwargDefaults = {
    'initialThreshold': -14,
    'reductionScale': 180,
    'gridSize': 0.1,
    'bmaxThreshold': 0.65,
    'maxN': 100,
    'bandName': "VV",
    'maxBuckets': 255,
    'minBucketWidth': 0.001,
    'maxRaw': 1e6,
    'invert': false,
    'verbose': false
};

// get a few sentinel1 images to run algorithms on
var s1 = ee.ImageCollection("COPERNICUS/S1_GRD")
                    .filterBounds(geometry)
                    .filterDate(start,end);

// apply slope correction and speckle filter
s1 = s1.map(slopeCorrection)
       .map(PeronaMalik);
          

var img = ee.Image(s1.first()).select("VV");

Map.addLayer(img,{min:-25,max:0},"s1 image VV");

var otsu = bmaxOtsu(img,kwargDefaults);

Leave a Reply