// Example workflow from Boothroyd et al., Science of the Total Environment (post-peer review) manuscript

// Code developed for extraction of the active river channel, active channel width quantification and calculation of similarity coefficients at bridges infrastructure in the Philippines using Landsat imagery.

// Code uses RivWidthCloud for quantification of centreline and active river width: Yang, X., Pavelsky, T. M., Allen, G. H. and Donchyts, G. 2019. RivWidthCloud: An Automated Google Earth Engine Algorithm for River Width Extraction From Remotely Sensed Imagery. IEEE Geoscience and Remote Sensing Letters. 17(2), pp. 217-221. https://doi.org/10.1109/LGRS.2019.2920225  
// Code uses water classification from Zou, Z., Xiao, X., Dong, J., Qin, Y., Doughty, R.B., Menarguez, M.A., Zhang, G. and Wang, J. 2018. Divergent trends of open-surface water body area in the contiguous United States from 1984 to 2016. Proceedings of the National Academy of Sciences. 115(15), 3810-3815. https://doi.org/10.1073/pnas.1719275115   

// Code writes binary active channel image (.tiff) and RivWidthCloud data (.csv) to Google Drive
// Code outputs similarity coefficient data to GEE console for analysis.

// T1, T2, T3 and T4 refer to time periods of analysis (1988-89 (T1), 1998-99 (T2), 2008-09 (T3) and 2018-19 (T4))


// User-defined inputs 
var lat = 120.6674653;
var long = 18.13764236;

var bridge_name = ('Sarrat_Bridge');
var buffer_length = (4500);

var sDate_T1 = "1988-01-01"; 
var eDate_T1 = "1990-01-01"; 

var sDate_T2 = "1998-01-01"; 
var eDate_T2 = "2000-01-01"; 

var sDate_T3 = "2008-01-01"; 
var eDate_T3 = "2010-01-01"; 

var sDate_T4 = "2018-01-01"; 
var eDate_T4 = "2020-01-01"; 

var scale = 30;
var cleaning_pixels = 100;
var mndwi_param = -0.40;
var ndvi_param = 0.20;
var kernel = ee.Kernel.circle({radius: 3});

var folder = ('GEE_Bridges_45km');
var T1_widths_name = (bridge_name+'_T1_'+'_Widths');
var T1_binary_name = (bridge_name+'_T1_'+'_Binary_Image');
var T2_widths_name = (bridge_name+'_T2_'+'_Widths');
var T2_binary_name = (bridge_name+'_T2_'+'_Binary_Image');
var T3_widths_name = (bridge_name+'_T3_'+'_Widths');
var T3_binary_name = (bridge_name+'_T3_'+'_Binary_Image');
var T4_widths_name = (bridge_name+'_T4_'+'_Widths');
var T4_binary_name = (bridge_name+'_T4_'+'_Binary_Image');


// Runs GEE code with a standard circular buffer around the bridge point: 
var geometry = ee.Geometry.Point(lat,long);
var roi = geometry.buffer(buffer_length);
Map.centerObject(roi,12)

// Runs GEE code with a clipped version of the buffer (e.g. for locations near the ocean).
// User input needed: draw an additional geometry called roi_clip
// var diff1 = roi.difference(roi_clip, ee.ErrorMargin(1));
// Map.addLayer(diff1, {color: 'FFFF00'}, 'diff1');
// var roi = diff1;


// Standardise band names, merge Landsat data/
var bn8 = ['B1', 'B2', 'B3', 'B4', 'B6', 'pixel_qa', 'B5', 'B7'];
var bn7 = ['B1', 'B1', 'B2', 'B3', 'B5', 'pixel_qa', 'B4', 'B7'];
var bn5 = ['B1', 'B1', 'B2', 'B3', 'B5', 'pixel_qa', 'B4', 'B7'];
var bns = ['uBlue', 'Blue', 'Green', 'Red', 'Swir1', 'BQA', 'Nir', 'Swir2'];

var ls5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR").select(bn5, bns);
var ls7 = (ee.ImageCollection("LANDSAT/LE07/C01/T1_SR")
  .filterDate('1999-04-15', '2003-05-30') // exclude LE7 images that affected by failure of Scan Line Corrector (SLC) 
  .select(bn7, bns));
var ls8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR").select(bn8, bns);
var merged = ls5.merge(ls7).merge(ls8);


// Functions for active river belt classification
var Ndvi = function(image) {
  // calculate normalized difference vegetation index
  var ndvi = image.normalizedDifference(['Nir', 'Red']).rename('ndvi');
  return(ndvi);
};

var Lswi = function(image) {
  // calculate land surface water index
  var lswi = image.normalizedDifference(['Nir', 'Swir1']).rename('lswi');
 return(lswi);
};

var Mndwi = function(image) {
  // calculate modified normalized difference water index
  var mndwi = image.normalizedDifference(['Green', 'Swir1']).rename('mndwi');
  return(mndwi);
};

var Evi = function(image) {
  // calculate the enhanced vegetation index
  var evi = image.expression('2.5 * (Nir - Red) / (1 + Nir + 6 * Red - 7.5 * Blue)', {
    'Nir': image.select(['Nir']),
    'Red': image.select(['Red']),
    'Blue': image.select(['Blue'])
    });
  return(evi.rename(['evi']));
};


