Discontinuity error in Sentinel-1 water map

Question

I tried to generate the water occurrence map from Sentinel-1 to compare with the JRC global dataset. Here is my problem, after calculating the Occurrence of water through sum() and divided it for each image by using size(). My result have some discontinuity like the image.

Does somebody have any idea to fix it?

Here is my code: https://code.earthengine.google.com/e000026102c2bc09376c35b9ea8cde3b?hideCode=true

 var aoi = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level0").filter(ee.Filter.eq('ADM0_NAME', 'Pakistan'));
    Map.addLayer(aoi, {color: 'black'}, 'Study Area',1);
    var aoi = aoi.geometry().bounds();
    // Zoom to regions of interest
    Map.centerObject(aoi);
    
    // import sentinel 1 and filter data series
    var s1 =  ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
    .filterBounds(aoi)
    //.filterBounds(Map.getBounds(true))
    .filterDate('2022-1-01','2022-12-31')
    //.filter(ee.Filter.contains({leftField: ".geo", rightValue: aoi})) // Filter partial S1-Images of AOI
    .map(function(image){return image.clip(Map.getBounds(true))})
    .map(function(image){return image.addBands(image.select('VV').focal_median(parseFloat('50'),'circle','meters').rename('VV_smoothed'))}); // Smooth S1-Images
    
    print(s1);
    
    // Return the DN that maximizes interclass variance in S1-band (in the region).
    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]);
    };
    
    // return image with water mask as additional band
    var add_waterMask = function(image){
      // Compute histogram
      var histogram = image.select('VV').reduceRegion({
        reducer: ee.Reducer.histogram(255, 2)
          .combine('mean', null, true)
          .combine('variance', null, true), 
        geometry: aoi, 
        scale: 10,
        bestEffort: true
      });
      // Calculate threshold via function otsu (see before)
      var threshold = otsu(histogram.get('VV_histogram'));
      
      // get watermask
      var waterMask = image.select('VV_smoothed').lt(threshold).rename('waterMask');
      waterMask = waterMask.updateMask(waterMask); //Remove all pixels equal to 0
      return image.addBands(waterMask);
    };
    
    s1 = s1.map(add_waterMask);
    
    //Calculating water occurrence
    var min_occurence = 10;
    var water_sum = s1.select('waterMask').reduce(ee.Reducer.mean());
    var water_frequency = water_sum.divide(s1.select('waterMask').count()).multiply(100);
    var water_frequency_masked = water_frequency.updateMask(water_frequency.gt(min_occurence));
    
    //Add color bar
    //base code adapted from: 
    //https://gis.stackexchange.com/questions/290713/adding-map-key-to-map-or-console-in-google-earth-engine
    //https://code.earthengine.google.com/9f890c110e98fa3391480543009c8028
    
    function ColorBar(palette) {
      return ui.Thumbnail({
        image: ee.Image.pixelLonLat().select(0),
        params: {
          bbox: [0, 0, 1, 0.1],
          dimensions: '300x15',
          format: 'png',
          min: 0,
          max: 1,
          palette: palette,
        },
        style: {stretch: 'horizontal', margin: '0px 22px'},
      });
    }
    function makeLegend(lowLine, midLine, highLine,lowText, midText, highText, palette) {
      var  labelheader = ui.Label('Water occurrence during investigation period',{margin: '5px 17px', textAlign: 'center', stretch: 'horizontal', fontWeight: 'bold'});
      var labelLines = ui.Panel(
          [
            ui.Label(lowLine, {margin: '-4px 21px'}),
            ui.Label(midLine, {margin: '-4px 0px', textAlign: 'center', stretch: 'horizontal'}),
            ui.Label(highLine, {margin: '-4px 21px'})
          ],
          ui.Panel.Layout.flow('horizontal'));
          var labelPanel = ui.Panel(
          [
            ui.Label(lowText, {margin: '0px 14.5px'}),
            ui.Label(midText, {margin: '0px 0px', textAlign: 'center', stretch: 'horizontal'}),
            ui.Label(highText, {margin: '0px 1px'})
          ],
          ui.Panel.Layout.flow('horizontal'));
        return ui.Panel({
          widgets: [labelheader, ColorBar(palette), labelLines, labelPanel], 
          style: {position:'bottom-left'}});
    }
    Map.add(makeLegend('|', '|', '|', "0 %", '50 %', '100%', ['orange','yellow','lightblue','darkblue']));
    
    // time-lapse animation
    var timelapse = {
      bands: ["VV","VV","VV"],
      region: aoi,
      min: -20,
      max: 0,
      framesPerSecond: 5};
      
    var animation = ui.Thumbnail({
      image: s1,
      params: timelapse,
      style: {
        position: 'bottom-left',
        width: '360px',
      }});
    
    //Add layers ans animation to map
    Map.addLayer(s1.median().clip(aoi),{bands: ['VV','VV','VV'],min: -20,max: 0,},'S1-image [median]');
    Map.addLayer(water_frequency_masked.clip(aoi),{min:0,max:100,palette:['orange','yellow','lightblue','darkblue']},'Percentage of annual water occurence');

Answer 1

First, unrelated to your "discontinuities": You are calculating the percentage by taking the mean and dividing it by the count then multiplying by 100. You will want to skip the division, and make sure not to mask out non-water pixels.

var add_waterMask = function(image) {
  ...
  var waterMask = image.select('VV_smoothed').lt(threshold).rename('waterMask')
  return image.addBands(waterMask)
}

s1 = s1.map(add_waterMask);
var water_mean = s1.select('waterMask').reduce(ee.Reducer.mean())
var water_frequency = water_mean.multiply(100)

You are also using VV_smoothed for the water mask, but calculating your threshold using VV. You will want to be consistent there.

Then, about that threshold. You calculate a new threshold for every individual image. That will lead to border effects. Perhaps you can create a median composite of your collection and calculate the threshold based on that instead?

You probably also want to do a bit more processing of the imagery. At the very least convert to gamma0, maybe do some masking of the bad imagery at the borders of the images.

var gamma0 = image.expression('i/(cos(angle * pi / 180))', {
  'i': image.select(['VV', 'VH']),
  'angle': image.select('angle'),
  'pi': Math.PI
})

var sideMask = angle.gt(30.88).and(angle.lt(45.35)) // Tweak angles for your AOI

Another cause of your "discontinuities" is probably coming from the fact that you have more imagery where orbits overlap. You can get rid of by this by masking based on the dominant orbit for each pixel.

The below script implements these improvements. It takes time to calculate the threshold and you still have some issues, but it's better.

https://code.earthengine.google.com/d431ed691d76cf52de70ecdd4c9cccec