////////////////////////////////////////
///////////// PREPARE IMAGE ////////////
////////////////////////////////////////

// "Zoom in" the map and center it on Mexico study area. Map must be "zoomed in" to level 10 or greater to function.
// If map is too "zoomed out", the computation is too intensive and the script does not run.
// After the script has been run and the area estimate has been calculated, you can safely
// zoom in, zoom out, and pan around the map.
Map.centerObject(mexico_mar, 10);

// Establish function for masking clouds from Sentinel-2 image collection
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

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

  return image.updateMask(mask).divide(10000);
}

// Load image collection of Sentinel-2 surface reflectance data
var sentinel_collection = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate('2020-01-01', '2020-05-31')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .map(maskS2clouds);

// Establish variable for visualizing data on screen (true color)
var median_vis = {
  min: 0.0,
  max: 0.2,
  bands: ['B4_median', 'B3_median', 'B2_median'],
};

// Reduce Sentinel image collection to single image
var sentinel_image = sentinel_collection.reduce(ee.Reducer.median());

// Clip Sentinel image to Mexico boundary
var mexico_image = sentinel_image.clip(mexico_mar);

// Display Mexico image
Map.addLayer(mexico_image, median_vis, "Mexico Image");


////////////////////////////////////////
///////////// FILTER IMAGE /////////////
////////////////////////////////////////

// Create "NDWI" layer to mask out water features
var ndwi_layer = mexico_image.normalizedDifference(['B3_median', 'B8_median']);

// Mask water out of Sentinel image using NDWI layer
var mexico_ndwi = mexico_image.updateMask(ndwi_layer.lte(0.0));
//Map.addLayer(mexico_ndwi, median_vis, 'NDWI mask applied');

// Create "elevation" layer to mask out high-elevation features
var srtm = ee.Image('USGS/SRTMGL1_003');
var elevation_layer = srtm.select('elevation');

// Mask high elevations out of Sentinel image using elevation layer
var mexico_elevation = mexico_ndwi.updateMask(elevation_layer.lte(20.0));
//Map.addLayer(mexico_elevation, median_vis, 'Elevation mask applied');

// Create "NDMI" layer to mask out non-mangrove candidates
var ndmi_layer = mexico_image.normalizedDifference(['B12_median', 'B3_median']);

// Mask non-mangrove candidates out of Sentinel image using NDMI layer
var mexico_ndmi = mexico_elevation.updateMask(ndmi_layer.lte(0.1));
//Map.addLayer(mexico_ndmi, median_vis, "NDMI mask applied");


////////////////////////////////////////
/////////// CLASSIFY IMAGE /////////////
////////////////////////////////////////

// Designate which spectral bands to include in the classification
var bands = ['B8_median', 'B4_median', 'B3_median', 'B2_median'];

// Train the classifier using the training points
var training = mexico_image.select(bands).sampleRegions({
  collection: mexico_training,
  properties: ['class'],
  scale: 10
});

// Establish the classifier
var classifier = ee.Classifier.smileRandomForest(200)
    .train({
      features: training,
      classProperty: 'class',
      inputProperties: bands
    });

// Classify the image
var mexico_classified = mexico_ndmi.select(bands).classify(classifier);

// Mask out all non-mangrove pixels from the classified thematic map
var mexico_class_mask = mexico_classified.updateMask(mexico_classified.lte(0));

// Display the preliminary classified thematic map
//Map.addLayer(mexico_class_mask, {palette: ('FF0000')}, 'Preliminary Mexico Mangrove Classification');


////////////////////////////////////////
//////// REFINE CLASSIFICATION /////////
////////////////////////////////////////

// Display Bunting mangrove layer for reference
var mexico_bunting_2016 = mexico_bunting.reduceToImage({
    properties: ['pxlval'],
    reducer: ee.Reducer.first()
});
//Map.addLayer(mexico_bunting_2016, {palette: ('0000FF')}, 'Bunting 2016');

// Display only the mangrove training points for reference
//Map.addLayer(mexico_mang_training, {palette: ('00FF00')}, 'Mexico Mangrove Training Points');

// Make sure that the manual adjustment geometry is of type FeatureCollection
var mexico_include = ee.FeatureCollection(mexico_include);

// Dissolve the polygons in the Feature Collection, so the overlapping polygons
// don't "cancel each other out"
var mexico_include = mexico_include.union();

