// This example demonstrates how to use Cloud Score+ QA bands to generate a 
// clear median composite for a specified date range.

// Harmonized Sentinel-2 Level 2A collection.
var s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED');

// Cloud Score+ image collection. Note Cloud Score+ is produced from Sentinel-2
// Level 1C data and can be applied to either L1C or L2A collections.
var csPlus = ee.ImageCollection('GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED');

//https://code.earthengine.google.com/?asset=projects/ee-danielelacecilia/assets/catchment_04326
// Region of interest. Sele
var ROI = ee.Geometry.Polygon(
	[[[14.890431483581969,40.58248757792104],
[14.934720118835875,40.51387691768456],
[15.001668055847594,40.54466910153016],
[14.945706446960875,40.60360445837052]]], null, false);

var name_export = 'GEE_Sele';
var name_export_Corine = ee.String(name_export).cat('_OPAC_Corine_2501');
var name_export_RGB = ee.String(name_export).cat('_RGB_2501');
var name_export_NDVI = ee.String(name_export).cat('_NDVI_2501');



// Date for mapping with random forest
var start_date = '2025-01-01'
var end_date   = '2025-01-30'



// Use 'cs' or 'cs_cdf', depending on your use-case; see docs for guidance.
var QA_BAND = 'cs';

// The threshold for masking; values between 0.50 and 0.65 generally work well.
// Higher values will remove thin clouds, haze & cirrus shadows.
var CLEAR_THRESHOLD = 0.60;

// Make composite.
var composite = s2
    .filterBounds(ROI)
    .filterDate(start_date, end_date)
    .linkCollection(csPlus, [QA_BAND])
    .map(function(img) {
      return img.updateMask(img.select(QA_BAND).gte(CLEAR_THRESHOLD))
    })
    .median();

// Sentinel-2 visualization parameters.
var s2Viz = {bands: ['B4', 'B3', 'B2'], min: 0, max: 2500};

Map.addLayer(composite, s2Viz, 'median composite',false);
Map.centerObject(ROI, 11);
//Map.addLayer(ROI, {color: 'red'}, false)
//print(composite)










var composite = composite.select(['B1','B2','B3','B4','B5','B6','B7','B8','B11','B12']);

//Add bands
var ndvi = composite.normalizedDifference(['B8','B4']).rename('NDVI');
var ndwi = composite.normalizedDifference(['B8','B11']).rename('NDWI');
var ni = composite.expression(
    '((2 * B8 - B4 - B12)/(2 * B8 + B4 + B12))', {
      'B8': composite.select('B8'),
      'B4': composite.select('B4'),
      'B12': composite.select('B12')
      }).rename('NI');
var newBands = ee.Image([ndvi,ndwi,ni]);
composite = composite.addBands(newBands);

print(composite)






var trainingSample = ee.FeatureCollection("projects/ee-danielelacecilia/assets/cicci");

var bands = ['B1','B2','B3','B4','B5','B6','B7','B8','B11','B12','NDVI','NDWI','NI'];
var label = 'Class';

// Train a 10-tree random forest classifier from the training sample.
var trainedClassifier = ee.Classifier.smileRandomForest(200,2).train({
  features: trainingSample,
  classProperty: label,
  inputProperties: bands
});

// Get information about the trained classifier.
// print('Results of trained classifier', trainedClassifier.explain());

// Classify the reflectance image from the trained classifier.
var imgClassified = composite.classify(trainedClassifier);

// GEE starts from 0
imgClassified = imgClassified.where(imgClassified.eq(0),10) // soil goes to 10, 0 is used for masking

// OPAC cannot distinguish between different structures, so classify as protected agriculture
imgClassified = imgClassified.where(imgClassified.eq(3),2)
imgClassified = imgClassified.where(imgClassified.eq(6),2)

// Add the layers to the map.
var classVis = {
  min: 1,
  max: 10,
  palette: ['#22B14C', '#FAADEE' ,'B98183', '#7E7E7E', '#FFFFFF', '#FF7980',
            '#D7ECC3', '0096a0', '#FFFFFF', '#FFD9AB']
};
// tunnel and plastic mulch reclassified as glasshouse
// forest (1:dark green), glasshouse (2:pink), tunnel (3:pale red), dark gray (4), light gray (5), salmon (6:plastic mulch),
// vegetated (7:pale green), blue (8:water), white (9:snow), baresoil (10:sand)
Map.addLayer(imgClassified, classVis, 'Classified',false);










// plot but first water mask
var waterMask = ee.ImageCollection('MODIS/006/MOD44W')
                  .filter(ee.Filter.date('2015-01-01', '2015-05-01'));
