Extracting index/position of maximum value in annual time series in Google Earth Engine?

Question

I am trying to extract the integer index/position of the maximum value in an annual MODIS NDVI time series in Earth Engine but I am having a difficult time to figure out a plan to do this. It is quite easy to calculate the maximum value of the image collection using the following javascript code:

var bbox = ee.Geometry.Polygon([
  [-8.44299,42.20410],
  [-8.43200,41.65239],
  [-7.75634,41.65239],
  [-7.76733,42.20817],
  [-8.44299,42.20410]]);

var getNDVI = function(img){
  return img.normalizedDifference(['sur_refl_b02','sur_refl_b01']);
};

var MOD = ee.ImageCollection('MODIS/006/MOD09A1').
          filterDate('2015-01-01','2015-12-31').
          filterBounds(bbox);

var NDVI = MOD.map(getNDVI);

var NDVImax = NDVI.max();

For example, in R raster package this would correspond to an operation with calc() or stackApply() with which.max() as a function to retrieve the position of the maximum value in the raster stack time series. However, I would prefer to do this directly in GEE without having to download the data.

Suggestions?

Answer 1

A bit tricky, but you can use arrays to achieve this: https://code.earthengine.google.com/7f79588ccca7a31b9aeee7710aa1779e. For example, you can add a timestamp band to associate NDVI values with images and then sort array to find maximum values, keeping both NDVI and timestamp bands together.

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

// sort array by the first band (NDVI)
var axes = { image:0, band:1 }
var bandIndex = 0 // index of NDVI band
var sort = array.arraySlice(axes.band, bandIndex, bandIndex + 1); 
var sorted = array.arraySort(sort);

// take the last image only (slice)
var length = sorted.arrayLength(axes.image)
var values = sorted.arraySlice(axes.image, length.subtract(1), length);

// convert back to an image  
var max = values.arrayProject([axes.band]).arrayFlatten([['ndvi', 'time']])

// visualize
Map.addLayer(max, {}, 'NDVI and time stamp', false)

var ndviMax = max.select(0)
var time = max.select(1)
Map.addLayer(ndviMax, {min: 0, max: 1}, 'NDVI' )
Map.addLayer(time.randomVisualizer(), {}, 'time stamp')

Note, that the maximum image value may vary from pixel to pixel. Also, it is better to think about image collection as a set, and not as a list, that is why I've used timestamps instead of image indices (non-existing). Image order may vary unless you sort image collection.

In case if you need to find only a single image with the overall maximum value for a given geometry, this becomes a bit easier: https://code.earthengine.google.com/2536c8d3ef38539de8274054e6c8fa3f

// compute maximum value for a given polygon for all images
NDVI = NDVI.map(function(img) {
  var max = img.reduceRegion({ reducer: ee.Reducer.max(), geometry: bbox, scale: 300 })
  return img.set({max: max.get('ndvi')})
})

// find image with a maximum value
var NDVImax = ee.Image(NDVI.sort('max', false).first())

// visualize
print('date of an image with the maximum value over bbox: ', NDVImax.date())

Map.addLayer(NDVImax, {min: 0, max: 1}, 'max')

However, in the latter case, results can be confusing, because there are multiple images with the same max value, so it might be better to use some high percentile for regional reduction instead of max.

Answer 2

I used a rather simple approach: https://code.earthengine.google.com/867c57e5f893df64e81e3caa683568d7

In this case, I am first computing an image with the minimum NDPI values over 2 months, then creating an image that shows the index of the image that produced a minimal value of the NDPI using .eq(), and then use .updateMask() to create a new mosaic containing all the bands that lead to the minimum value at a given pixel:

var index_min_img = imgs.map(function (i) {
var index_nr = ee.Number(i.get('index_nr'));
return ee.Image(0).where(i.select('ndpi').eq(min_img), index_nr).float();
}).max();

var mosaic_nr_list = ee.List.sequence(1, imgs.size(), 1).map(function (nr) {
var nr_list = ee.Number(nr).subtract(1);
var mosaic_nr = ee.Image(list.get(nr_list)).updateMask(index_min_img.eq(ee.Number(nr)));
return mosaic_nr;
});
var mosaic_responsible = ee.ImageCollection.fromImages(mosaic_nr_list).mosaic();