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
});
//Add image collection
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
// Fmask classification values var FMASK_CLEAR_GROUND = 0; var FMASK_WATER = 2; var FMASK_CLOUD_SHADOW = 3; var FMASK_SNOW = 4; var FMASK_CLOUD = 5;
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);
};
// A function to mask out cloudy shadow pixels.
var cloud_shadows = function(image) {
// Select the QA band.
var QA = image.select(['pixel_qa']);
// Get the internal_cloud_algorithm_flag bit.
return getQABits(QA, 3,3, 'Cloud_shadows').eq(0);
// 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 maskClouds = function(image) {
var cs = cloud_shadows(image);
var c = clouds(image);
image = image.updateMask(cs);
return image.updateMask(c);
};
//Will be used for frequency calculaton
var clouds_free = non_mosaic.map(maskClouds);
//display the mosaicked origional image and cloud free image
var mosaic_free = non_mosaic.map(maskClouds).median();
var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};
Map.addLayer(non_mosaic, visParams, 'With clouds');
Map.addLayer(mosaic_free, visParams, 'Cloud free');
//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');
return water.addBands(vegetation).copyProperties(image);
};
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);
Map.addLayer(freq,VisParamWar,'Water Frequency');
Map.addLayer(freq,VisParamVeg,'Vegetation Frequency');
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
});
//Add image collection
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
// Fmask classification values var FMASK_CLEAR_GROUND = 0; var FMASK_WATER = 2; var FMASK_CLOUD_SHADOW = 3; var FMASK_SNOW = 4; var FMASK_CLOUD = 5;
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);
};
// A function to mask out cloudy shadow pixels.
var cloud_shadows = function(image) {
// Select the QA band.
var QA = image.select(['pixel_qa']);
// Get the internal_cloud_algorithm_flag bit.
return getQABits(QA, 3,3, 'Cloud_shadows').eq(0);
// 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 maskClouds = function(image) {
var cs = cloud_shadows(image);
var c = clouds(image);
image = image.updateMask(cs);
return image.updateMask(c);
};
//Will be used for frequency calculaton
var clouds_free = non_mosaic.map(maskClouds);
//display the mosaicked origional image and cloud free image
var mosaic_free = non_mosaic.map(maskClouds).median();
var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};
Map.addLayer(non_mosaic, visParams, 'With clouds');
Map.addLayer(mosaic_free, visParams, 'Cloud free');
// 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);
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 threshold1 = otsu(histogram.get('B2_histogram'));
print('Threshold for AWEIs', threshold1);
var classA = aweis.gt(threshold1);
Map.addLayer(classA.mask(classA), {palette: 'blue'}, 'Dec 31st AWEIs applied OTSU');
So far, the above code works well. It gives me a map of water pixels of December 31st 2017, the threshold was given by Otsu's method.
Now I wish to insert this Otsu's method into the decision tree so that I can get the time-series water frequency map based on the whole year's Landsat image collection, just like the output given by the first piece of code in this post. But I get stuck now.
Could anyone gives me some hints?
