GCI and NDVI Uncorrelated in Same Region and Time Period

Question

I am attempting to calculate NDVI and GCI values during the 2019 growing season within districts in Mali. Specifically, I want to take the max of NDVI/GCI values for each pixel over the growing season for each district and then average all of these maximum values within each district to get a "mean of the max" for each district. I am using MODIS for NDVI values and LANDSAT for GCI. While I know the resolutions are different I am getting almost no correlation between the values for different districts, which seems odd given that the scholarly articles I've read say LANDSAT and MODIS are very similar. I know NDVI measures the fraction of red light reflected and GCI is more a measure of chlorophyll, but I would expect some correlation. Here is my code for each measure. Any explanation relating to coding errors, or the nature of NDVI and GCI?

NDVI

 var districts = ee.FeatureCollection("users/geerootfold/gadm36_MLI_4")
  .filter(ee.Filter.eq('NAME_0', 'Mali'))



var year = 2019

 
 var modis = ee.ImageCollection('MODIS/006/MOD13Q1')
    .filterDate('2019-05-01','2019-09-30')
    .select('NDVI')
    
var mosaic = modis.max().divide(10000).clip(districts)

// Add reducer output to the Features in the collection.
var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});


var featureCollection = mosaic.reduceRegions({
  collection: districts,
  reducer: reducers,
  scale: 250,
});

print(featureCollection)


var size=featureCollection.size()
var toList = featureCollection.toList(size)
print(toList, 'toList')

var idList = ee.List.sequence(0, size.subtract(1))


var newList = idList.map(function(x){
  
  var index = idList.get(x)
  var feat = toList.get(x)
  
  return ee.Feature(feat).set('year', 2019).copyProperties(feat)
})

print(newList, 'newList')

var MeansOfFeatures = ee.FeatureCollection(newList)


Export.table.toDrive({
  collection: MeansOfFeatures, 
  description: 'meanOfMaxNdvi'+year,
  folder: 'GEE NDVI2',
  // Explicitly specify your columns to exclude .geo
  fileFormat: 'CSV',
  selectors: ['NAME_4', 'year', 'mean','stdDev'] 
})


 var districts = ee.FeatureCollection("users/geerootfold/gadm36_MLI_4")
  .filter(ee.Filter.eq('NAME_0', 'Mali'))
  
  
  var colour = {
 min: 0.0,
 max: 1.0,
 palette: [
 'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
 '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
 '012E01', '011D01', '011301'
 ],
};

Map.addLayer(mosaic, colour, 'spatial mosaic')

GCI

 var districts = ee.FeatureCollection("users/geerootfold/gadm36_MLI_4")
  .filter(ee.Filter.eq('NAME_0', 'Mali'))



var year = 2019

    

// Add reducer output to the Features in the collection.
var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});

 
var modis1 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
    .filterDate('2019-05-01','2019-09-30')
var modis = modis1.max().clip(districts)

print(modis)
var B5 = modis.select('B5')
var B3 = modis.select('B3')
var gci = function(modis){
  var nir = modis.select('B5')
  var green = modis.select('B3')
  return nir.divide(green).subtract(1).rename('GCI')
}
// Add reducer output to the Features in the collection.
var modisGCI = gci(modis)
print('modisgci',modisGCI)



var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});

print(modisGCI)

var featureCollection = modisGCI.reduceRegions({
  collection: districts,
  reducer: reducers,
  scale: 30,
});

print(featureCollection)


var size=featureCollection.size()
var toList = featureCollection.toList(size)
print(toList, 'toList')

var idList = ee.List.sequence(0, size.subtract(1))


var newList = idList.map(function(x){
  
  var index = idList.get(x)
  var feat = toList.get(x)
  
  return ee.Feature(feat).set('year', 2019).copyProperties(feat)
})

print(newList, 'newList')

var MeansOfFeatures = ee.FeatureCollection(newList)


