Machine Learning-Supervised classification

/*
Copyright (c) 2022 Desmond Lartey.
This work is licensed under the terms of the MIT license.  
For a copy, see https://opensource.org/licenses/MIT
*/
var s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED');
var urbanAreas = ee.FeatureCollection('projects/ee-desmond/assets/ne_10m_urban_areas');
var alos = ee.Image('JAXA/ALOS/AW3D30/V2_2');
//Map.addLayer(urbanAreas)

// Perform supervised classification for your city
// Find the feature id by adding the layer to the map and using Inspector.
var city = urbanAreas.filter(ee.Filter.eq('system:index', '00000000000000001ec5'));
var geometry = city.geometry();
Map.centerObject(geometry);


var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
  .filter(ee.Filter.date('2019-01-01', '2020-01-01'))
  .filter(ee.Filter.bounds(geometry))
  .select('B.*');

// Function to remove cloud and snow pixels from Sentinel-2 SR image

function maskCloudAndShadowsSR(image) {
  var cloudProb = image.select('MSK_CLDPRB');
  var snowProb = image.select('MSK_SNWPRB');
  var cloud = cloudProb.lt(10);
  var scl = image.select('SCL'); 
  var shadow = scl.eq(3); // 3 = cloud shadow
  var cirrus = scl.eq(10); // 10 = cirrus
  // Cloud probability less than 10% or cloud shadow classification
  var mask = cloud.and(cirrus.neq(1)).and(shadow.neq(1));
  return image.updateMask(mask);
}


var composite = filtered.median().clip(geometry);

// Function to add Indices
var addIndices = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename(['ndvi']);
  var ndbi = image.normalizedDifference(['B11', 'B8']).rename(['ndbi']);
  var mndwi = image.normalizedDifference(['B3', 'B11']).rename(['mndwi']); 
  var bsi = image.expression(
      '(( X + Y ) - (A + B)) /(( X + Y ) + (A + B)) ', {
        'X': image.select('B11'), //swir1
        'Y': image.select('B4'),  //red
        'A': image.select('B8'), // nir
        'B': image.select('B2'), // blue
  }).rename('bsi');
  return image.addBands(ndvi).addBands(ndbi).addBands(mndwi).addBands(bsi);
};

var composite = addIndices(composite);


// Calculate Slope and Elevation
var elev = alos.select('AVE_DSM').rename('elev');
var slope = ee.Terrain.slope(alos.select('AVE_DSM')).rename('slope');

var composite = composite.addBands(elev).addBands(slope);

var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000, gamma: 1.2};
Map.addLayer(composite, visParams, 'RGB');

// Normalize the image 

// Machine learning algorithms work best on images when all features have
// the same range

// Function to Normalize Image
// Pixel Values should be between 0 and 1
// Formula is (x - xmin) / (xmax - xmin)
//************************************************************************** 
function normalize(image){
  var bandNames = image.bandNames();
  // Compute min and max of the image
  var minDict = image.reduceRegion({
    reducer: ee.Reducer.min(),
    geometry: geometry,
    scale: 10,
    maxPixels: 1e9,
    bestEffort: true,
    tileScale: 16
  });
  var maxDict = image.reduceRegion({
    reducer: ee.Reducer.max(),
    geometry: geometry,
    scale: 10,
    maxPixels: 1e9,
    bestEffort: true,
    tileScale: 16
  });
  var mins = ee.Image.constant(minDict.values(bandNames));
  var maxs = ee.Image.constant(maxDict.values(bandNames));

  var normalized = image.subtract(mins).divide(maxs.subtract(mins));
  return normalized;
}

// Display the input composite.

var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
Map.addLayer(composite, rgbVis, 'image');

// Exercise
// Add training points for 4 classes
// Assign the 'landcover' property as follows

// urban: 0
// bare: 1
// water: 2
// vegetation: 3

// After adding points, uncomments lines below

var gcp = urban.merge(bare).merge(water).merge(vegetation);

// Overlay the point on the image to get training data.
var training = composite.sampleRegions({
  collection: gcp, 
  properties: ['landcover'], 
  scale: 10,
  tileScale: 16
});
print(training);


// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50).train({
  features: training,  
  classProperty: 'landcover', 
  inputProperties: composite.bandNames()
});
// // Classify the image.
var classified = composite.classify(classifier);
Map.addLayer(classified, {min: 0, max: 3, palette: ['red', 'brown', 'blue', 'green']}, '2019'); 

var gcp = gcp.randomColumn();

var trainingGcp = gcp.filter(ee.Filter.lt('random', 0.6));
var validationGcp = gcp.filter(ee.Filter.gte('random', 0.6));

var training = composite.sampleRegions({
  collection: trainingGcp,
  properties: ['landcover'],
  scale: 10,
  tileScale: 16
});

// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50)
.train({
  features: training,  
  classProperty: 'landcover',
  inputProperties: composite.bandNames()
});

// Classify the image.
var classified = composite.classify(classifier);

//************************************************************************** 
// Accuracy Assessment
//************************************************************************** 

// Use classification map to assess accuracy using the validation fraction
// of the overall training set created above.
var test = classified.sampleRegions({
  collection: validationGcp,
  properties: ['landcover'],
  tileScale: 16,
  scale: 10,
});

var testConfusionMatrix = test.errorMatrix('landcover', 'classification');
print('Confusion Matrix', testConfusionMatrix);
print('Test Accuracy', testConfusionMatrix.accuracy());

print(testConfusionMatrix.producersAccuracy())
print(testConfusionMatrix.consumersAccuracy())
print(testConfusionMatrix.kappa())

// identify area do each of the classes

var vegetation = classified.eq(3);
var areaImage = vegetation.multiply(ee.Image.pixelArea());

var area = areaImage.reduceRegion({
  reducer: ee.Reducer.sum(),//.group({
    //groupField: 1, 
    //groupName: 'classified'),
  geometry: geometry,
  scale: 10,
  maxPixels: 1e10
});

var vegetationAreaSqKm = ee.Number(area.get('classification')).divide(1e6).round();

// Compute and Print the percentage of green cover of the city

var cityAreaSqKm = geometry.area()
var areaPercent = (vegetationAreaSqKm.divide(cityAreaSqKm)).multiply(100).round();
print(areaPercent)