Sample your composite to create the land cover map

Create training data, calculate covariates and apply machine learning

Step 1: import the image you exported in the previous step and add it to the map

Map.addLayer(image,{min:0,max:0.6,bands:"??,??,??"})

Step 2: Create a new feature collection for water and add markers on water pixels

firstImage8.png

Step 3: Do the same for forest
firstImage9.png

Step 4: Add as many categories as you like

Step 5: Combine all data

var TrainingData = water.merge(forest).merge(..)

Save the script

Step 6: Open a new tab and create a script called covariate_module with the code below.

var elevation = ee.Image("USGS/SRTMGL1_003")
var jrcImage = ee.Image("JRC/GSW1_0/GlobalSurfaceWater")

var ndCovariatesList = [
  ['blue', 'green'],
  ['blue', 'red'],
  ['blue', 'nir'],
  ['blue', 'swir1'],
  ['blue', 'swir2'],
  ['green', 'red'],
  ['green', 'nir'],
  ['green', 'swir1'],
  ['green', 'swir2'],
  ['red', 'swir1'],
  ['red', 'swir2'],
  ['nir', 'red'],
  ['nir', 'swir1'],
  ['nir', 'swir2'],
  ['swir1', 'swir2']
];

var rCovariatesList = [
  ['swir1', 'nir'],
  ['red', 'swir1']
];

var ComputeNDCovariatesList = function () {
  var list = [];
  for (var index in ndCovariatesList) {
    var list_ = [ndCovariatesList[index][0], ndCovariatesList[index][1]];
    list.push(list_);
  }
  return list;
};

var addNDCovariates = function (image){
  var list = ComputeNDCovariatesList();
  print(list)
  for (var index in list) {
    image = image.addBands(image.normalizedDifference(list[index]).rename('ND_'+ ndCovariatesList[index][0] + '_' + ndCovariatesList[index][1]));
  }
  return image;
};

var ComputeRCovariatesList = function () {
  var list = [];
  for (var index in rCovariatesList) {
    var list_ = [rCovariatesList[index][0], rCovariatesList[index][1]];
    list.push(list_);
  }
  return list;
};

var addRCovariates = function (image) {
  var list = ComputeRCovariatesList();
  for (var index in list) {
    image = image.addBands(image.select(list[index][0]).divide(image.select(list[index][1]))
            .rename('_R_' + rCovariatesList[index][0] + '_' + rCovariatesList[index][1]));
  }
  return image;
};

var addEVI = function (image) {
  var evi = image.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
    'NIR' : image.select('nir'),
    'RED' : image.select('red'),
    'BLUE': image.select('blue')
  }).float();
  return image.addBands(evi.rename('EVI'));
};

var addSAVI = function (image) {
  // Add Soil Adjust Vegetation Index (SAVI)
	// using L = 0.5;
	var savi = image.expression('(NIR - RED) * (1 + 0.5)/(NIR + RED + 0.5)', {
    'NIR': image.select('nir'),
    'RED': image.select('red')
	}).float();
	return image.addBands(savi.rename('SAVI'));
};

var addIBI = function (image) {
  // Add Index-Based Built-Up Index (IBI)
  var ibiA = image.expression('2 * SWIR1 / (SWIR1 + NIR)', {
    'SWIR1': image.select('swir1'),
    'NIR'  : image.select('nir')
  }).rename(['IBI_A']);

  var ibiB = image.expression('(NIR / (NIR + RED)) + (GREEN / (GREEN + SWIR1))', {
    'NIR'  : image.select('nir'),
    'RED'  : image.select('red'),
    'GREEN': image.select('green'),
    'SWIR1': image.select('swir1')
  }).rename(['IBI_B']);

  var ibiAB = ibiA.addBands(ibiB);
  var ibi = ibiAB.normalizedDifference(['IBI_A', 'IBI_B']);
  return image.addBands(ibi.rename(['IBI']));
};

// Function to compute the Tasseled Cap transformation and return an image
var getTassledCapComponents = function (image) {
  var coefficients = ee.Array([
    [0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863],
    [-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800],
    [0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572],
    [-0.8242, 0.0849, 0.4392, -0.0580, 0.2012, -0.2768],
    [-0.3280, 0.0549, 0.1075, 0.1855, -0.4357, 0.8085],
    [0.1084, -0.9022, 0.4120, 0.0573, -0.0251, 0.0238]
  ]);
  var bands = ee.List(['blue', 'green', 'red', 'nir', 'swir1', 'swir2']);

  // Make an Array Image, with a 1-D Array per pixel.
  var arrayImage1D = image.select(bands).toArray();

  // Make an Array Image with a 2-D Array per pixel, 6 x 1
  var arrayImage2D = arrayImage1D.toArray(1);

  var componentsImage = ee.Image(coefficients).matrixMultiply(arrayImage2D).arrayProject([0])
                        .arrayFlatten([['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth']]).float();
  // Get a multi-band image with TC-named bands
  return image.addBands(componentsImage);
};

