[GIS] Calculate MSAVI (Modified Soil-adjusted Vegetation Index) in Google Earth Engine

google-earth-enginejavascriptremote sensingsavivegetation-index

How do you write the MSAVI formula in JavaScript function for Google Earth Engine?

I am trying to calculate the MSAVI as a function. The aim is to map the result over a Sentinel-2 composite and not a single image.

I tried the same thing for the OSAVI Index and it worked fine:
OSAVI = (NIR-Red)/ (NIR+Red+0.16)

function addOSAVI(input) {
  var nir = input.select(['B8']);
  var red = input.select(['B4']);
  var osavi = nir.subtract(red).divide(nir.add(red).add(0.16)).rename('OSAVI');
return input.addBands(osavi);

}

I need a function in the same manner as the OSAVI function, but I have problem with implementing sqrt() and pow() correctly.

enter image description here

function addMSAVI(input) {
  var nir = input.select(['B8']);
  var red = input.select(['B4']);
  var msavi = (nir.add(1).multiply(2).subtract(sqrt(pow(nir.add(1).multiply(2),2)).subtract((nir.subtract(red)).multiply(2)))).divide(2).rename('MSAVI');
return input.addBands(msavi);

}

More information about MSAVI:
http://wiki.landscapetoolbox.org/doku.php/remote_sensing_methods:modified_soil-adjusted_vegetation_index

Best Answer

Both ways work just fine and give the same results: https://code.earthengine.google.com/3f766772035a8890b0caf231eb652a1f.

There is also a small mistake in the formula in your code - wrong bracket.

// compute composite for all images aquired in the last 3 days
var now = ee.Date(Date.now())
var composite = s2
  .filterDate(now.advance(-3, 'day'), now)
  .mosaic()
  .divide(10000)

Map.addLayer(composite, {bands: ['B8', 'B8', 'B4'], min: 0, max: 0.4}, 'composite')

// compute MSAVI2 using expression
var msavi2 = composite.expression(
  '(2 * NIR + 1 - sqrt(pow((2 * NIR + 1), 2) - 8 * (NIR - RED)) ) / 2', 
  {
    'NIR': composite.select('B8'), 
    'RED': composite.select('B4')
  }
);

Map.addLayer(msavi2, {min: -0.1, max: 0.5}, 'msavi2 (expression)')

// compute MSAVI2 using fluent API
msavi2 = composite.select('B8').multiply(2).add(1)
  .subtract(composite.select('B8').multiply(2).add(1).pow(2)
    .subtract(composite.select('B8').subtract(composite.select('B4')).multiply(8)).sqrt()
  ).divide(2)

Map.addLayer(msavi2, {min: -0.1, max: 0.5}, 'msavi2 (fluent)')
Related Question