# Otsu's method combined with water index to calculate the time-series data?

This question is a continuation of Calculate water frequency in Google Earth Engine based on time series?

The previous question used a set of criteria (the numerical relationship like "larger than", "less than" etc between some indices to extract the water pixels. Then applied a decision tree to get the time series map of water frequency: the map shows the frequency of each pixel appear as water among a collection Landsat images of the whole year.

To get the better result, I slightly revised the code and used a new index: AWEIs, to extract the water pixels.

The difference is, in the previous question, the pixels are filtered by numerical relationship, so the output is simply "Yes" or "No", so the binary map is easy to generated.

However, now I use an index, which is a formula and the output are numbers. So, if I want to get a water-land binary map, I need to find a cutoff value to set the threshold.

To easy demonstrate it, here I attached my revised code, which uses the index of AWEIs:

``````//Choose Area
var region = ee.Geometry.Polygon({
coords: [[[-80.79, 24.98], [-80.08, 26.43], [-80.71, 26.16],[-81.07, 26.20]]],
geodesic: false
});

var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');

var non_mosaic = landsat8 .filterBounds(region) .filterDate('2017-01-01', '2017-12-31');

var getQABits = function(image, start, end, newName) {
// Compute the bits we need to extract.
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += Math.pow(2, i);
}
// Return a single band image of the extracted QA bits, giving the band
// a new name.
return image.select([0], [newName])
.bitwiseAnd(pattern)
.rightShift(start);
};

// Select the QA band.
var QA = image.select(['pixel_qa']);
// Get the internal_cloud_algorithm_flag bit.
// Return an image masking out cloudy areas.
};

// A function to mask out cloud pixels.
var clouds = function(image) {
// Select the QA band.
var QA = image.select(['pixel_qa']);
// Get the internal_cloud_algorithm_flag bit.
return getQABits(QA, 5,5, 'Cloud').eq(0);
// Return an image masking out cloudy areas.
};

var c = clouds(image);
};

//Will be used for frequency calculaton

//display the mosaicked origional image and cloud free image
var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};

//design the decision tree for water and vegetation coverage detection
var decision_tree = function(image){
var mndwi = image.normalizedDifference(['B3', 'B6']);
var evi = image.expression(
'(NIR - RED)/(NIR+6*RED-7.5*BLUE+1)*2.5', {
'BLUE': image.select('B2'),
'RED': image.select('B4'),
'NIR': image.select('B5'),
});
var aweis = image.expression(
'BLUE+2.5*GREEN-1.5*(NIR+SWIR1)-0.25*SWIR2', {
'BLUE' : image.select('B2'),
'GREEN': image.select('B3'),
'SWIR1': image.select('B6'),
'NIR': image.select('B5'),
'SWIR2': image.select('B7'),
});
var msavi = mosaic_free.expression(
'(2*NIR+1-((2*NIR+1)**2-8*(NIR-RED))**(1/2))/2', {
'RED' : mosaic_free.select('B4'),
'NIR': mosaic_free.select('B5'),
});
var ndvi = image.normalizedDifference(['B5', 'B4']);
var lswi = image.normalizedDifference(['B5', 'B6']);
var water = aweis.gte(0).rename('water');
var vegetation = evi.gte(0.1).and (ndvi.gte(0.2)).and (lswi.gt(0)).rename('vegetation');
};

var result = clouds_free.map(decision_tree);
//calculate the frequency
var freq = result.sum().divide(result.count());

var VisParamWar = {"bands":["water"],"min":0,"max":1, palette: ['0000ff', '00ffff','ffffff', 'ffff00', 'ff0000']};
var VisParamVeg = {"bands":["vegetation"],"min":0,"max":1, palette: ['0000ff', '00ffff', 'ffffff', 'ffff00', 'ff0000']};

//display the frequencies
Map.setCenter(-80.702898,25.679568, 8);
``````

The above code works well. As you can see, I selected 0 as the threshold of AWEIs. The pixels with the values equal or larger than 0 are considered as water pixels.

Now, instead of the constant value, I want to use Otsu's method to find the optimized threshold and apply it to the water frequency calculation. Here I share my code as follow:

``````//Choose Area
var region = ee.Geometry.Polygon({
coords: [[[-80.79, 24.98], [-80.08, 26.43], [-80.71, 26.16],[-81.07, 26.20]]],
geodesic: false
});

var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');

var non_mosaic = landsat8 .filterBounds(region) .filterDate('2017-01-01', '2017-12-31');

var getQABits = function(image, start, end, newName) {
// Compute the bits we need to extract.
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += Math.pow(2, i);
}
// Return a single band image of the extracted QA bits, giving the band
// a new name.
return image.select([0], [newName])
.bitwiseAnd(pattern)
.rightShift(start);
};

// Select the QA band.
var QA = image.select(['pixel_qa']);
// Get the internal_cloud_algorithm_flag bit.
// Return an image masking out cloudy areas.
};

// A function to mask out cloud pixels.
var clouds = function(image) {
// Select the QA band.
var QA = image.select(['pixel_qa']);
// Get the internal_cloud_algorithm_flag bit.
return getQABits(QA, 5,5, 'Cloud').eq(0);
// Return an image masking out cloudy areas.
};

var c = clouds(image);
};

//Will be used for frequency calculaton

//display the mosaicked origional image and cloud free image
var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};

// Display the AWEIs.
// Define the visualization parameters.
var vizParams = {bands: ['B6', 'B5', 'B2'], min: 0, max: 0.5,gamma: [0.95, 1.1, 1]};
// Create an AWEIs image, define visualization parameters and display.
var aweis = mosaic_free.expression(
'BLUE+2.5*GREEN-1.5*(NIR+SWIR1)-0.25*SWIR2', {
'BLUE' : mosaic_free.select('B2'),
'GREEN': mosaic_free.select('B3'),
'SWIR1': mosaic_free.select('B6'),
'NIR': mosaic_free.select('B5'),
'SWIR2': mosaic_free.select('B7'),
});

// Compute the histogram of the AWEIs
var histogram = aweis.reduceRegion({
reducer: ee.Reducer.histogram(255, 2)
.combine('mean', null, true)
.combine('variance', null, true),
geometry: region,
scale: 30,
bestEffort: true
});
print(histogram);

// Chart the histogram
print(Chart.image.histogram(aweis, region, 30));

// Return the DN that maximizes interclass variance in B5 (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);
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 threshold1 = otsu(histogram.get('B2_histogram'));
print('Threshold for AWEIs', threshold1);

var classA = aweis.gt(threshold1);