// Google Earth Engine example code for creating binary masks 
//  Extract the wetted and the active river channel and quantify change from collections of Landsat imagery (30 m).
//  v1.1 - June 2023

//  User defines the region of interest (ROI) - Step 1
//  User defines the time period for analysis - Step 2
//  User defines a set of parameters to classify and clean the active river channel mask - Step 3

//  Outputs are displayed as layers on the map and can be exported to Google Drive

//  Adapted from Boothroyd et al., 2020  (https://doi.org/10.1002/wat2.1496)
// modified to have a smaller area and extract wetted channel as well as other outputs.

//11 jan 24 calculate active channel subtract wetted channel added
// 5 feb 24 calcuations and output for active channel added 

// INPUTS

// Step 1 = Set the region of interest (ROI)
// (you can also draw a polygon in the console and rename it to "roi"

// Karnali River example roi:

var roi = ee.Geometry.Polygon([
  [[81.25174024228566,28.65473982399883],
[81.1714027178716,28.570650600075638],
[81.26341321591848,28.529938112248203],
[81.29499890927785,28.63877106609554]]]);



// Step 2 = define the time windows (TW's) for the analysis:
// (six month period selected to coincide with the post-monsoon period each year)

var TW1_start = "2001-10-01"; 
var TW1_end = "2002-01-01";

var TW2_start = "2016-10-01"; 
var TW2_end = "2017-01-01";


// Step 3 = define parameters for the riverscape classification and cleaning:

var cloudy_score =                50;           // default = 50         Value between 0-100, higher value includes more cloudy images.
var percentile_reducer =          33;           // default = 33         Value between 0-100, values in the range 10-50 usually work well but depends on cloudyness of the site. 
var ndvi_thresh =                 0.17         // default = 0.225      Sets an NDVI threshold, suggest 0.15 to 0.3. 
var val_thresh =                  0.2         // default = 0.125      Sets a threshold for the first brightness indicator, suggest 0.1 - 0.2. 
var red_thresh =                  0.2     // default = 0.125      Sets a threshold for the second brightness indicator, suggest 0.1 - 0.2.
var morph_radius =                15;           // default = 15         Sets radius of morphological operator, suggest 0.5 pixel = 15 m. 
var disconnectedSizeMeters =      6000;         // default = 6000       Sets a threshold to remove disconnected/misclassified pixels not connected to the main channel, higher value removes larger disconnected objects. 
var islandSizeMeters =            100;          // default = 100        Sets a threshold to fill small holes/islands within the active channel mask.
var exportResolutionMeters =      30;           // default = 30         Sets the pixel size.

// Step 4 = set the filenames/folders to save the outputs:

var folder = ('Karnali_Landsat_Lowflow_Masks');
var river = ('Karnali');



// *** Code below this point shouldn't need modifying ***

// This sets the view:
Map.addLayer(roi, {}, 'ROI', false);
Map.centerObject(roi,13);

// This defines the image collection to use for the analysis (Landsat Tier 1 Level 2)
// Gets Landsat 5,7,8,9 atmospherically corrected surface reflectance datasets
var ls5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
var ls7 = (ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
  .filterDate('1999-04-15', '2013-05-30')) // note LE7 images affected by failure of Scan Line Corrector (SLC) after May03 
var ls8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
var ls9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")

// Function to apply scaling factor (updated 27/01/2023)
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  return image.addBands(opticalBands, null, true)
}

// Standardises band names
var bn9 = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B6', 'QA_PIXEL', 'SR_B5', 'SR_B7'];
var bn8 = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B6', 'QA_PIXEL', 'SR_B5', 'SR_B7'];
var bn7 = ['SR_B1', 'SR_B1', 'SR_B2', 'SR_B3', 'SR_B5', 'QA_PIXEL', 'SR_B4', 'SR_B7'];
var bn5 = ['SR_B1', 'SR_B1', 'SR_B2', 'SR_B3', 'SR_B5', 'QA_PIXEL', 'SR_B4', 'SR_B7'];
var bns = ['uBlue', 'Blue', 'Green', 'Red', 'Swir1', 'QA_PIXEL', 'Nir', 'Swir2'];

var ls5_scaled = ls5.map(applyScaleFactors).select(bn5, bns);
var ls7_scaled = ls7.map(applyScaleFactors).select(bn7, bns);
var ls8_scaled = ls8.map(applyScaleFactors).select(bn8, bns);
var ls9_scaled = ls9.map(applyScaleFactors).select(bn9, bns);


