0

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)
2
  • 1
    Please provide a script where you're trying to solve this. – Daniel Wiell May 19 '20 at 8:13
  • Hi Daniel, I added my coding efforts. I end up getting a feature with the NDVI values for each available dates. But then I need to "select" the original Sentinel2 image with the highest NDVI average value. I added the date as a band, but I got really lost with this last issue – Alfonso de Lara May 19 '20 at 16:50
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()
    )
}
4
  • Outstanding! I was so far from doing this.. Pretty hard to understand your code, but I tried it and compare it manually and does the job. Thank you so much – Alfonso de Lara May 20 '20 at 17:58
  • I'm trying to do some tune ups for this script, but I found this is like "backwards". You declare variables after calling them. Could you explain how that works? Or maybe just re-organize it making it easier to understand. Cheers – Alfonso de Lara Jun 14 '20 at 18:12
  • I use the functions before declaring them so the code goes from abstract to more specific, as you read a file. You can get a general idea of an algorithm/workflow by a composition of functions with descriptive names. Then you can delve into each function, to find out how each step of the algorithm works in detail. – Daniel Wiell Jun 16 '20 at 9:22
1

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)
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.