// Sets visualisation parameters
var params_false = {crs: 'EPSG:4326', region: roi, min: 0.0, max: 0.3, bands: ["Swir1", "Red", "Green"], dimensions: 1000};




// T1

// Filter date range, roi and apply simple cloud processing
var imgCol = merged.filterDate(sDate_T1, eDate_T1) 
                 .filterBounds(roi)
                 .map(function(image) { 
                    var cloudShadowBitMask = 1 << 3; 
                    var cloudsBitMask = 1 << 5; 
                    var qa = image.select('BQA'); 
                    var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) 
                                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0)); 
                    return image.updateMask(mask).multiply(0.0001).clip(roi)
                   }); 
                   

// Define and rename quantiles of interest, display false colour
var bnp50 = ['uBlue_p50', 'Blue_p50', 'Green_p50', 'Red_p50', 'Swir1_p50', 'BQA_p50', 'Nir_p50', 'Swir2_p50'];
var p50 = imgCol.reduce(ee.Reducer.percentile([50])).select(bnp50, bns);
Map.addLayer(p50.clip(roi), params_false, 'T1 (1988 - 1990) false colour composite image', true);

// Apply to each percentile
var mndwi_p50 = Mndwi(p50);
var ndvi_p50 = Ndvi(p50);
var evi_p50 = Evi(p50);
var lswi_p50 = Lswi(p50);

// Water classification from Zou 2018:
var water_p50 = (mndwi_p50.gt(ndvi_p50).or(mndwi_p50.gt(evi_p50))).and(evi_p50.lt(0.1));

// Active channel classification:
var activebelt_p50 = (mndwi_p50.gte(mndwi_param)).and(ndvi_p50.lte(ndvi_param));

// Binary image for active channel classication 
var active_p50 = (water_p50).or(activebelt_p50).rename('active');

// Morphological closing, noise removal and area (# pix and km2) calculation 
var T1_orig = active_p50.select(['active']);

var T1_noise_removal = T1_orig
          .updateMask(T1_orig.connectedPixelCount(cleaning_pixels, false).gte(cleaning_pixels))
          .unmask()
          .clip(roi)
          .reproject('EPSG:4326', null, scale);

var T1_noise_removal_active_only = T1_noise_removal.updateMask(T1_noise_removal.gt(0)).reproject('EPSG:4326', null, scale);

var T1_noise_removal_closed = T1_noise_removal
          .focal_max({kernel: kernel, iterations: 1})
          .focal_min({kernel: kernel, iterations: 1})
          .reproject('EPSG:4326', null, scale)
          .clip(roi);

var T1_noise_removal_closed_active_only = T1_noise_removal_closed.updateMask(T1_noise_removal_closed.gt(0)).reproject('EPSG:4326', null, scale);

var T1_noise_removal_closed_active_stats = T1_noise_removal_closed_active_only.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: roi,
  scale: 30,
}).get('active');

print(T1_noise_removal_closed_active_stats,'T1 final area (# of pix)');

var T1_area_km = T1_noise_removal_closed_active_only.multiply(ee.Image.pixelArea()).divide(1000*1000);

var T1_stats_km = T1_area_km.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
});

var T1_stats_km_value = T1_stats_km.get('active');
print(T1_stats_km_value,'T1 final area (km^2)');

var T1_final = T1_noise_removal_closed.reproject('EPSG:4326', null, scale);
Map.addLayer(T1_final, {}, 'T1 (1988 - 1990) active channel binary mask');




// T2

// Filter date range, roi and apply simple cloud processing
var imgCol = merged.filterDate(sDate_T2, eDate_T2) 
                 .filterBounds(roi)
                 .map(function(image) { 
                    var cloudShadowBitMask = 1 << 3; 
                    var cloudsBitMask = 1 << 5; 
                    var qa = image.select('BQA'); 
                    var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) 
                                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0)); 
                    return image.updateMask(mask).multiply(0.0001).clip(roi)
                   }); 
                   

// Define and rename quantiles of interest, display false colour
var bnp50 = ['uBlue_p50', 'Blue_p50', 'Green_p50', 'Red_p50', 'Swir1_p50', 'BQA_p50', 'Nir_p50', 'Swir2_p50'];
var p50 = imgCol.reduce(ee.Reducer.percentile([50])).select(bnp50, bns);
Map.addLayer(p50.clip(roi), params_false, 'T2 (1998 - 2000) false colour composite image', true);

// Apply to each percentile
var mndwi_p50 = Mndwi(p50);
var ndvi_p50 = Ndvi(p50);
var evi_p50 = Evi(p50);
var lswi_p50 = Lswi(p50);

// Water classification from Zou 2018:
var water_p50 = (mndwi_p50.gt(ndvi_p50).or(mndwi_p50.gt(evi_p50))).and(evi_p50.lt(0.1));

// Active channel classification:
var activebelt_p50 = (mndwi_p50.gte(mndwi_param)).and(ndvi_p50.lte(ndvi_param));

// Binary image for active channel classication 
var active_p50 = (water_p50).or(activebelt_p50).rename('active');

// Morphological closing, noise removal and area (# pix and km2) calculation 

var T2_orig = active_p50.select(['active']);

var T2_noise_removal = T2_orig
          .updateMask(T2_orig.connectedPixelCount(cleaning_pixels, false).gte(cleaning_pixels))
          .unmask()
          .clip(roi)
          .reproject('EPSG:4326', null, scale);