// Merges Landsat 5,7,8,9 datasets and sets date filter for analysis 
var Lcol = ls5_scaled.merge(ls7_scaled).merge(ls8_scaled).merge(ls9_scaled).filterBounds(roi);

// Defines the cloud masking procedure 
var maskcloud = function(image) {
                    var cloudShadowBitMask = 1 << 3; 
                    var cloudsBitMask = 1 << 5; 
                    var qa = image.select('QA_PIXEL'); 
                    var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) 
                                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0)); 
                    return image.updateMask(mask).clip(roi).resample('bicubic')
                    .set('system:time_start', image.get('system:time_start'));
                    }; 

// Defines the multispectral indices used to classify the active channel
var Mndwi = function(image) {
  var mndwi = image.normalizedDifference(['Green', 'Swir1'])
  return(mndwi.rename(['mndwi']));
};

var Ndvi = function(image) {
  var ndvi = image.normalizedDifference(['Nir', 'Red'])
  return(ndvi.rename(['ndvi']));
};

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


// Defines the function to remove disconnected pixels and fill small islands:

function removeDisconnected(activeImage) {
  // reduceConnectedComponents expects an interger image
  activeImage = activeImage.int();
   
  // Define neighborhood based on user parameters
  var connectedPixels1 = ee.Number(disconnectedSizeMeters)
    .divide(exportResolutionMeters).int();
    
  var connectedPixels2 = ee.Number(islandSizeMeters)
    .divide(exportResolutionMeters).int();
 
  // Remove disconnected pixels
  var landFilled = activeImage.addBands(activeImage)
   .reduceConnectedComponents(
     ee.Reducer.median(), 'mndwi', connectedPixels1)
   .unmask(99).eq(99).and(activeImage.neq(0));

  // Remove small islands  
  var activeFilled = landFilled.addBands(landFilled)
    .reduceConnectedComponents(
      ee.Reducer.median(), 'mndwi_1', connectedPixels2)
    .unmask(99).eq(99).and(landFilled.neq(1));    
   
  // Output
  return activeFilled;
}


// Sets the visualisation parameters
var params_true = {crs: 'EPSG:4326', region: roi, min: 0.0, max: 0.3, bands: ["Red", "Green", "Blue"]};
var activechannelViz_TW1 = {min: 0, max: 1, palette: ['white','b550ab']};
var activechannelViz_TW2 = {min: 0, max: 1, palette: ['white','71b850']};
var waterViz_TW1 = {min: 0, max: 1, palette: ['white','00FFFF']};
var waterViz_TW2 = {min: 0, max: 1, palette: ['white','0000FF']};
var channelViz = {min:0, max: 1, palette: ['white', 'dbdbdb']}
var erosionViz = {min: 0, max: 1, palette: ['white','e03d34']};
var accretionViz = {min: 0, max: 1, palette: ['white','346ae0']};




// Applies the GEE code to the time window one (TW1).

// Filters the image collection and applies cloud masking procedure
var imgCol_TW1 = Lcol.filterDate(TW1_start, TW1_end) 
                 .filter(ee.Filter.lt('CLOUD_COVER', cloudy_score))
                 .map(maskcloud);
                 
//Map.addLayer(imgCol_TW1, params_true, 'Image Collection ' + TW1_year , true); //trying to see what it looks like before cloud masking

// Applies a percentile reducer to composite the cloud-masked images so that the RGB image is displayed
var percentileComposite_TW1 = imgCol_TW1.select(bns)
var percentile_rename_TW1 = ('_p' + percentile_reducer)
var percentileComposite_TW1 = percentileComposite_TW1.reduce(ee.Reducer.percentile([percentile_reducer])).regexpRename(percentile_rename_TW1, '');
var percentileComposite_TW1_true = percentileComposite_TW1.select(['Red','Green','Blue']);
var percentileComposite_TW1_red = percentileComposite_TW1.select(['Red']);

var hsv_rgb_TW1 = percentileComposite_TW1.select(['Red','Green','Blue']).rgbToHsv();
var mndwi_pC_TW1 = Mndwi(percentileComposite_TW1);
var ndvi_pC_TW1 = Ndvi(percentileComposite_TW1);
var evi_pC_TW1 = Evi(percentileComposite_TW1);

