Google Earth Engine - Get Image with highest max Scene NDVI from a time period

Question

I have a shapefile with several fields, corresponding to different wheat fields of a certain year. Using as Input an image collection converted to NDVIs, I would like to calculate the mean NDVI value for each field for every available date. Then, export the raster of the maximum "mean value" of each field. The key here is not to do a pixel-wise maximum NDVI selection, but by FIELD.

In other words, I need to get the Sentinel-2 NDVI image from a field that corresponds to the highest mean NDVI value (of the field polygon) registered within a time period. This is not the max() function over a ImageCollection (which synthesizes a image with pixels from different dates). I need the image of a REAL scene, which actually happened at one DATE.  

var geometry = /* color: #d63000 */ee.FeatureCollection(
    [ee.Feature(
        ee.Geometry.Polygon(
            [[[-61.58584513408002, -35.77369642524422],
              [-61.58035197001752, -35.778222714936696],
              [-61.5727130387431, -35.772373307284425],
              [-61.56704821330365, -35.77683003785912],
              [-61.55262865763959, -35.76506094262748],
              [-61.5643016312724, -35.75586730183619],
              [-61.58472933512982, -35.77286077435719]]]),
        {
          "system:index": "0"
        }),
    ee.Feature(
        ee.Geometry.Polygon(
            [[[-60.64656168176598, -35.72030669380027],
              [-60.62716394617028, -35.70358079489495],
              [-60.645703374881215, -35.68963986472292],
              [-60.652913152713246, -35.69214941201044],
              [-60.66733270837731, -35.70427777742714]]]),
        {
          "system:index": "1"
        })]);

Map.addLayer(geometry,{},'ROIs')


function getNDVI10 (image){
  return image.normalizedDifference(['B8','B4'])
              .select([0],['NDVI'])
              .addBands(image.metadata('system:time_start').rename('DATE'))
}
function getStats (image){
  return image.select(['NDVI']).reduceRegions({collection: geometry,reducer: ee.Reducer.mean(),scale: 10,
})}

var nbrParams = {min: 0.4, max: 0.8, 'palette': ['red','yellow','green']};
var s2 = ee.ImageCollection("COPERNICUS/S2")
  .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))
  .filterBounds(geometry);
var NDVI10_2017 = s2.filterDate('2017-10-01','2017-12-01').map(getNDVI10)

////////////////////////////////////////////////////////////////////////////////////////////////////////////
var NDVIpixel_2017 = NDVI10_2017.qualityMosaic('NDVI').clip(geometry)
Map.addLayer(NDVIpixel_2017.select(['NDVI']),{}, "maxNDVI S2",true)
////////////////////////////////////////////////////////////////////////////////////////////////////////////

var NDVI10_2017_stats = NDVI10_2017.map(getStats)

Answer 1

If I understood you correct, the below script might do what you want. It creates an image per field, with the NDVI from the date with the highest mean NDVI. Each image is created like this:

1. Create daily mosaics
2. Calculate mean for each mosaic
3. Return the mosaic with the highest mean

https://code.earthengine.google.com/013a113d3f01c192ff7370d236197a8a

var startDate = ee.Date('2017-10-01')
var endDate = ee.Date('2017-12-01')
var deltaDays = 1
var fields = ee.FeatureCollection(
  [ee.Feature(
      ee.Geometry.Polygon(
        [
          [
            [-61.58584513408002, -35.77369642524422],
            [-61.58035197001752, -35.778222714936696],
            [-61.5727130387431, -35.772373307284425],
            [-61.56704821330365, -35.77683003785912],
            [-61.55262865763959, -35.76506094262748],
            [-61.5643016312724, -35.75586730183619],
            [-61.58472933512982, -35.77286077435719]
          ]
        ]), {
        "system:index": "0"
      }),
    ee.Feature(
      ee.Geometry.Polygon(
        [
          [
            [-60.64656168176598, -35.72030669380027],
            [-60.62716394617028, -35.70358079489495],
            [-60.645703374881215, -35.68963986472292],
            [-60.652913152713246, -35.69214941201044],
            [-60.66733270837731, -35.70427777742714]
          ]
        ]), {
        "system:index": "1"
      })
  ])