Export.table.toDrive({
  collection: MeansOfFeatures, 
  description: 'meanOfMaxGCI'+year,
  folder: 'GEE GCI2',
  // Explicitly specify your columns to exclude .geo
  fileFormat: 'CSV',
  selectors: ['NAME_4', 'year', 'mean','stdDev'] 
})


 var districts = ee.FeatureCollection("users/geerootfold/gadm36_MLI_4")
  .filter(ee.Filter.eq('NAME_0', 'Mali'))
var colour = {
 min: 0.0,
 max: 1.0,
 palette: [
 'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
 '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
 '012E01', '011D01', '011301'
 ],
};

Map.addLayer(modisGCI, colour, 'spatial mosaic')
    
// Add reducer output to the Features in the collection.
var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.stdDev(),
  sharedInputs: true
});

Link to District Shapefile: https://code.earthengine.google.com/?asset=users/geerootfold/gadm36_MLI_4

Link to NDVI code: https://code.earthengine.google.com/6fc9b38ab45af9e8dafc9a993a455a9e

Link to GCI code: https://code.earthengine.google.com/605e5de51cd8df0cffd2027f3b227c6c

Answer 1

For starters, you should try to learn how to use create function to reduce duplication in your scripts. You'd cut your 900 line of code GCI script down to 100 with that. It makes life a lot easier!

The main reason your GCI values are strange are because of clouds. You calculate the max values for green and NIR without masking out clouds, then calculate GCI based on these max values. That will more or less guarantee that if there are any clouds or haze, they will be included - these pixels have high green and NIR. I would suggest that you mask clouds, calculate the GCI for every Landsat image, then calculate the max GCI value:

var districts = ee.FeatureCollection("users/geerootfold/gadm36_MLI_4")
  .filter(ee.Filter.eq('NAME_0', 'Mali'))
  .select(['NAME_4'], ['region']) // Rename the property you're interested to something cleaner

Map.centerObject(districts)

var years = sequence(2015, 2019)
years.map(process)

function process(year) {
  var start = ee.Date.fromYMD(year, 5, 1)
  var end = ee.Date.fromYMD(year, 10, 1) // end is exclusive

  var maxGci = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterDate(start, end)
    .filterBounds(districts)
    .map(function (image) {
      var normalized = image
        .multiply(0.0000275).subtract(0.2)
      var cloudFree = bitwiseExtract(image.select('QA_PIXEL'), 0, 5).eq(0)
      return ee.Image()
        .expression('nir / green - 1', {
          green: normalized.select('SR_B3'),
          nir: normalized.select('SR_B5')
        })
        .updateMask(cloudFree)
        .rename('GCI')
    })
    .max()
    .clip(districts)
  
  var meanOfMaxGci = maxGci
    .reduceRegions({
      collection: districts,
      reducer: ee.Reducer.mean().combine({
        reducer2: ee.Reducer.stdDev(),
        sharedInputs: true
      }),
      scale: 3000
    })
    .map(function (feature) {
      return feature
        .set('year', year)
    })
  
  print(meanOfMaxGci)
  
  Export.table.toDrive({
    collection: meanOfMaxGci,
    description: 'meanOfMaxGCI' + year,
    folder: 'GEE GCI2',
    // Explicitly specify your columns to exclude .geo
    fileFormat: 'CSV',
    selectors: ['region', 'year', 'mean', 'stdDev']
  })
  
  
  var visParams = {
    min: 0,
    max: 6,
    palette: [
      'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
      '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
      '012E01', '011D01', '011301'
    ],
  }
  Map.addLayer(maxGci, visParams, 'max GCI ' + year)    
}


function bitwiseExtract(value, fromBit, toBit) {
  if (toBit === undefined)
    toBit = fromBit
  var maskSize = ee.Number(1).add(toBit).subtract(fromBit)
  var mask = ee.Number(1).leftShift(maskSize).subtract(1)
  return value.rightShift(fromBit).bitwiseAnd(mask)
}


function sequence(start, end) {
    return Array
      .apply(null, Array(end - start + 1))
      .map(function (_, i) {
        return i + start
      })
  }

https://code.earthengine.google.com/142e998c27800b51f9615753bff4e4bd