var water_TW1 = (mndwi_pC_TW1.gt(ndvi_pC_TW1).or(mndwi_pC_TW1.gt(evi_pC_TW1))).and(evi_pC_TW1.lt(0.1)); 
var non_veg_TW1 = ndvi_pC_TW1.lte(ndvi_thresh)
var bright_TW1 = hsv_rgb_TW1.select('value').gte(val_thresh).and(percentileComposite_TW1_red.gte(red_thresh));
var active_channel_TW1 = water_TW1.or(non_veg_TW1).or(bright_TW1).rename('mndwi')

// Applies standard image processing procedures to clean/remove noise from the active channel mask
var active_channel_original_TW1 = active_channel_TW1.select(['mndwi'])

var active_channel_opened_TW1 = active_channel_original_TW1
    .focalMax({
    'radius': morph_radius,
    'units': 'meters',
    'kernelType': 'circle'})
    .focalMin({
    'radius': morph_radius,
    'units': 'meters',
    'kernelType': 'circle'});

var active_channel_output_TW1 = removeDisconnected(active_channel_opened_TW1).not();


// Display TW1 outputs in console and on the map:

var TW1_year = (TW1_start.slice(-10, 4));
print(imgCol_TW1.size(),('Number of Landsat images used ' + TW1_year));
print(imgCol_TW1);
Map.addLayer(percentileComposite_TW1.clip(roi), params_true, 'RGB image ' + TW1_year , true);
Map.addLayer(active_channel_output_TW1.clip(roi).selfMask(), activechannelViz_TW1, 'Active channel ' + TW1_year , false);
//add wetted channel
Map.addLayer(water_TW1.clip(roi).selfMask(), waterViz_TW1, 'Wetted channel ' + TW1_year , false);



// Applies the GEE code to the time window two (TW2).

// Filters the image collection and applies cloud masking procedure
var imgCol_TW2 = Lcol.filterDate(TW2_start, TW2_end) 
                 .filter(ee.Filter.lt('CLOUD_COVER', cloudy_score))
                 .map(maskcloud);

// Applies a percentile reducer to composite the cloud-masked images so that the RGB image is displayed
var percentileComposite_TW2 = imgCol_TW2.select(bns)
var percentile_rename_TW2 = ('_p' + percentile_reducer)
var percentileComposite_TW2 = percentileComposite_TW2.reduce(ee.Reducer.percentile([percentile_reducer])).regexpRename(percentile_rename_TW2, '');
var percentileComposite_TW2_true = percentileComposite_TW2.select(['Red','Green','Blue']);
var percentileComposite_TW2_red = percentileComposite_TW2.select(['Red']);

var hsv_rgb_TW2 = percentileComposite_TW2.select(['Red','Green','Blue']).rgbToHsv();
var mndwi_pC_TW2 = Mndwi(percentileComposite_TW2);
var ndvi_pC_TW2 = Ndvi(percentileComposite_TW2);
var evi_pC_TW2 = Evi(percentileComposite_TW2);

var water_TW2 = (mndwi_pC_TW2.gt(ndvi_pC_TW2).or(mndwi_pC_TW2.gt(evi_pC_TW2))).and(evi_pC_TW2.lt(0.1)); 
var non_veg_TW2 = ndvi_pC_TW2.lte(ndvi_thresh)
var bright_TW2 = hsv_rgb_TW2.select('value').gte(val_thresh).and(percentileComposite_TW2_red.gte(red_thresh));
var active_channel_TW2 = water_TW2.or(non_veg_TW2).or(bright_TW2).rename('mndwi')

// Applies standard image processing procedures to clean/remove noise from the active channel mask
var active_channel_original_TW2 = active_channel_TW2.select(['mndwi'])

var active_channel_opened_TW2 = active_channel_original_TW2
    .focalMax({
    'radius': morph_radius,
    'units': 'meters',
    'kernelType': 'circle'})
    .focalMin({
    'radius': morph_radius,
    'units': 'meters',
    'kernelType': 'circle'});

var active_channel_output_TW2 = removeDisconnected(active_channel_opened_TW2).not();


// Display TW2 outputs in console and on the map:

var TW2_year = (TW2_start.slice(-10, 4));
print(imgCol_TW2.size(),('Number of Landsat images used ' + TW2_year));
print(imgCol_TW2);
Map.addLayer(percentileComposite_TW2.clip(roi), params_true, 'RGB image ' + TW2_year , false);
Map.addLayer(active_channel_output_TW2.clip(roi).selfMask(), activechannelViz_TW2, 'Active channel ' + TW2_year , false);
//add wetted channel
Map.addLayer(water_TW2.clip(roi).selfMask(), waterViz_TW2, 'Wetted channel ' + TW2_year , false);


