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

Advertisements

2 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

    Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s