// Function to add Tasseled Cap angles and distances to an image. Assumes image has bands: 'brightness', 'greenness', and 'wetness'.
var getTassledCapAngleAndDistance = function (image) {
  var brightness = image.select('brightness');
  var greenness = image.select('greenness');
  var wetness = image.select('wetness');

  // Calculate tassled cap angles and distances
  var tcAngleBG = brightness.atan2(greenness).divide(Math.PI).rename(['tcAngleBG']);
  var tcAngleGW = greenness.atan2(wetness).divide(Math.PI).rename(['tcAngleGW']);
  var tcAngleBW = brightness.atan2(wetness).divide(Math.PI).rename(['tcAngleBW']);

  var tcDistanceBG = brightness.hypot(greenness).rename(['tcDistanceBG']);
  var tcDistanceGW = greenness.hypot(wetness).rename(['tcDistanceGW']);
  var tcDistanceBW = brightness.hypot(wetness).rename(['tcDistanceBW']);

  image = image.addBands(tcAngleBG).addBands(tcAngleGW).addBands(tcAngleBW).addBands(tcDistanceBG).addBands(tcDistanceGW).addBands(tcDistanceBW);

  return image;
};

var computeTassledCap = function (image) {
  image = getTassledCapComponents(image);
  image = getTassledCapAngleAndDistance(image);
  return image;
};

var addTopography = function (image) {

  // Calculate slope, aspect and hillshade
  var topo = ee.Algorithms.Terrain(elevation);

  // From aspect (a), calculate eastness (sin a), northness (cos a)
  var deg2rad = ee.Number(Math.PI).divide(180);
  var aspect = topo.select(['aspect']);
  var aspect_rad = aspect.multiply(deg2rad);
  var eastness = aspect_rad.sin().rename(['eastness']).float();
  var northness = aspect_rad.cos().rename(['northness']).float();

  // Add topography bands to image
  topo = topo.select(['elevation','slope','aspect']).addBands(eastness).addBands(northness);
  image = image.addBands(topo);
  return image;
};

var addJRCDataset = function (image) {
  // Update the mask.
  jrcImage = jrcImage.unmask(0);

  image = image.addBands(jrcImage.select(['occurrence']).rename(['occurrence']));
  image = image.addBands(jrcImage.select(['change_abs']).rename(['change_abs']));
  image = image.addBands(jrcImage.select(['change_norm']).rename(['change_norm']));
  image = image.addBands(jrcImage.select(['seasonality']).rename(['seasonality']));
  image = image.addBands(jrcImage.select(['transition']).rename(['transition']));
  image = image.addBands(jrcImage.select(['max_extent']).rename(['max_extent']));

  return image;
};
exports.addCovariates = function (image) {
  image = addNDCovariates(image);
  image = addEVI(image);
  image = addSAVI(image);
  image = addIBI(image);
  image = computeTassledCap(image);
  return image;
};

var addJRCAndTopo = function (image) {
  image = addTopography(image);
  image = addJRCDataset(image);
  return image;
};

exports.addJRCAndTopo = addJRCAndTopo;

firstImage11.png

Step 6: go back to the previous script and add the code below to create your landcover map

var covariates = require("users/servirmekong/mrc:covariate_module");
var cambodia = ee.FeatureCollection("users/servirmekong/countries/KHM_adm1");
var KampongThum = cambodia.filter(ee.Filter.eq("NAME_1","Kâmpóng Thum")).geometry();

image = covariates.addJRCAndTopo(image);
image = covariates.addCovariates(image);

var trainingSample = image.sampleRegions(TrainingData,["land_class"],10);

var bandNames = image.bandNames();
var classifier = ee.Classifier.randomForest(30,0).train(trainingSample,"land_class",bandNames);

var classification = image.classify(classifier,'Mode');

Map.addLayer(classification.clip(KampongThum),{min:0,max:5,palette:"aec3d4,152106,387242,cc0013,3bc3b2,8dc33b"},"landcover")

example

4 comments

  1. This script was very interesting and well done, I never see one divided in multiple parts. In this way it runs smoother.
    I successfully applied it to my area and the results are good!. The only problem is that the classifier is confusing between
    the bare soils and urban areas, do you have a suggestion to improve the final accuracy of the script?
    My idea is to use other indices other than the original one.
    You can see my results here: https://code.earthengine.google.com/ef9ec929ca9571704646117db4cdc05f
    Thank you and good job!
    Alex

  2. Thank you very much for the code! Unfortunately, the “.randomForest” function is deprecated, I’m trying to use “.smileRandomForest” but I’m not an expert on this, could you help us how to solve this?

Leave a Reply