// Calculate difference between TW1 and TW2 (not used)
// 11 jan 24 calculate active channel subtract wetted channel

var Active_TW1 = active_channel_output_TW1.unmask();
var Active_TW2 = active_channel_output_TW2.unmask();

var TW1_TW2_non_channel_no_change = Active_TW1.eq(0).and(Active_TW2.eq(0));
var TW1_TW2_channel_no_change = Active_TW1.eq(1).and(Active_TW2.eq(1));
var TW1_TW2_erosion = Active_TW1.eq(0).and(Active_TW2.eq(1));
var TW1_TW2_accretion = Active_TW1.eq(1).and(Active_TW2.eq(0));
var TW1_active_only = Active_TW1.eq(1).subtract(water_TW1.eq(1)).clip(roi).selfMask();
var TW2_active_only = Active_TW2.eq(1).subtract(water_TW2.eq(1)).clip(roi).selfMask();

//Map.addLayer(TW1_TW2_channel_no_change.clip(roi).selfMask(), channelViz, 'No change', true);
//Map.addLayer(TW1_TW2_erosion.clip(roi).selfMask(), erosionViz, 'Erosion', true);
//Map.addLayer(TW1_TW2_accretion.clip(roi).selfMask(), accretionViz, 'Accretion', true);
Map.addLayer(TW1_active_only.clip(roi).selfMask(), accretionViz, 'Active_only' + TW1_year, true);
Map.addLayer(TW2_active_only.clip(roi).selfMask(), erosionViz, 'Active_only' + TW2_year, true);


//print(ui.Thumbnail(TW1_active_only, erosionViz),"Active_only" + TW1_year); 

//print(ui.Thumbnail(TW2_active_only, erosionViz),"Active_only" + TW2_year); 

// Calculate # and area of pixels of each class

var areaImage_TW1_TW2_no_change = TW1_TW2_channel_no_change.selfMask().multiply(ee.Image.pixelArea())

var area_TW1_TW2_no_change = areaImage_TW1_TW2_no_change.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 10,
  maxPixels: 1e13
  })
  
var area_TW1_TW2_no_change = ee.Number(area_TW1_TW2_no_change.get('mndwi_1_1')).divide(1e6);


var areaImage_TW1_TW2_accretion = TW1_TW2_accretion.selfMask().multiply(ee.Image.pixelArea())

var area_TW1_TW2_accretion = areaImage_TW1_TW2_accretion.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 10,
  maxPixels: 1e13
  })

var area_TW1_TW2_accretion = ee.Number(area_TW1_TW2_accretion.get('mndwi_1_1')).divide(1e6);


var areaImage_TW1_TW2_erosion = TW1_TW2_erosion.selfMask().multiply(ee.Image.pixelArea())

var area_TW1_TW2_erosion = areaImage_TW1_TW2_erosion.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 10,
  maxPixels: 1e13
  })
  
  
var area_TW1_TW2_erosion = ee.Number(area_TW1_TW2_erosion.get('mndwi_1_1')).divide(1e6);

var areaImage_Active_TW1 = Active_TW1.selfMask().multiply(ee.Image.pixelArea())
var areaImage_Active_TW2 = Active_TW2.selfMask().multiply(ee.Image.pixelArea())

var area_Active_TW1 = areaImage_Active_TW1.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 10,
  maxPixels: 1e13
  })
  
var area_Active_TW2 = areaImage_Active_TW2.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 10,
  maxPixels: 1e13
  })
  
var area_Active_TW1 = ee.Number(area_Active_TW1.get('mndwi_1_1')).divide(1e6);
var area_Active_TW2 = ee.Number(area_Active_TW2.get('mndwi_1_1')).divide(1e6);

//print(area_TW1_TW2_no_change,('Area of no change [km2]'));
//print(area_TW1_TW2_accretion,('Area of accretion [km2]'));
//print(area_TW1_TW2_erosion,('Area of erosion [km2]'));

print(area_Active_TW1,('active channel '+ TW1_year));
print(area_Active_TW2,('active channel '+ TW2_year));

// Export outputs to Google Drive:

// Export.image.toDrive({  
// image: active_channel_output_TW1, 
//   description: (river + '_active_channel_' + TW1_year),
//   fileNamePrefix: (river + '_active_channel_' + TW1_year),
//   region: roi,
  // scale: 30,
  // crs: 'EPSG:4326',
  // fileFormat: 'GeoTIFF',  
  // folder: folder,
  // maxPixels: 1e13
