Using Otsu's Method to separate water and not-water in GEE

Question

I am trying to separate permanent waterbody and not-water from Sentinel-1 image in earth engine. To do this I am taking a median and mean image for a month from a range of years and then using Otsu's threshold to get a threshold to separate the two classes.

var roi = ee.Geometry.Polygon([[92.214, 24.121],
    [92.214, 25.165],
    [91.006, 25.165],
    [91.006, 24.121],
    [91.006, 24.121]], null, false);
Map.centerObject(roi);

// Function to get all the images for a month for a set year
var jan = function (year) {
    var startDate = ee.Date.fromYMD(year, 1, 1);
    var endDate = ee.Date.fromYMD(year, 1, 31);
    var filtered = collection.filter(ee.Filter.date(startDate, endDate));
    return filtered;
};

var jul = function (year) {
    var startDate = ee.Date.fromYMD(year, 7, 1);
    var endDate = ee.Date.fromYMD(year, 7, 31);
    var filtered = collection.filter(ee.Filter.date(startDate, endDate));
    return filtered;
};

// Otsu's Method

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)));
  });
  
  print(ui.Chart.array.values(ee.Array(bss), 0, means));
  
  // Return the mean value corresponding to the maximum BSS.
  return means.sort(bss).get([-1]);
};

var collection = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.or(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'), ee.Filter.eq('orbitProperties_pass', 'DESCENDING')))
    .filterBounds(roi)
    .select('VV');

var years1 = ee.List.sequence(2015, 2022); // As Sentinel-1 Was launched in April 4, 2014. Some months only have data from 2015
var years2 = ee.List.sequence(2014,2022);

var janCollection = years1.map(jan);
var julCollection = years2.map(jul);

var all_jan = ee.ImageCollection(ee.FeatureCollection(janCollection).flatten());
var all_jul = ee.ImageCollection(ee.FeatureCollection(julCollection).flatten());

var median_jan = all_jan.median().clip(roi);
var median_jul = all_jul.median().clip(roi);

Map.addLayer(median_jan, {min: -25, max: 0}, 'median_jan', 0);
Map.addLayer(median_jul, {min: -25, max: 0}, 'median_jul', 0);

var mean_jan = all_jan.mean().clip(roi);
var mean_jul = all_jul.mean().clip(roi);

Map.addLayer(mean_jan, {min: -25, max: 0}, 'mean_jan', 0);
Map.addLayer(mean_jul, {min: -25, max: 0}, 'mean_jul', 0);

var threshold_func = function(image) {
  var histogram = image.reduceRegion({
  reducer: ee.Reducer.histogram()
      .combine('mean', null, true)
      .combine('variance', null, true), 
  geometry: roi, 
  scale: 10,
  bestEffort: true
  });
  //print('DEBUG HIST', histogram);
  
  var plot_hist = ui.Chart.image.histogram({
  image: image,
  region: roi,
  scale: 10,
  maxPixels: 1e13
  });
  print('Histogram', plot_hist);
  var threshold = otsu(histogram.get('VV_histogram'));
  print('s1 threshold', threshold);
  
  var water = image.lt(threshold).selfMask();
  return water;
};

var jan_water = threshold_func(median_jan);
var jul_water = threshold_func(median_jul);

Map.addLayer(jan_water, {palette: 'blue'}, 'Jan Water', 0);
Map.addLayer(jul_water, {palette: 'blue'}, 'Jul Water', 0);

var jan_water_mean = threshold_func(mean_jan);
var jul_water_mean = threshold_func(mean_jul);

Map.addLayer(jan_water_mean, {palette: 'blue'}, 'Jan Water Mean', 0);
Map.addLayer(jul_water_mean, {palette: 'blue'}, 'Jul Water Mean', 0);

Here's the link to the code too: https://code.earthengine.google.com/efd45e08fed2b79dec04605f79e385fd

The method seems to work pretty well for the month of July, but It's performance is not good in January. It's a lot more noisier. Furthermore, in the mean image there's a line in the middle possibly because that's where two images from sentinel-1 joined in this area of interest.

How do I overcome the problems?

Answer 1

The Sentinel 1 images are noisy on the sides. A quick and dirty fix is to mask based on the angle band. You might have to play with the min/max angle based on your ROI. I'd also switch to a calenderRange() filter. It's simpler, and probably a bit more performant.

var maskBorder = function (image) {
  var minAngle = 30.88
  var maxAngle = 45.35
  var angle = image.select('angle')
  var sideMask = angle.gt(minAngle).and(angle.lt(maxAngle))
  return image.updateMask(sideMask)
}

var collection = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.or(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'), ee.Filter.eq('orbitProperties_pass', 'DESCENDING')))
    .filter(ee.Filter.inList('orbitProperties_pass', ['ASCENDING', 'DESCENDING']))
    .filterBounds(roi)

var all_jan = collection
  .filter(ee.Filter.calendarRange(1, 1, 'month'))
  .map(maskBorder)
  .select('VV')
var all_jul = collection
  .filter(ee.Filter.calendarRange(7, 7, 'month'))
  .map(maskBorder)
  .select('VV')

https://code.earthengine.google.com/05e317c5473ec353a5fccceffe664d2c