var T2_noise_removal_active_only = T2_noise_removal.updateMask(T2_noise_removal.gt(0)).reproject('EPSG:4326', null, scale);

var T2_noise_removal_closed = T2_noise_removal
          .focal_max({kernel: kernel, iterations: 1})
          .focal_min({kernel: kernel, iterations: 1})
          .reproject('EPSG:4326', null, scale)
          .clip(roi);

var T2_noise_removal_closed_active_only = T2_noise_removal_closed.updateMask(T2_noise_removal_closed.gt(0)).reproject('EPSG:4326', null, scale);

var T2_noise_removal_closed_active_stats = T2_noise_removal_closed_active_only.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: roi,
  scale: 30,
}).get('active');

print(T2_noise_removal_closed_active_stats,'T2 final area (# of pix)');

var T2_area_km = T2_noise_removal_closed_active_only.multiply(ee.Image.pixelArea()).divide(1000*1000);

var T2_stats_km = T2_area_km.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
});

var T2_stats_km_value = T2_stats_km.get('active');
print(T2_stats_km_value,'T2 final area (km^2)');

var T2_final = T2_noise_removal_closed.reproject('EPSG:4326', null, scale);
Map.addLayer(T2_final, {}, 'T2 (1998 - 2000) active channel binary mask');




// T3

// Filter date range, roi and apply simple cloud processing
var imgCol = merged.filterDate(sDate_T3, eDate_T3) 
                 .filterBounds(roi)
                 .map(function(image) { 
                    var cloudShadowBitMask = 1 << 3; 
                    var cloudsBitMask = 1 << 5; 
                    var qa = image.select('BQA'); 
                    var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) 
                                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0)); 
                    return image.updateMask(mask).multiply(0.0001).clip(roi)
                   }); 
                   

// Define and rename quantiles of interest, display false colour
var bnp50 = ['uBlue_p50', 'Blue_p50', 'Green_p50', 'Red_p50', 'Swir1_p50', 'BQA_p50', 'Nir_p50', 'Swir2_p50'];
var p50 = imgCol.reduce(ee.Reducer.percentile([50])).select(bnp50, bns);
Map.addLayer(p50.clip(roi), params_false, 'T3 (2008 - 2010) false colour composite image', true);

// Apply to each percentile
var mndwi_p50 = Mndwi(p50);
var ndvi_p50 = Ndvi(p50);
var evi_p50 = Evi(p50);
var lswi_p50 = Lswi(p50);

// Water classification from Zou 2018:
var water_p50 = (mndwi_p50.gt(ndvi_p50).or(mndwi_p50.gt(evi_p50))).and(evi_p50.lt(0.1));

// Active channel classification:
var activebelt_p50 = (mndwi_p50.gte(mndwi_param)).and(ndvi_p50.lte(ndvi_param));

// Binary image for active channel classication 
var active_p50 = (water_p50).or(activebelt_p50).rename('active');

// Morphological closing, noise removal and area (# pix and km2) calculation 

var T3_orig = active_p50.select(['active']);

var T3_noise_removal = T3_orig
          .updateMask(T3_orig.connectedPixelCount(cleaning_pixels, false).gte(cleaning_pixels))
          .unmask()
          .clip(roi)
          .reproject('EPSG:4326', null, scale);

var T3_noise_removal_active_only = T3_noise_removal.updateMask(T3_noise_removal.gt(0)).reproject('EPSG:4326', null, scale);

var T3_noise_removal_closed = T3_noise_removal
          .focal_max({kernel: kernel, iterations: 1})
          .focal_min({kernel: kernel, iterations: 1})
          .reproject('EPSG:4326', null, scale)
          .clip(roi);

var T3_noise_removal_closed_active_only = T3_noise_removal_closed.updateMask(T3_noise_removal_closed.gt(0)).reproject('EPSG:4326', null, scale);

var T3_noise_removal_closed_active_stats = T3_noise_removal_closed_active_only.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: roi,
  scale: 30,
}).get('active');

print(T3_noise_removal_closed_active_stats,'T3 final area (# of pix)');

var T3_area_km = T3_noise_removal_closed_active_only.multiply(ee.Image.pixelArea()).divide(1000*1000);

var T3_stats_km = T3_area_km.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
});

var T3_stats_km_value = T3_stats_km.get('active');
print(T3_stats_km_value,'T3 final area (km^2)');

var T3_final = T3_noise_removal_closed.reproject('EPSG:4326', null, scale);
Map.addLayer(T3_final, {}, 'T3 (2008 - 2010) active channel binary mask');




// T4

// Filter date range, roi and apply simple cloud processing
var imgCol = merged.filterDate(sDate_T4, eDate_T4) 
                 .filterBounds(roi)
                 .map(function(image) { 
                    var cloudShadowBitMask = 1 << 3; 
                    var cloudsBitMask = 1 << 5; 
                    var qa = image.select('BQA'); 
                    var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) 
                                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0)); 
                    return image.updateMask(mask).multiply(0.0001).clip(roi)
                   }); 
                   

// Define and rename quantiles of interest, display false colour
var bnp50 = ['uBlue_p50', 'Blue_p50', 'Green_p50', 'Red_p50', 'Swir1_p50', 'BQA_p50', 'Nir_p50', 'Swir2_p50'];
var p50 = imgCol.reduce(ee.Reducer.percentile([50])).select(bnp50, bns);
Map.addLayer(p50.clip(roi), params_false, 'T4 (2018 - 2020) false colour composite image', true);