// }); 

// Export.image.toDrive({  
//   image: active_channel_output_TW2, 
//   description: (river + '_active_channel_' + TW2_year),
//   fileNamePrefix: (river + '_active_channel_' + TW2_year),
//   region: roi,   
//   scale: 30,
//   crs: 'EPSG:4326',
//   fileFormat: 'GeoTIFF',  
//   folder: folder,
//   maxPixels: 1e13
// }); 
 
Export.image.toDrive({  
image: Active_TW1, 
  description: (river + '_active_channel' + TW1_year),
  fileNamePrefix: (river + '_active_channel' + TW1_year),
  region: roi,
  scale: 30,
  crs: 'EPSG:4326',
  fileFormat: 'GeoTIFF',  
  folder: folder,
  maxPixels: 1e13
}); 

Export.image.toDrive({  
  image: Active_TW2, 
  description: (river + '_active_channel' + TW2_year),
  fileNamePrefix: (river + '_active_channel' + TW2_year),
  region: roi,   
  scale: 30,
  crs: 'EPSG:4326',
  fileFormat: 'GeoTIFF',  
  folder: folder,
  maxPixels: 1e13
}); 

// Export.image.toDrive({  
//   image: percentileComposite_TW1_true, 
//  description: (river + '_compsite_RGB_' + TW1_year),
//  fileNamePrefix: (river + '_composite_RGB_' + TW1_year),
//   region: roi,
//   scale: 30,
//   crs: 'EPSG:4326',
//   fileFormat: 'GeoTIFF',  
//   folder: folder,
//   maxPixels: 1e13
// }); 

// Export.image.toDrive({  
//   image: percentileComposite_TW2_true, 
//   description: (river + '_composite_RGB_' + TW2_year),
//   fileNamePrefix: (river + '_composite_RGB_' + TW2_year),
//   region: roi,
//   scale: 30,
//   crs: 'EPSG:4326',
//   fileFormat: 'GeoTIFF',  
//   folder: folder,
//   maxPixels: 1e13
// }); 

// Export.image.toDrive({  
//   image: TW1_TW2_channel_no_change, 
//   description: (river + '_active_channel_no_change_' + TW1_year + '_' + TW2_year),
//   fileNamePrefix: (river + '_active_channel_no_change_' + TW1_year + '_' + TW2_year),
//   region: roi,
//   scale: 30,
//   crs: 'EPSG:4326',
//   fileFormat: 'GeoTIFF',  
//   folder: folder,
//   maxPixels: 1e13
// }); 

// Export.image.toDrive({  
//   image: TW1_TW2_accretion, 
//   description: (river + '_active_channel_accretion_' + TW1_year + '_' + TW2_year),
//   fileNamePrefix: (river + '_active_channel_accretion_' + TW1_year + '_' + TW2_year),
//   region: roi,
//   scale: 30,
//   crs: 'EPSG:4326',
//   fileFormat: 'GeoTIFF',  
//   folder: folder,
//   maxPixels: 1e13
// }); 

// Export.image.toDrive({  
//   image: TW1_TW2_erosion, 
//   description: (river + '_active_channel_erosion_' + TW1_year + '_' + TW2_year),
//   fileNamePrefix: (river + '_active_channel_erosion_' + TW1_year + '_' + TW2_year),
//   region: roi,
//   scale: 30,
//   crs: 'EPSG:4326',
//   fileFormat: 'GeoTIFF',  
//   folder: folder,
//   maxPixels: 1e13
// }); 

// Export.image.toDrive({  
//  image: water_TW1, 
//  description: (river + '_wetted_channel_' + TW1_year),
//  fileNamePrefix: (river + '_wetted_channel_' + TW1_year),
//  region: roi,
//  scale: 30,
//  crs: 'EPSG:4326',
//  fileFormat: 'GeoTIFF',  
//  folder: folder,
//  maxPixels: 1e13
//}); 

//Export.image.toDrive({  
//  image: water_TW2, 
//  description: (river + '_wetted_channel_' + TW2_year),
//  fileNamePrefix: (river + '_wetted_channel_' + TW2_year),
//  region: roi,
//  scale: 30,
//  crs: 'EPSG:4326',
//  fileFormat: 'GeoTIFF',  
//  folder: folder,
//  maxPixels: 1e13
// }); 

