var fromList = ee.FeatureCollection("FAO/GAUL/2015/level0");
var name_export = 'U.K. of Great Britain and Northern Ireland';
var ROI = fromList.filter(ee.Filter.inList('ADM0_NAME', [name_export]));


// OR
/*
// Region of interest with problems in FAO. Estonia
var ROI = ee.Geometry.Polygon(
	[[[22.00481269779213,57.849968872150654],
[23.27922676029213,57.7269883335885],
[24.41081855716713,57.80316854346273],
[25.24577949466713,57.931723313811595],
[26.56413886966713,57.43248781612245],
[27.44304511966713,57.47976815199024],
[28.00334785404213,57.83827448647697],
[27.55290840091713,58.756192595531374],
[28.37688301029213,59.3885767729706],
[28.00334785404213,59.55042601983533],
[25.68523261966713,59.72258143280693],
[24.44377754154213,59.6726917680118],
[21.88396308841713,58.989042848638725],
[21.62029121341713,58.33194683228729],
[22.00481269779213,57.849968872150654]]], null, false);
var name_export = 'Estonia';
*/

////////////////////////////////////////////////////////////////////////////////////////////////////

//print(ROI);
//Map.setCenter(25.1, 45.5, 6);
//Map.addLayer(ROI, {}, 'Regions');


var table = ee.FeatureCollection("users/ee-dlacecilia-unipd/Validation");

Map.centerObject(table, 3);
Map.addLayer(table, {color: 'black'}, 'myPOI');

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


// Region of interest. Austria, Hungary, Romania, Slovenia, Croatia, Serbia, Bulgaria, Albania, Montenegro, Kosovo, North Macedonia
//var ROI = filteredArea;
//var name_export = 'GEE_UK';
var name_export = ee.String(name_export).cat('_GEE');
var name_export_Corine = ee.String(name_export).cat('_OPAC_Corine');
var name_export_RGB = ee.String(name_export).cat('_RGB');
//var table = ee.FeatureCollection("users/ee-dlacecilia-unipd/a_DF_validation_11_copy"); // Daniele









// Date for mapping with random forest
var start_date = '2021-05-01'
var end_date   = '2021-08-31'



// 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))
    })
    .select(['B3','B1','B2','B4','B5','B6','B7','B8','B11','B12']);


var addDate = function(image) {
  // return image.addBands(ee.Image.constant(image.date().millis()).long().rename('time'))
  return image.addBands(ee.Image.constant(ee.Number.parse(image.date().format("YYYYMMdd"))).rename('time').float()); 
};

var composite = composite.map(addDate)
.select(['B3','time','B1','B2','B4','B5','B6','B7','B8','B11','B12']); // green band in the first position



//var compositesave = composite
// turn image collection into an array
var array = composite.toArray()

// sort array by the first band, keeping other bands
var axes = { image:0, band:1 }

var sort = array.arraySlice(axes.band, 0, 1);  // select bands from index 0 to 1 (NDVI)

var sorted = array.arraySort(sort);


// take the first image only
// for the min value sorted
var values = sorted.arraySlice(axes.image, 0, 1);
print(values)
// convert back to an image  
composite = values.arrayProject([axes.band]).arrayFlatten([['B3','time','B1','B2','B4','B5','B6','B7','B8','B11','B12']])
composite = composite.select(['B1','B2','B3','B4','B5','B6','B7','B8','B11','B12'])
//print(max)
// add to map
// Map.addLayer(max, {}, 'NDVI and time stamp')




// clip for country statistics
composite = composite.clip(ROI)




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

Map.addLayer(composite, s2Viz, 'median composite',false);
Map.centerObject(ROI, 6);
//Map.addLayer(ROI, {color: 'red'}, false)
//Map.addLayer(compositesave, s2Viz, 'median composite'); // this is to check in the inspector that the minimum B3 is taken and all band values of this time snapshot are part of this
//Map.addLayer(table);


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








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
// if water at high altitude, sometimes went to snow. so avoid reclassifying water pixels
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)).and(imgClassified.neq(8)),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');










/*
// 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: 'OPAC2',
  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: 'OPAC2',
  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);
}



/*


imgClassified = imgClassified.eq(2);
// Mask 0 pixel values and count remaining pixels.
imgClassified = imgClassified.selfMask();

var statsMasked = imgClassified.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: ROI,
    scale: 20,
    maxPixels: 1e13,
    tileScale: 12 // not enough
    });
var imgClassifiedPixels = statsMasked.get('classification');
print(imgClassifiedPixels);

*/

/*
// memory intensive
var pixelCountStats = imgClassified.reduceRegion({
    reducer: ee.Reducer.frequencyHistogram().unweighted(),
    geometry: ROI,
    scale: 4000,
    maxPixels: 1e13
    });

var pixelCounts = ee.Dictionary(pixelCountStats.get('classification'));
print(pixelCounts);
*/

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

/*
// Format the results to make it more readable.
var classLabels = ee.List([
    'forest', 'glasshouse', 'tunnel', 'dark', 'light',
    'mulch', 'vegetated', 'water', 'snow', 'baresoil'
    ]);
*/

/*
// Rename keys with class names.
var pixelCountsFormatted = pixelCounts.rename(
  pixelCounts.keys(), classLabels);
print(pixelCountsFormatted);
*/

/*
// Create a Feature Collection.
var exportFc = ee.FeatureCollection(
  //ee.Feature(null, pixelCountsFormatted));
  ee.Feature(null, pixelCounts));

// Export the results as a CSV file.
Export.table.toDrive({
  collection: exportFc,
  description: 'pixel_counts_export',
  folder: 'OPAC2',
  fileNamePrefix: 'pixel_counts',
  fileFormat: 'CSV',
});
*/