// Apply to each percentile
var mndwi_p50 = Mndwi(p50);
var ndvi_p50 = Ndvi(p50);
var evi_p50 = Evi(p50);
var lswi_p50 = Lswi(p50);

// Water classification from Zou 2018:
var water_p50 = (mndwi_p50.gt(ndvi_p50).or(mndwi_p50.gt(evi_p50))).and(evi_p50.lt(0.1));

// Active channel classification:
var activebelt_p50 = (mndwi_p50.gte(mndwi_param)).and(ndvi_p50.lte(ndvi_param));

// Binary image for active channel classication 
var active_p50 = (water_p50).or(activebelt_p50).rename('active');

// Morphological closing, noise removal and area (# pix and km2) calculation 

var T4_orig = active_p50.select(['active']);

var T4_noise_removal = T4_orig
          .updateMask(T4_orig.connectedPixelCount(cleaning_pixels, false).gte(cleaning_pixels))
          .unmask()
          .clip(roi)
          .reproject('EPSG:4326', null, scale);

var T4_noise_removal_active_only = T4_noise_removal.updateMask(T4_noise_removal.gt(0)).reproject('EPSG:4326', null, scale);

var T4_noise_removal_closed = T4_noise_removal
          .focal_max({kernel: kernel, iterations: 1})
          .focal_min({kernel: kernel, iterations: 1})
          .reproject('EPSG:4326', null, scale)
          .clip(roi);

var T4_noise_removal_closed_active_only = T4_noise_removal_closed.updateMask(T4_noise_removal_closed.gt(0)).reproject('EPSG:4326', null, scale);

var T4_noise_removal_closed_active_stats = T4_noise_removal_closed_active_only.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: roi,
  scale: 30,
}).get('active');

print(T4_noise_removal_closed_active_stats,'T4 final area (# of pix)');

var T4_area_km = T4_noise_removal_closed_active_only.multiply(ee.Image.pixelArea()).divide(1000*1000);

var T4_stats_km = T4_area_km.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
});

var T4_stats_km_value = T4_stats_km.get('active');
print(T4_stats_km_value,'T4 final area (km^2)');

var T4_final = T4_noise_removal_closed.reproject('EPSG:4326', null, scale);
Map.addLayer(T4_final, {}, 'T4 (2018 - 2020) active channel binary mask');



// Sets up RWC analysis - from RivWidthCloud

var AssignDefault = function(x, dv) {
  return(typeof x !== 'undefined' ? x : dv);
};

// here is the modified rwc function, see the example below for usage
exports.rwGenSR_waterMask = function(MAXDISTANCE, FILL_SIZE, MAXDISTANCE_BRANCH_REMOVAL, AOI) {
   
  // assign default values
  // WATER_METHOD = AssignDefault(WATER_METHOD, 'Jones2019');
  MAXDISTANCE = AssignDefault(MAXDISTANCE, 50000);
  FILL_SIZE = AssignDefault(FILL_SIZE, 10);
  MAXDISTANCE_BRANCH_REMOVAL = AssignDefault(MAXDISTANCE_BRANCH_REMOVAL, 250);
  AOI = AssignDefault(AOI, null);
  
  
  var grwl = ee.FeatureCollection("users/eeProject/grwl");
  var lsFun = require('users/eeProject/RivWidthCloudPaper:functions_Landsat578/functions_landsat.js');
  var riverFun = require('users/eeProject/RivWidthCloudPaper:functions_river.js');
  var clWidthFun = require('users/eeProject/RivWidthCloudPaper:functions_centerline_width.js');
  
  exports.CalculateWidth_waterMask = function(imgIn) {
    var crs = imgIn.get('crs');
    var scale = imgIn.get('scale');
    var imgId = imgIn.get('image_id');
    var bound = imgIn.select('riverMask').geometry();
    var angle = imgIn.select('orthDegree');
    var dem = ee.Image("users/eeProject/MERIT");
    var infoEnds = imgIn.select('riverMask');
    var infoExport = imgIn.select('channelMask')
    // .addBands(imgIn.select('^flag.*'))
    .addBands(dem.rename('flag_elevation'));
    var dm = imgIn.select('distanceMap');
    
    var widths = clWidthFun.GetWidth(angle, infoExport, infoEnds, dm, crs, bound, scale, imgId)
    .map(clWidthFun.prepExport);
    return(widths);
  };
  
  // generate function based on user choice
  var tempFUN = function(img) {
    
    // AOI = ee.Algorithms.If(AOI, AOI, img.geometry());
    // img = img.clip(AOI);
    
    // derive water mask and masks for flags
    // var imgOut = lsFun.CalculateWaterAddFlagsSR(image, WATER_METHOD);
    // calculate river mask
    var imgOut = riverFun.ExtractRiver(ee.Image(img), grwl, MAXDISTANCE, FILL_SIZE);
    // calculate centerline
    imgOut = clWidthFun.CalculateCenterline(imgOut);
    // calculate orthogonal direction of the centerline
    imgOut = clWidthFun.CalculateOrthAngle(imgOut);
    // export widths
    var widthOut = exports.CalculateWidth_waterMask(imgOut);
    Map.addLayer(imgOut.selfMask(), {bands: ['cleanedCL'], min: 0, max: 1, palette: ['blue']}, "cleanedCL", true); // RJB added
    
    return(widthOut);
  };
  
  return(tempFUN);
};