waterMask = waterMask.select('water_mask').toBands();

imgClassified = imgClassified.where(waterMask.eq(1),8) // OPAC can say that sediments in water or foam is soil. set it to water
Map.addLayer(imgClassified, classVis, 'Classified with water MODIS',false);

/*
print(waterMask)
var waterMaskVis = {
  min: 0,
  max: 1,
  palette: ['bcba99', '2d0491'],
};
Map.addLayer(waterMask, watersnMaskVis, 'Water Mask');
*/









// SNOW
var ndsiMask = composite.normalizedDifference(['B3','B11']).rename('NDSI');
// in Almeria, tunnels like snow (NDSI up to 0.25). so replace values NDSI at low altitude
var srtm = ee.Image('USGS/SRTMGL1_003');
var elevation = srtm.select('elevation');
//var slope = ee.Terrain.slope(srtm.select('elevation')).lt(10);
//Map.addLayer(slope, {min: 0, max: 60}, 'slope');

//ndsiMask = ndsiMask.where(elevation.lt(1000).and(ndsiMask.gt(0)), -99)
//imgClassified = imgClassified.updateMask(ndsiMask.lte(0));
//Map.addLayer(imgClassified, classVis, 'Classified with mask for clouds, water, NDVI and NDSI_elevation',false);

imgClassified = imgClassified.where(elevation.gt(1000).and(ndsiMask.gt(0.3)),9)
Map.addLayer(imgClassified, classVis, 'Classified with water_MODIS and NDSI_DTM',false);








// Corine mask
var corine = ee.Image('COPERNICUS/CORINE/V20/100m/2018');
//Map.addLayer(corine)
//print(corine);
//print(corine.propertyNames());
//var lc_value = corine.get('landcover_class_values');
//print(lc_value);

var agri = corine.updateMask(corine.gte(211).and(corine.lte(244))) // agriculture is in this range

// Update the composite mask with the water mask.
imgClassified = imgClassified.updateMask(agri);
Map.addLayer(imgClassified, classVis, 'Classified with water_MODIS, NDSI_DTM and Corine',false);











// Export a Cloud Optimized GeoTIFF (COG) by setting the "cloudOptimized"
// parameter to true.
var imageRGB = composite.visualize({bands: ['B4', 'B3', 'B2'], min: 0, max: 2500});
Export.image.toDrive({
  image: imageRGB,
  description: name_export_RGB.getInfo(),
  folder: 'OPAC_CATHY',
  region: ROI,
  scale: 20, // meters per pixel
  crs: 'EPSG:4326',
  maxPixels: 1e13
//  formatOptions: {
//    cloudOptimized: true
//  }
});





// Export a Cloud Optimized GeoTIFF (COG) by setting the "cloudOptimized"
// parameter to true.
// var name_export_Corine = ee.String(name_export).cat('_Corine');
// var imgClassified_exp = imgClassified.where(imgClassified.eq(0),10) // so that I can make them transparent in QGIS
Export.image.toDrive({
  image: imgClassified,
  description: name_export_Corine.getInfo(),
  folder: 'OPAC_CATHY',
  region: ROI,
  scale: 20, // meters per pixel
  crs: 'EPSG:4326',
  maxPixels: 1e13,
  formatOptions: {
    cloudOptimized: true
  }
});

// NDVI
var greenest = composite.select('NDVI');
Export.image.toDrive({
  image: greenest,
  description: name_export_NDVI.getInfo(),
  folder: 'OPAC_CATHY',
  region: ROI,
  scale: 20, // meters per pixel
  crs: 'EPSG:4326',
  maxPixels: 1e13,
  formatOptions: {
    cloudOptimized: true
  }
});









//Add NDVI band
function addNDVI(image){
  var ndvi = image.normalizedDifference(['B8','B4']).rename('NDVI');
  return image.addBands(ndvi);
}

//Add NDWI band
function addNDWI(image){
  var ndwi = image.normalizedDifference(['B8','B11']).rename('NDWI');
  return image.addBands(ndwi);
}

//Add NI band
function addNI(image){
  var ni = image.expression(
    '((2 * B8 - B4 - B12)/(2 * B8 + B4 + B12))', {
      'B8': image.select('B8'),
      'B4': image.select('B4'),
      'B12': image.select('B12')
      }).rename('NI');
  return image.addBands(ni);
}

//Add NDSI band
function addNDSI(image){
  var ndsi = image.normalizedDifference(['B3','B11']).rename('NDSI');
  return image.addBands(ndsi);
}