Automatic water extent mapping from Sentinel-1 SAR using Otsu’s method

A number of techniques have been proposed for SAR-based water extent mapping. They are supervised and unsupervised classification approaches, thresholding, object-based image analysis, and hybrid approaches. Among them, the most commonly used method in SAR image analysis to distinguish water and non-water areas is image thresholding method where the image is segmented into water/non water class based on the threshold value. Most of the time these threshold values are chosen from the image manually. Finding an optimal threshold value in an automated and data-driven way is crucial for effective mapping.

Otsu’s method is a means of automatically finding an optimal threshold based on the observed distribution of pixel values (Otsu. 1979). The algorithm assumes that the histogram of the image is bimodal (i.e., two classes). In the simplest form, the algorithm returns a single intensity threshold that separates pixels into two classes, foreground and background. This threshold is determined by minimizing intra-class intensity variance, or equivalently, by maximizing inter-class variance. The detailed explanation of the algorithm with example is given in this link.

Interestingly, SAR image exhibits a bimodal pixel distribution due to the low radar returns from water bodies and high radar returns from the surrounding terrain. Hence, we can apply the Otsu’s method on SAR image to find an optimal threshold to delineate water and non-water areas.The following is the sentinel-1 VV image of Mettur Stanley reservoir used for the analysis.A number of techniques have been proposed for SAR-based water extent mapping. They are supervised and unsupervised classification approaches, thresholding, object-based image analysis, and hybrid approaches. Among them, the most commonly used method in SAR image analysis to distinguish water and non-water areas is image thresholding method where the image is segmented into water/non water class based on the threshold value. Most of the time these threshold values are chosen from the image manually. Finding an optimal threshold value in an automated and data-driven way is crucial for effective mapping.

Otsu’s method is a means of automatically finding an optimal threshold based on the observed distribution of pixel values (Otsu. 1979). The algorithm assumes that the histogram of the image is bimodal (i.e., two classes). In the simplest form, the algorithm returns a single intensity threshold that separates pixels into two classes, foreground and background. This threshold is determined by minimizing intra-class intensity variance, or equivalently, by maximizing inter-class variance. The detailed explanation of the algorithm with example is given in this link.

Interestingly, SAR image exhibits a bimodal pixel distribution due to the low radar returns from water bodies and high radar returns from the surrounding terrain. Hence, we can apply the Otsu’s method on SAR image to find an optimal threshold to delineate water and non-water areas.The following is the sentinel-1 VV image of Mettur Stanley reservoir used for the analysis.

Sentinel-1 VV image

Implementation

Below is the earth engine implementation of Otsu’s algorithm with step by step explanation. The script is largely inspired from this article.

First, the histogram for Sentinel-1 VV band is computed using ee.Reducer.histogram() function.

var histogram =S1.select('VV').reduceRegion({
  reducer: ee.Reducer.histogram(255, 2),
  geometry: aoi, 
  scale: 20,
  // maxPixels: 1e13,
  bestEffort: true
});
print(histogram);
print(Chart.image.histogram(S1.select('VV'), aoi, 20));

The histogram shows the bimodal pixel distribution of sentinel-1 VV image.

Once the histogram of the image is computed, we need to use it in Otsu’s algorithm to get the DN value that maximizes interclass variance in sentinel -1 VV band.The following snippet is the functional implementation of Otsu’s algorithm which takes the output of the computed histogram as input.

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

After applying the threshold, the region having slope greater than 7 is masked out, in order to remove the noise due to shadow effect in the SAR image.

var alos = alos.select('AVE_DSM').rename('elev')
var slope = alos.addBands(ee.Terrain.slope(alos), ['slope']);
var water = water .updateMask(slope.select('slope').lt(7))

Below is the final water mask delineated from the SAR Image with a decent accuracy. This appoarch can also be used for flood mapping.You can get the complete code from this link.