// Applies RWC to T1

// calculation starts here
// generate the watermask
// this water mask can be replaced by any user-defined water mask, as long as
// it has three essential properties ['crs', 'scale', 'image_id']
var T1_mask = T1_final
            .reproject('EPSG:32651', null, scale)
            .clip(roi);
// Append crs, scale and image_id for RivWidthCloud:
var T1_mask_for_rwc = T1_mask
.select(['active'], ['waterMask'])
.setMulti({
  crs: T1_mask.projection().crs(),
  scale: 30,
  image_id: '001'
});  
T1_mask_for_rwc = ee.Image(T1_mask_for_rwc);

// Map.addLayer(T1_mask_for_rwc, {}, 'T1 mask for rwc', true);

// initiate and apply the rivwidthcloud function
var T1_rwc = exports.rwGenSR_waterMask(50000, 10, 250);
var T1_widths = T1_rwc(T1_mask_for_rwc);

// examine the first record
//print(T1_widths.limit(100).aggregate_array('width'), 'T1 example widths');
//print(T1_widths.aggregate_stats('width'), 'T1 width stats');
print(T1_widths.aggregate_count('width'), 'T1 count');
print(T1_widths.aggregate_min('width'), 'T1 min');
print(T1_widths.aggregate_max('width'), 'T1 max');
print(T1_widths.aggregate_mean('width'), 'T1 mean');
print(T1_widths.aggregate_total_sd('width'), 'T1 std dev');

// // remove the geometry before exporting the width as CSV file
T1_widths = T1_widths.map(function(f) {return(f.setGeometry(null))});

// export the result as a CSV file into Google drive
Export.table.toDrive({
  collection: T1_widths,
  description: T1_widths_name,
  folder: folder,
  fileNamePrefix: T1_widths_name,
  fileFormat: "CSV"});

// Export Binary Image to Google Drive:
Export.image.toDrive({  
  image: T1_mask, 
  description: T1_binary_name,
  region: roi,
  scale: 30,
  fileFormat: 'GeoTIFF',  
  folder: folder,
  fileNamePrefix: T1_binary_name,
  maxPixels: 1e12,  
}); 

// Applies RWC to T2

// calculation starts here
// generate the watermask
// this water mask can be replaced by any user-defined water mask, as long as
// it has three essential properties ['crs', 'scale', 'image_id']
var T2_mask = T2_final
            .reproject('EPSG:32651', null, scale)
            .clip(roi);
// Append crs, scale and image_id for RivWidthCloud:
var T2_mask_for_rwc = T2_mask
.select(['active'], ['waterMask'])
.setMulti({
  crs: T2_mask.projection().crs(),
  scale: 30,
  image_id: '001'
});  
T2_mask_for_rwc = ee.Image(T2_mask_for_rwc);

//Map.addLayer(T2_mask_for_rwc, {}, 'T2 mask for rwc', true);

// initiate and apply the rivwidthcloud function
var T2_rwc = exports.rwGenSR_waterMask(50000, 10, 250);
var T2_widths = T2_rwc(T2_mask_for_rwc);

// examine the first record
//print(T2_widths.limit(100).aggregate_array('width'), 'T2 example widths');
//print(T2_widths.aggregate_stats('width'), 'T2 width stats');
print(T2_widths.aggregate_count('width'), 'T2 count');
print(T2_widths.aggregate_min('width'), 'T2 min');
print(T2_widths.aggregate_max('width'), 'T2 max');
print(T2_widths.aggregate_mean('width'), 'T2 mean');
print(T2_widths.aggregate_total_sd('width'), 'T2 std dev');

// // remove the geometry before exporting the width as CSV file
T2_widths = T2_widths.map(function(f) {return(f.setGeometry(null))});

// export the result as a CSV file into Google drive
Export.table.toDrive({
  collection: T2_widths,
  description: T2_widths_name,
  folder: folder,
  fileNamePrefix: T2_widths_name,
  fileFormat: "CSV"});

// Export Binary Image to Google Drive:
Export.image.toDrive({  
  image: T2_mask, 
  description: T2_binary_name,
  region: roi,
  scale: 30,
  fileFormat: 'GeoTIFF',  
  folder: folder,
  fileNamePrefix: T2_binary_name,
  maxPixels: 1e12,  
}); 


// Applies RWC to T3

// calculation starts here
// generate the watermask
// this water mask can be replaced by any user-defined water mask, as long as
// it has three essential properties ['crs', 'scale', 'image_id']
var T3_mask = T3_final
            .reproject('EPSG:32651', null, scale)
            .clip(roi);
// Append crs, scale and image_id for RivWidthCloud:
var T3_mask_for_rwc = T3_mask
.select(['active'], ['waterMask'])
.setMulti({
  crs: T3_mask.projection().crs(),
  scale: 30,
  image_id: '001'
});  
T3_mask_for_rwc = ee.Image(T3_mask_for_rwc);

//Map.addLayer(T3_mask_for_rwc, {}, 'T3 mask for rwc', true);

// initiate and apply the rivwidthcloud function
var T3_rwc = exports.rwGenSR_waterMask(50000, 10, 250);
var T3_widths = T3_rwc(T3_mask_for_rwc);

