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; }; //Edge Otsu Algorithms function edgeOtsu(img,kwargs) { var geom = img.geometry() // get list of band names used later var bandList = img.bandNames(); 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)) ); } // parameters for all methods var initialThreshold = ee.Number( params.get('initialThreshold') ), reductionScale = ee.Number( params.get('reductionScale') ), smoothing = ee.Number( params.get('smoothing') ), bandName = ee.String( params.get('bandName') ), connectedPixels = ee.Number( params.get('connectedPixels') ), edgeLength = ee.Number( params.get('edgeLength') ), smoothEdges = ee.Number( params.get('smoothEdges') ), cannyThreshold = ee.Number( params.get('cannyThreshold') ), cannySigma = ee.Number( params.get('cannySigma') ), cannyLt = ee.Number( params.get('cannyLt') ), 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 preliminary water var binary = img.lt(initialThreshold).rename('binary'); Map.addLayer(binary,{min:0,max:1},"binary threshold") // get canny edges var canny = ee.Algorithms.CannyEdgeDetector(binary, cannyThreshold, cannySigma); Map.addLayer(canny,{},"Canny edge detection") // process canny edges var connected = canny.updateMask(canny).lt(cannyLt).connectedPixelCount(connectedPixels, true); var edges = connected.gte(edgeLength); edges = edges.updateMask(edges); Map.addLayer(edges,{},"edges") var edgeBuffer = edges.focal_max(smoothEdges, 'square', 'meters'); Map.addLayer(edgeBuffer,{},"Edge buffer") // get histogram for Otsu var histogram_image = img.updateMask(edgeBuffer); // histogram_image = histogram_image.clip(geometry2) var histogram = ee.Dictionary(histogram_image.reduceRegion({ reducer:ee.Reducer.histogram(maxBuckets, minBucketWidth,maxRaw) .combine('mean', null, true).combine('variance', null,true), geometry: geom, scale: reductionScale, maxPixels: 1e13, tileScale:16 }).get(bandName.cat('_histogram'))); var threshold = otsu(histogram); var chart = constructHistChart(histogram,threshold) .setOptions({ title: 'Edge Search Histogram', hAxis: { title: 'Values', }, vAxis:{ title:'Count' } }); print('Algorithm parameters:',params); print("Calculated threshold:",threshold); print('Thresholding histogram:',chart); // segment image and mask 0 values (not water) var waterImg = ee.Image(ee.Algorithms.If(invert,img.gt(threshold),img.lt(threshold))); Map.addLayer(waterImg,{palette:"white,blue"},"water image") 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, 'smoothing': 100, 'bandName': "VV", 'connectedPixels': 100, 'edgeLength': 20, 'smoothEdges': 20, 'cannyThreshold': 1, 'cannySigma': 1, 'cannyLt': 0.05, '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 = edgeOtsu(img,kwargDefaults);