Calculate water frequency in Google Earth Engine based on time series?

Question

I have the Landsat 8 images cover southern Florida during 2017 (from Jan 1st to Dec 31st), with the cloud, scan edges and other bad observation pixels eliminated.

I have also calculated several indices like NDVI, EVI, LSVI and mNDWI. For each pixel, I filtered, selected and rendered different colors so that the map coverage area could be divided into different areas based on its land cover information.

Based on those indices and the generated maps, I further examine each pixel by comparing the mathematical relationships (equal, larger than, less than etc.) between different indices for that pixel's area. In this way, I further filtered the pixels and generated the maps for water and vegetation areas. For example, in the water map, the water pixels have value of 1, while the non-water pixels have value of 0, the map is hence divided into two areas by the value of the pixel.

The above-mentioned map is just for a sea-land division on a single day. What I want is, since I have a whole year's image collection, I want to summarize it and work out a statistics map. That is, a pixel is 1 if it is water, while it is 0 if it is land (non-water). For different images in a single year's collection, that pixel can sometimes be 0 and sometimes be 1, so we can get a year-long frequency value, which is between 0 and 1. I want to visualize it on a map, with the gradient color indicates the frequency.

Here's my code for better reference, the layers named 'water' and 'vegetation' are the layers that I want to generate and visualize the frequencies.

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

// Function to cloud mask from the pixel_qa band of Landsat 8 SR data.
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;

  // Get the pixel QA band.
  var qa = image.select('pixel_qa');

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));

  // Return the masked image, scaled to reflectance, without the QA bands.
  return image.updateMask(mask).divide(10000)
      .select("B[0-9]*")
      .copyProperties(image, ["system:time_start"]);
}

// Map the function over one year of data.
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
    .filterBounds(region)
    .filterDate('2017-01-01', '2017-12-31')
    .map(maskL8sr)

var composite = collection.median();

// Display the results.
Map.addLayer(composite, {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.3}, "Original");

// Display the mNDWI.
Map.setCenter(-80.7200, 25.7874, 8);
// Define the visualization parameters.
var vizParams = {bands: ['B6', 'B5', 'B2'], min: 0, max: 0.5,gamma: [0.95, 1.1, 1]};
// Create an mNDWI image, define visualization parameters and display.
var mndwi = composite.normalizedDifference(['B3', 'B6']);
var mndwiViz = {min: 0.5, max: 1, palette: ['00FFFF', '0000FF']};
Map.addLayer(mndwi, mndwiViz, 'mNDWI', false);

// Display the EVI.
// Define the visualization parameters.
var vizParams = {bands: ['B6', 'B5', 'B2'], min: 0, max: 0.5,gamma: [0.95, 1.1, 1]};
// Create an EVI image, define visualization parameters and display.
var evi = composite.expression(
    '(NIR - RED)/(NIR+6*RED-7.5*BLUE+1)*2.5', {
      'BLUE': composite.select('B2'),
      'RED': composite.select('B4'),
      'NIR': composite.select('B5'),
});
var eviViz = {min: 0.5, max: 1, palette: ['00FFFF', '0000FF']};
Map.addLayer(evi, eviViz, 'EVI', false);

// Display the NDVI.
// Define the visualization parameters.
var vizParams = {bands: ['B6', 'B5', 'B2'], min: 0, max: 0.5,gamma: [0.95, 1.1, 1]};
// Create an NDVI image, define visualization parameters and display.
var ndvi = composite.normalizedDifference(['B5', 'B4']);
var ndviViz = {min: 0.5, max: 1, palette: ['00FFFF', '0000FF']};
Map.addLayer(ndvi, ndviViz, 'NDVI', false);

// Display the LSWI.
// Define the visualization parameters.
var vizParams = {bands: ['B6', 'B5', 'B2'], min: 0, max: 0.5,gamma: [0.95, 1.1, 1]};
// Create an LSWI image, define visualization parameters and display.
var lswi = composite.normalizedDifference(['B5', 'B6']);
var lswiViz = {min: 0.5, max: 1, palette: ['00FFFF', '0000FF']};
Map.addLayer(ndvi, ndviViz, 'LSWI', false);

var water = evi.lt(0.1).and(mndwi.gt(evi).or(mndwi.gt(ndvi)))
var waterViz = {min: 0, max: 1, palette: ['ffff00', '0000FF']};
Map.addLayer(water, waterViz, 'water', false);

var vegetation = evi.gte(0.1).and (ndvi.gte(0.2)).and (lswi.gt(0))
var vegetationViz = {min: 0, max: 1, palette: ['ffffff00', '1aff1a']};
Map.addLayer(vegetation, vegetationViz, 'vegetation', false);

Answer 1

You can use sum() reducer to compute frequency. I took your code to reduce it and change the approach. As I said to you in another post, reducers need to be applied as the last process.

The process itself is just:

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 ndvi = image.normalizedDifference(['B5', 'B4']);
  var lswi = image.normalizedDifference(['B5', 'B6']);
  var water = evi.lt(0.1).and(mndwi.gt(evi).or(mndwi.gt(ndvi))).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);
var freq = result.sum().divide(result.count());

var VisParamWar = {"bands":["water"],"min":0,"max":1};
var VisParamVeg = {"bands":["vegetation"],"min":0,"max":1};

Map.addLayer(freq,VisParamWar,'Water freq');
Map.addLayer(freq,VisParamVeg,'Vegetation freq');

Link to full code: https://code.earthengine.google.com/c9272582da740ca0a4b247bcc4469970