// examine the first record
//print(T3_widths.limit(100).aggregate_array('width'), 'T3 example widths');
//print(T3_widths.aggregate_stats('width'), 'T3 width stats');
print(T3_widths.aggregate_count('width'), 'T3 count');
print(T3_widths.aggregate_min('width'), 'T3 min');
print(T3_widths.aggregate_max('width'), 'T3 max');
print(T3_widths.aggregate_mean('width'), 'T3 mean');
print(T3_widths.aggregate_total_sd('width'), 'T3 std dev');

// // remove the geometry before exporting the width as CSV file
T3_widths = T3_widths.map(function(f) {return(f.setGeometry(null))});

// export the result as a CSV file into Google drive
Export.table.toDrive({
  collection: T3_widths,
  description: T3_widths_name,
  folder: folder,
  fileNamePrefix: T3_widths_name,
  fileFormat: "CSV"});

// Export Binary Image to Google Drive:
Export.image.toDrive({  
  image: T3_mask, 
  description: T3_binary_name,
  region: roi,
  scale: 30,
  fileFormat: 'GeoTIFF',  
  folder: folder,
  fileNamePrefix: T3_binary_name,
  maxPixels: 1e12,  
}); 


// Applies RWC to T4

// calculation starts here
// generate the watermask
// this water mask can be replaced by any user-defined water mask, as long as
// it has three essential properties ['crs', 'scale', 'image_id']
var T4_mask = T4_final
            .reproject('EPSG:32651', null, scale)
            .clip(roi);
// Append crs, scale and image_id for RivWidthCloud:
var T4_mask_for_rwc = T4_mask
.select(['active'], ['waterMask'])
.setMulti({
  crs: T4_mask.projection().crs(),
  scale: 30,
  image_id: '001'
});  
T4_mask_for_rwc = ee.Image(T4_mask_for_rwc);

//Map.addLayer(T4_mask_for_rwc, {}, 'T4 mask for rwc', true);

// initiate and apply the rivwidthcloud function
var T4_rwc = exports.rwGenSR_waterMask(50000, 10, 250);
var T4_widths = T4_rwc(T4_mask_for_rwc);

// examine the first record
//print(T4_widths.limit(100).aggregate_array('width'), 'T1 example widths');
//print(T4_widths.aggregate_stats('width'), 'T4 width stats');
print(T4_widths.aggregate_count('width'), 'T4 count');
print(T4_widths.aggregate_min('width'), 'T4 min');
print(T4_widths.aggregate_max('width'), 'T4 max');
print(T4_widths.aggregate_mean('width'), 'T4 mean');
print(T4_widths.aggregate_total_sd('width'), 'T4 std dev');

// // remove the geometry before exporting the width as CSV file
T4_widths = T4_widths.map(function(f) {return(f.setGeometry(null))});

// export the result as a CSV file into Google drive
Export.table.toDrive({
  collection: T4_widths,
  description: T4_widths_name,
  folder: folder,
  fileNamePrefix: T4_widths_name,
  fileFormat: "CSV"});

// Export Binary Image to Google Drive:
Export.image.toDrive({  
  image: T4_mask, 
  description: T4_binary_name,
  region: roi,
  scale: 30,
  fileFormat: 'GeoTIFF',  
  folder: folder,
  fileNamePrefix: T4_binary_name,
  maxPixels: 1e12,  
}); 



// Calculates similarity coefficients for the  active channel mask for each time period

// T1 - T2

// a (intersection)

var T1_T2_a = T1_final.and(T2_final).reproject('EPSG:4326', null, scale);
var T1_T2_a_active_only = T1_T2_a.updateMask(T1_T2_a.gt(0)).reproject('EPSG:4326', null, scale).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T1_T2_a_active_stats = T1_T2_a_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');


// union 
var T1_T2_union = T1_final.add(T2_final).reproject('EPSG:4326', null, scale);
var T1_T2_union = T1_T2_union.gt(0).reproject('EPSG:4326', null, scale).clip(roi);
var T1_T2_union_active_only = T1_T2_union.updateMask(T1_T2_union.gt(0)).reproject('EPSG:4326', null, scale).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T1_T2_union_active_stats = T1_T2_union_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

// b_c (mismatch)
var T1_T2_b_c = T1_final.subtract(T2_final).reproject('EPSG:4326', null, scale);
var T1_T2_b = T1_T2_b_c.lt(0).reproject('EPSG:4326', null, scale);
var T1_T2_c = T1_T2_b_c.gt(0).reproject('EPSG:4326', null, scale);
var T1_T2_b_active_only = T1_T2_b.updateMask(T1_T2_b.gt(0)).reproject('EPSG:4326', null, scale).clip(roi).multiply(ee.Image.pixelArea()).divide(1000*1000);
var T1_T2_c_active_only = T1_T2_c.updateMask(T1_T2_c.gt(0)).reproject('EPSG:4326', null, scale).clip(roi).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T1_T2_b_active_stats = T1_T2_b_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

var T1_T2_c_active_stats = T1_T2_c_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

var T1_T2_a_value = ee.Number(T1_T2_a_active_stats);
var T1_T2_b_value = ee.Number(T1_T2_b_active_stats);
var T1_T2_c_value = ee.Number(T1_T2_c_active_stats);