// Clip the classification to the Feature Collection of manual adjustments, and display
var mexico_mangroves_without_sian_kaan = mexico_class_mask.clip(mexico_include);
//Map.addLayer(mexico_mangroves_without_sian_kaan, {palette: ('00FF00')}, 'Mexico Mangroves without Sian Kaan');


////////////////////////////////////////
////////// SIAN KA'AN ADDITION//////////
////////////////////////////////////////

// Clip Sentinel image to Sian Ka'an boundary
var sian_kaan_image = sentinel_image.clip(sian_kaan);

// Mask water out of Sian Ka'an image using NDWI layer
var sian_kaan_ndwi = sian_kaan_image.updateMask(ndwi_layer.lte(0.0));
//Map.addLayer(sian_kaan_ndwi, median_vis, 'NDWI mask applied');

// Mask high elevations out of Sian Ka'an image using elevation layer
var sian_kaan_elevation = sian_kaan_ndwi.updateMask(elevation_layer.lte(20.0));
//Map.addLayer(sian_kaan_elevation, median_vis, 'Elevation mask applied');

// Mask non-mangrove candidates out of Sian Ka'an image using NDMI layer
var sian_kaan_ndmi = sian_kaan_elevation.updateMask(ndmi_layer.lte(0.4));
//Map.addLayer(sian_kaan_ndmi, median_vis, "NDMI mask applied");

// Designate which spectral bands to include in the classification
var bands = ['B8_median', 'B4_median', 'B3_median', 'B2_median'];

// Train the classifier using the training points
var training = sian_kaan_image.select(bands).sampleRegions({
  collection: mexico_training,
  properties: ['class'],
  scale: 10
});

// Establish the classifier
var classifier = ee.Classifier.smileRandomForest(200)
    .train({
      features: training,
      classProperty: 'class',
      inputProperties: bands
    });

// Classify the image
var sian_kaan_classified = sian_kaan_ndmi.select(bands).classify(classifier);

// Mask out all non-mangrove pixels from the classified thematic map
var sian_kaan_class_mask = sian_kaan_classified.updateMask(sian_kaan_classified.lte(0));

var sian_kaan_include = sian_kaan_include.union();
var sian_kaan_mangroves = sian_kaan_class_mask.clip(sian_kaan_include);

// Display the preliminary classified thematic map
//Map.addLayer(sian_kaan_class_mask, {palette: ('0000FF')}, 'Preliminary Sian Kaan Mangrove Classification');

// Put the "Mexico without Sian Ka'an" and "Sian Ka'an" images together in image collection
// Then reduce the image collection to a single merged "Mexico Total" image
var mexico_collection = ee.ImageCollection([mexico_mangroves_without_sian_kaan, sian_kaan_mangroves]);
var mexico_mangroves = mexico_collection.reduce(ee.Reducer.max());
Map.addLayer(mexico_mangroves, {palette: '00FF00'}, 'Mexico MAR Mangroves');

var total_include = mexico_include.merge(sian_kaan_include);
var total_include = total_include.union();

// Calculate mangrove area for Mexico, and print to console
var mexico_mangrove_area = mexico_mangroves.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: total_include,
    scale: 10,
    maxPixels: 1e9
});
var mexico_mangrove_area = ee.Number(mexico_mangrove_area.get('classification_max')).divide(10000);
var mexico_mangrove_area = ee.Number(mexico_mangrove_area).format();
print('Mexico MAR mangrove area (km2) as of 2020:', mexico_mangrove_area);

////////////////////////////////////////
////////////// EXPORT MAP //////////////
////////////////////////////////////////

// Convert Mexico mangrove raster to polygon
var mexico_mangroves_polygon = mexico_mangroves.toInt().reduceToVectors({
  geometry: total_include,
  crs: mexico_mangroves.projection(),
  scale: 10,
  geometryType: 'polygon',
  eightConnected: false,
  maxPixels: 1e9
});

// Export Mexico mangrove polygon to Google Drive
//Export.table.toDrive({
//  collection: mexico_mangroves_polygon,
//  description:'mexico_mar_mangroves',
//  fileFormat: 'KMZ'
//});


////////////////////////////////////////
////////////// VALIDATION //////////////
////////////////////////////////////////

// Generate random points for "mangrove presence" validation, and export to Google Drive
// I will generate "mangrove absence" buffer and points in GIS program
var mexico_validation_points_presence = ee.FeatureCollection.randomPoints(mexico_mangroves_polygon, 2000);
//Export.table.toDrive({collection: mexico_validation_points_presence, fileFormat: 'GeoJSON'});