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

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)

// slope steepness in azimuth (eq 3)

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

// Volumetric Model
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());

// 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

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

output.select(['VV', 'VH', 'slope_1', 'slope_2'], ['VV', 'VH', 'layover', 'shadow']),
null,
true
}

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

/**
* Perona-Malik (anisotropic diffusion) convolution
*
* 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();
break;
case 2:
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));

// 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 out = ee.Algorithms.If(geom.contains(box.geometry()),box,null);
return ee.Feature(out);
});
return out;
}).flatten()
);

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

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