var T1_T2_jaccard_index = T1_T2_a_value.divide((T1_T2_a_value.add(T1_T2_b_value).add(T1_T2_c_value)));
print(T1_T2_jaccard_index,'T1 - T2 Jaccard index');

var T1_T2_kulczynski_similarity_coefficient = (T1_T2_a_value.divide(T1_T2_a_value.add(T1_T2_b_value))).add((T1_T2_a_value.divide(T1_T2_a_value.add(T1_T2_c_value)))).divide(2)
print(T1_T2_kulczynski_similarity_coefficient,'T1 - T2 Kulczynski similarity coefficient');

var T1_T2_dice_similarity_coefficient = ((T1_T2_a_value.multiply(2)).divide(((T1_T2_a_value.multiply(2)).add(T1_T2_b_value).add(T1_T2_c_value))));
print(T1_T2_dice_similarity_coefficient,'T1 - T2 Dice similarity coefficient');

var T1_T2_ochiai_similarity_coefficient = T1_T2_a_value.divide(((T1_T2_a_value.add(T1_T2_b_value)).multiply((T1_T2_a_value.add(T1_T2_c_value)))).sqrt());
print(T1_T2_ochiai_similarity_coefficient,'T1 - T2 Ochiai similarity coefficient');



// T2 - T3

// a (intersection)

var T2_T3_a = T2_final.and(T3_final).reproject('EPSG:4326', null, scale);
var T2_T3_a_active_only = T2_T3_a.updateMask(T2_T3_a.gt(0)).reproject('EPSG:4326', null, scale).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T2_T3_a_active_stats = T2_T3_a_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');


// union 
var T2_T3_union = T2_final.add(T3_final).reproject('EPSG:4326', null, scale);
var T2_T3_union = T2_T3_union.gt(0).reproject('EPSG:4326', null, scale).clip(roi);
var T2_T3_union_active_only = T2_T3_union.updateMask(T2_T3_union.gt(0)).reproject('EPSG:4326', null, scale).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T2_T3_union_active_stats = T2_T3_union_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

// b_c (mismatch)
var T2_T3_b_c = T2_final.subtract(T3_final).reproject('EPSG:4326', null, scale);
var T2_T3_b = T2_T3_b_c.lt(0).reproject('EPSG:4326', null, scale);
var T2_T3_c = T2_T3_b_c.gt(0).reproject('EPSG:4326', null, scale);
var T2_T3_b_active_only = T2_T3_b.updateMask(T2_T3_b.gt(0)).reproject('EPSG:4326', null, scale).clip(roi).multiply(ee.Image.pixelArea()).divide(1000*1000);
var T2_T3_c_active_only = T2_T3_c.updateMask(T2_T3_c.gt(0)).reproject('EPSG:4326', null, scale).clip(roi).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T2_T3_b_active_stats = T2_T3_b_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

var T2_T3_c_active_stats = T2_T3_c_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

var T2_T3_a_value = ee.Number(T2_T3_a_active_stats);
var T2_T3_b_value = ee.Number(T2_T3_b_active_stats);
var T2_T3_c_value = ee.Number(T2_T3_c_active_stats);

var T2_T3_jaccard_index = T2_T3_a_value.divide((T2_T3_a_value.add(T2_T3_b_value).add(T2_T3_c_value)));
print(T2_T3_jaccard_index,'T2 - T3 Jaccard index');

var T2_T3_kulczynski_similarity_coefficient = (T2_T3_a_value.divide(T2_T3_a_value.add(T2_T3_b_value))).add((T2_T3_a_value.divide(T2_T3_a_value.add(T2_T3_c_value)))).divide(2)
print(T2_T3_kulczynski_similarity_coefficient,'T2 - T3 Kulczynski similarity coefficient');

var T2_T3_dice_similarity_coefficient = ((T2_T3_a_value.multiply(2)).divide(((T2_T3_a_value.multiply(2)).add(T2_T3_b_value).add(T2_T3_c_value))));
print(T2_T3_dice_similarity_coefficient,'T2 - T3 Dice similarity coefficient');

var T2_T3_ochiai_similarity_coefficient = T2_T3_a_value.divide(((T2_T3_a_value.add(T2_T3_b_value)).multiply((T2_T3_a_value.add(T2_T3_c_value)))).sqrt());
print(T2_T3_ochiai_similarity_coefficient,'T2 - T3 Ochiai similarity coefficient');



// T3 - T4

// a (intersection)

var T3_T4_a = T3_final.and(T4_final).reproject('EPSG:4326', null, scale);
var T3_T4_a_active_only = T3_T4_a.updateMask(T3_T4_a.gt(0)).reproject('EPSG:4326', null, scale).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T3_T4_a_active_stats = T3_T4_a_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');


// union 
var T3_T4_union = T3_final.add(T4_final).reproject('EPSG:4326', null, scale);
var T3_T4_union = T3_T4_union.gt(0).reproject('EPSG:4326', null, scale).clip(roi);
var T3_T4_union_active_only = T3_T4_union.updateMask(T3_T4_union.gt(0)).reproject('EPSG:4326', null, scale).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T3_T4_union_active_stats = T3_T4_union_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