var greenest = ee.ImageCollection(
  fields.map(greenestDateMosaic)
)

print(greenest)
Map.addLayer(greenest)
Map.centerObject(fields)

function greenestDateMosaic(field) {
  var mosaics = createMosaics(field)
  var mosaicsWithMean = mosaics.map(function (image) {
    return setMean(field, image)
  })
  var highestMean = mosaicsWithMean.aggregate_max('mean')
  return  mosaicsWithMean
    .filterMetadata('mean', 'equals', highestMean)
    .first()
    .clip(field.geometry())
}

function setMean(field, image) {
  var mean = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: field.geometry(),
    scale: 10,
    maxPixels: 1e13
  }).get('ndvi')
  return image
    .set('mean', mean)
}

function createMosaics(field) {
  var days = endDate.difference(startDate, 'days')
  var dateOffsets = ee.List.sequence(0, days.subtract(1), deltaDays)
  return ee.ImageCollection(dateOffsets
    .map(function (dateOffset) {
      return createMosaic(field, dateOffset)
    }))
    .filterMetadata('empty', 'equals', 0)  
}

function createMosaic(field, dateOffset) {
  var date = startDate.advance(dateOffset, 'days')
  var mosaic = ee.ImageCollection("COPERNICUS/S2")
    .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))
    .filterDate(date, date.advance(deltaDays, 'days'))
    .filterBounds(field.geometry())
    .map(toNdvi)
    .mosaic()  
  return mosaic
    .set('system:time_start', date.millis())
    .set('empty', mosaic.bandNames().size().eq(0))
}

function toNdvi(image) {
  return image
    .normalizedDifference(['B8', 'B4'])
    .rename('ndvi')
    .updateMask(
      image.select('QA60').not()
    )
}

Answer 2

I think you want the sentinel image/date which returns the highest NDVI for each region. This is how I would do it:

This algorithm first maps over all your geometries, and filters the bounds of the imageCollection to it. Then it maps over all Images and calculates the NDVI for the current geometry. Finally it sorts the feature Collection by highest NDVI and returns only that feature. That leaves you with one Feature per Geometry with NDVI and the corresponding Sentinel Image for the NDVI value.

Code

function getNDVI10 (image){
  return image.normalizedDifference(['B8','B4'])
              .select([0],['NDVI'])
              .set('system:time_start', image.get('system:time_start'))
}
function getStats (image){
  return image.select(['NDVI']).reduceRegions({collection: geometry,reducer: ee.Reducer.mean(),scale: 10,
})}

var nbrParams = {min: 0.4, max: 0.8, 'palette': ['red','yellow','green']};
var s2 = ee.ImageCollection("COPERNICUS/S2")
  .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))
  .filterBounds(geometry);
var NDVI10_2017 = s2.filterDate('2017-10-01','2017-12-01').map(getNDVI10)

////////////////////////////////////////////////////////////////////////////////////////////////////////////
var NDVIpixel_2017 = NDVI10_2017.qualityMosaic('NDVI').clip(geometry)
Map.addLayer(NDVIpixel_2017.select(['NDVI']),{}, "maxNDVI S2",true)
////////////////////////////////////////////////////////////////////////////////////////////////////////////

var NDVI10_2017_stats = NDVI10_2017.map(getStats)

var highestNDVI = geometry.map(function(geom){
  var singleGeom = ee.Feature(geom).geometry()
  var filterCollection = NDVI10_2017.filterBounds(singleGeom)

  var meanFeatureCollection = filterCollection.map(function(image){
    var means = image.reduceRegion({
        reducer: ee.Reducer.mean(),
        geometry: singleGeom,
        scale: 10
      })
    var feature = ee.Feature(singleGeom, means)
    return feature.set({'system:time_start': image.get('system:time_start'),
      'sentinel_id': image.get('system:index')})
  })

  return ee.Feature(meanFeatureCollection.sort('mean').first())
})

print(highestNDVI)

Map.addLayer(highestNDVI, {}, "Date maxNDVI S2",true)