// b_c (mismatch)
var T3_T4_b_c = T3_final.subtract(T4_final).reproject('EPSG:4326', null, scale);
var T3_T4_b = T3_T4_b_c.lt(0).reproject('EPSG:4326', null, scale);
var T3_T4_c = T3_T4_b_c.gt(0).reproject('EPSG:4326', null, scale);
var T3_T4_b_active_only = T3_T4_b.updateMask(T3_T4_b.gt(0)).reproject('EPSG:4326', null, scale).clip(roi).multiply(ee.Image.pixelArea()).divide(1000*1000);
var T3_T4_c_active_only = T3_T4_c.updateMask(T3_T4_c.gt(0)).reproject('EPSG:4326', null, scale).clip(roi).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T3_T4_b_active_stats = T3_T4_b_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

var T3_T4_c_active_stats = T3_T4_c_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

var T3_T4_a_value = ee.Number(T3_T4_a_active_stats);
var T3_T4_b_value = ee.Number(T3_T4_b_active_stats);
var T3_T4_c_value = ee.Number(T3_T4_c_active_stats);

var T3_T4_jaccard_index = T3_T4_a_value.divide((T3_T4_a_value.add(T3_T4_b_value).add(T3_T4_c_value)));
print(T3_T4_jaccard_index,'T3 - T4 Jaccard index');

var T3_T4_kulczynski_similarity_coefficient = (T3_T4_a_value.divide(T3_T4_a_value.add(T3_T4_b_value))).add((T3_T4_a_value.divide(T3_T4_a_value.add(T3_T4_c_value)))).divide(2)
print(T3_T4_kulczynski_similarity_coefficient,'T3 - T4 Kulczynski similarity coefficient');

var T3_T4_dice_similarity_coefficient = ((T3_T4_a_value.multiply(2)).divide(((T3_T4_a_value.multiply(2)).add(T3_T4_b_value).add(T3_T4_c_value))));
print(T3_T4_dice_similarity_coefficient,'T3 - T4 Dice similarity coefficient');

var T3_T4_ochiai_similarity_coefficient = T3_T4_a_value.divide(((T3_T4_a_value.add(T3_T4_b_value)).multiply((T3_T4_a_value.add(T3_T4_c_value)))).sqrt());
print(T3_T4_ochiai_similarity_coefficient,'T3 - T4 Ochiai similarity coefficient');



// T1 - T4

// a (intersection)

var T1_T4_a = T1_final.and(T4_final).reproject('EPSG:4326', null, scale);
var T1_T4_a_active_only = T1_T4_a.updateMask(T1_T4_a.gt(0)).reproject('EPSG:4326', null, scale).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T1_T4_a_active_stats = T1_T4_a_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');


// union 
var T1_T4_union = T1_final.add(T4_final).reproject('EPSG:4326', null, scale);
var T1_T4_union = T1_T4_union.gt(0).reproject('EPSG:4326', null, scale).clip(roi);
var T1_T4_union_active_only = T1_T4_union.updateMask(T1_T4_union.gt(0)).reproject('EPSG:4326', null, scale).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T1_T4_union_active_stats = T1_T4_union_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

// b_c (mismatch)
var T1_T4_b_c = T1_final.subtract(T4_final).reproject('EPSG:4326', null, scale);
var T1_T4_b = T1_T4_b_c.lt(0).reproject('EPSG:4326', null, scale);
var T1_T4_c = T1_T4_b_c.gt(0).reproject('EPSG:4326', null, scale);
var T1_T4_b_active_only = T1_T4_b.updateMask(T1_T4_b.gt(0)).reproject('EPSG:4326', null, scale).clip(roi).multiply(ee.Image.pixelArea()).divide(1000*1000);
var T1_T4_c_active_only = T1_T4_c.updateMask(T1_T4_c.gt(0)).reproject('EPSG:4326', null, scale).clip(roi).multiply(ee.Image.pixelArea()).divide(1000*1000);

var T1_T4_b_active_stats = T1_T4_b_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

var T1_T4_c_active_stats = T1_T4_c_active_only.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e12
}).get('active');

var T1_T4_a_value = ee.Number(T1_T4_a_active_stats);
var T1_T4_b_value = ee.Number(T1_T4_b_active_stats);
var T1_T4_c_value = ee.Number(T1_T4_c_active_stats);

var T1_T4_jaccard_index = T1_T4_a_value.divide((T1_T4_a_value.add(T1_T4_b_value).add(T1_T4_c_value)));
print(T1_T4_jaccard_index,'T1 - T4 Jaccard index');

var T1_T4_kulczynski_similarity_coefficient = (T1_T4_a_value.divide(T1_T4_a_value.add(T1_T4_b_value))).add((T1_T4_a_value.divide(T1_T4_a_value.add(T1_T4_c_value)))).divide(2)
print(T1_T4_kulczynski_similarity_coefficient,'T1 - T4 Kulczynski similarity coefficient');

var T1_T4_dice_similarity_coefficient = ((T1_T4_a_value.multiply(2)).divide(((T1_T4_a_value.multiply(2)).add(T1_T4_b_value).add(T1_T4_c_value))));
print(T1_T4_dice_similarity_coefficient,'T1 - T4 Dice similarity coefficient');

var T1_T4_ochiai_similarity_coefficient = T1_T4_a_value.divide(((T1_T4_a_value.add(T1_T4_b_value)).multiply((T1_T4_a_value.add(T1_T4_c_value)))).sqrt());
print(T1_T4_ochiai_similarity_coefficient,'T1 - T4 Ochiai similarity coefficient');