How to calculate VCI index in GEE?

Question

I would like to calculate the VCI index from MODIS/006/MCD43A4 NDVI over the months July and August for the period 2002-2017. I have the following code, but I do not know how to continue.

var imageCollection = ee.ImageCollection("MODIS/006/MCD43A4")
 .filterBounds(geometry)
    .map(function(image){return image.clip(geometry)});
   var modNDVI = imageCollection.select("Nadir_Reflectance_Band2","Nadir_Reflectance_Band1","BRDF_Albedo_Band_Mandatory_Quality_Band1","BRDF_Albedo_Band_Mandatory_Quality_Band2");

var qas = function(image){ 
  var mask1 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band1").eq(0);
  var mask2 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band2").eq(0);
  return image.updateMask(mask1).updateMask(mask2);
};

var merged = modNDVI.map(qas);

var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['Nadir_Reflectance_Band2', 'Nadir_Reflectance_Band1']).rename('NDVI');
  return image.addBands(ndvi);
};

var ndvi = merged.map(addNDVI);

var NDVI=ndvi.select('NDVI')

var reclassified = NDVI.map(function(img){
      return img.updateMask(img.gt(0.1))
    })

print(reclassified)

I am aware of the VCI form, however, I am not sure how to calculate it in GEE.

I tried to filter the months as folow, but I am not sure if I am in the right path.

var mod_JA = reclassified.filter(ee.Filter.calendarRange(1, 6, 'month'))
  .filter(ee.Filter.calendarRange(2002, 2017, 'year'))
  .map(function(img) {
    return img.set('year', img.date().get('year'));
  });
var mod_JAjoin = ee.Join.saveAll('same_year').apply({
  primary: mod_JA.distinct('year'),
  secondary: mod_JA,
  condition: ee.Filter.equals({leftField: 'year', rightField: 'year'})
});

Any help?

Noel Gorelick · Accepted Answer · 2021-07-30 12:12:38Z

Though you could eventually get there with a join, you don't strictly need it. You just need to figure out which historical values you want, compute the min/max from them and then do the math.

It's not clear if you want to compute VCI for each day or each month, and it's also not clear if you want to do that for one year, of for all years, but some combination/modification to this code which computes it daily for each day in July and August of 2017, using the 15 years before that to compute min/max, should get you there.

If you want monthly, then make monthly composites first by mapping over months and take a mean of all the images in that month first (don't forget to set the system:time_start for the composites).

var vci = function(image) {
  var date = image.date()
  var doy = date.getRelative('day', 'year')

  // Get a collection of historic values for this DOY.
  var history = reclassified
      // filter down to everything before this image.
      .filterDate('2000-01-01', ee.Date.fromYMD(date.get('year').subtract(1), 1, 1))
      // and the same DOY.
      .filter(ee.Filter.calendarRange(doy, doy.add(1), 'day_of_year'))
  var min = history.min()
  var max = history.max()
  var vci = image.subtract(min).divide(max.subtract(min))
  return vci.rename('vci')
}

// Compute the VCI for each day in July/Aug of 2017
var vci = reclassified.filterDate('2017-07-01', '2017-08-30').map(function(img) {
  return img.addBands(vci(img))
})
Map.addLayer(vci.first())

MrXsquared · Accepted Answer · 2023-02-16 10:50:59Z

You might want to look at this code...

/////////////////////////////  MODIS VCI Time Series /////////////////////////////////////////


// ---------------------------------------------------------------------------------------------
///// Preprocessing Functions /////

// Compute MODIS QA bits
var getQABits = function(image) { // Compute  bits to extract
    var pattern = 0;
    var start = 0;
    var end = 1;
    var newName = "cloudmask";
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i); 
    }
    return image.select([0], [newName]) // Return single band image of the extracted QA bits
                .bitwiseAnd(pattern)
                .rightShift(start);
};

// Cloud Mask Function
var maskClouds = function(image) {
  var QA = image.select("DetailedQA"); // Select the QA band
  var internalCloud = getQABits(QA); // Get the internal_cloud_algorithm_flag bit
  return image.updateMask(internalCloud.eq(0)); // Return image masking out cloudy areas
};

// Scale MODIS Vegetation Indices by a factor of 0.0001
var scaleMOD = function(image2) {
        return  image2
          .divide(10000) //divide image values by 10000
          .float() // convert to float
          .set("system:time_start", image2.get("system:time_start")); // keep time info
};


// Convert MODIS Land Cover Dynamics (time info) days since epoch to doy
var doy = function(img) {
  // get days since epoch for first day of year from properties
  var firstDayOfYear = ee.Number(img
    .get("system:time_start"))
    .divide(24*60*60*1000);
  var doy = img.subtract(firstDayOfYear).add(1);
  return doy
    .set("system:time_start", img.get("system:time_start"));
};

// ---------------------------------------------------------------------------------------------


// ---------------------------------------------------------------------------------------------
///// Dataset Imports & Preprocessing /////

// ROI als Layer für Thüringen
// Study Area - Region of Interest Thueringen

print (roi)
Map.addLayer(roi)
Map.centerObject(roi, 7);


// MODIS Image Collections: Einmal eine Vegetations Kollektion mit 250 Auflösung
// MODIS Image Collections: Einmal mit Wasser Vegetation mit 250m Auflösung
// MODIS Vegetation Data
var terra_vegetation = ee.ImageCollection("MODIS/006/MOD13Q1") // provides a Vegetation Index (VI) value at a per pixel basis
  .filterBounds(roi)
  .map(maskClouds)
  .map(scaleMOD);
var aqua_vegetation = ee.ImageCollection("MODIS/006/MYD13Q1") // provides a Vegetation Index (VI) value at a per pixel basis
  .filterBounds(roi)
  .map(maskClouds)
  .map(scaleMOD);
var modis_vegetation = terra_vegetation.merge(aqua_vegetation); 
print("MODIS Vegetation Indices", modis_vegetation);


// ---------------------------------------------------------------------------------------------
///// MODIS Monthly NDVI Anomaly Calculation I /////

// Einen Zeitraum festelgen, in welchem die Untersuchung stattfinden soll
// Define first & last year of interest + search months
var startyear = 2019; 
var endyear = 2020;
var startmonth = 1; 
var endmonth = 12; 

// Zeitrahmen berechnen und eine Sequenz erstellen mit MOnaten und Jahren -- keine Tage mehr
// Compute beginning & end of study period + sequence of months & years
var startdate = ee.Date.fromYMD(startyear, startmonth, 1);
var enddate = ee.Date.fromYMD(endyear , endmonth, 30);
var years = ee.List.sequence(startyear, endyear);
var months = ee.List.sequence(startmonth,endmonth);

// Nur NDVI aus der MODIS Kollektion nehmen
// filter MODIS vegetation indices collection
var ndviCollection = modis_vegetation
  .filterDate(startdate, enddate)
  .select("NDVI");

// ---------------------------------------------------------------------------------------------
///// MODIS Monthly NDVI Anomaly Calculation I /////


// Monatliche NDVI Kompisition erstellen von April bis September für jedes Jahr
// Create NDVI composite for every month
var monthlyNDVI =  ee.ImageCollection.fromImages(
  years.map(function (y) { 
    return months.map(function(m) {
      var monthly = ndviCollection
        .filter(ee.Filter.calendarRange(y, y, "year"))
        .filter(ee.Filter.calendarRange(m, m, "month"))
        .mean(); 
      return monthly
        .set("year", y) 
        .set("month", m) 
        .set("system:time_start", ee.Date.fromYMD(y, m, 1));}); })
  .flatten());
print (monthlyNDVI, 'monthly NDVI')





// Für jeden Monat werden die NDVI Maximal Werte ausgerechnet --> somit bekommen wir eine
// MAX Value Composition
// Calculate mean values for each month over all years
var MonthlyMAX =  ee.ImageCollection.fromImages(months
  .map(function (m) {
    var maxNDVI = monthlyNDVI
      .filter(ee.Filter.eq("month", m))
      .select("NDVI")
      .reduce(ee.Reducer.percentile({percentiles: [90]}))
      .rename("max_NDVI");
  return maxNDVI
    .set("month", m);})
  .flatten());
print (MonthlyMAX, 'MonthlyMAX');
Map.addLayer (MonthlyMAX.first().select('max_NDVI').clip(roi),  {min:0, max:1,  'palette': ['red','yellow', 'green']}, 'MonthlyMAX')

// Für jeden Monat werden die NDVI Minimal Werte ausgerechnet --> somit bekommen wir eine
// MIN Value Compisition
var MonthlyMIN =  ee.ImageCollection.fromImages(months
  .map(function (m) {
    var minNDVI = monthlyNDVI
      .filter(ee.Filter.eq("month", m))
      .select("NDVI")
      .reduce(ee.Reducer.percentile({percentiles: [10]}))
      .rename("min_NDVI");
  return minNDVI
    .set("month", m);})
  .flatten());
print (MonthlyMIN, 'MonthlyMIN');
Map.addLayer (MonthlyMIN.first().select('min_NDVI').clip(roi), {min:0, max:1,  'palette': ['red','yellow', 'green']}, 'MonthlyMIN')

var MonthlyMIN_img = MonthlyMIN.mean();
var MonthlyMAX_img = MonthlyMAX.mean();

/////////////////////////////////////////////////////// Calculation of VCI //////////////////////////////////////////////////////

var Mndvi2 = monthlyNDVI.map(
        function(img) {
        var vci = img.select('NDVI')
                          .subtract(MonthlyMIN_img).divide((MonthlyMAX_img.subtract(MonthlyMIN_img))).multiply(100)
                          .rename('vci');
                          
        return img.addBands(vci);
  }
);

var VCI = Mndvi2.select('vci')
print(VCI, 'VCI')


print(ui.Chart.image.series(VCI , roi , ee.Reducer.mean(), 5000));
Map.addLayer(VCI.first())

/*


//Create a list to calulate the Index for each image related to each month
var vci_list = [];

for(var j = startyear; j <=endyear; j++){
  for(var i = startmonth; i <=endmonth; i++) {
  
    // hier wird dann jeweils der value für i, j hinten an die Liste (push) gesetzt. diese list enthält dann die werte
    // hier sollen der vci mit min und max des jeweiligen Monat i für die 20 Jahre und in eine Liste gepusht werden.
    // NDVI liegen als ImageCollection vor (Bild)
    // VCI Formel VCI = 100 * (NDVI_monat - NDVI_monat_min) / (NDVI_monat_max - NDVI_monat_min)
    // also VCI = 100 * (monthlyNDVI - monthlyMIN) / (monthlyMAX - monthlyMIN)

  var current_vci = 100 * (monthlyNDVI[j][i] - MonthlyMIN[i]) / (MonthlyMAX[i] - MonthlyMIN[i]);
 // var current_vci = 100 * ('NDVI'[j][i] - 'min_NDVI'[i]) / ('max_NDVI'[i] - 'min_NDVI'[i]);
    vci_list.push(current_vci); 
  }
}

print(vci_list);


/////////////////////////////////// second attempt //////////////////////////////


var modis_list = [];
var modis_filtered_flattened = 
    ee.ImageCollection(ee.FeatureCollection(modis_list).flatten());
for(var i = startmonth; i <=endmonth; i++) {
  for(var j = startyear; j <=endyear; j++){
    modis_list.push(monthlyNDVI.filter(ee.Filter.calendarRange(i, i, 'month'))
                          .filter(ee.Filter.calendarRange(j, j, 'year')));
  }
}
print(modis_list);

//VCI calculation
var vci = modis_filtered_flattened.map(function(img){
 var id = img.id();
 var min = monthlyMIN;
 var max = monthlyMAX;
 return img.expression(
   "(NDVI-min)/(max-min)",{
     "NDVI" : img,
     "max" : ee.Number(max),
     "min" : ee.Number(min)
   });
});
print(vci)





var VCI = monthlyNDVI.map(function(image){
    return month.addBands(month.expression(
     "100*(NDVI-min)/(max-min)",{
       "NDVI": month
       "max": ee.Number(MonthlyMAX)
       "min": ee.Number(MonthlyMIN)
       ))
    });

print (VCI)
*/

/*
//VCI calculation (Stackexchange)
var vci = modis_list.map(function(img){
 var id = img.id();
 return img.expression(
   "(NDVI-min)/(max-min)",{
     "NDVI" : img,
     "max" : ee.Number(max),
     "min" : ee.Number(min)
   }).copyProperties(img,['system:time_start','system:time_end']);
});
print(vci)
*/

/*
var VCI2 = monthlyNDVI.map(function(image) {
  var max = month.select('NDVI').reduce(ee.Reducer.percentile({percentiles: [90]}));
  var min = month.select('NDVI').reduce(ee.Reducer.percentile({percentiles: [10]}));
  return month.addBands(month.expression(
   "100*(NDVI-min)/(max-min)",{
     "NDVI" : month,
     "max" : ee.Number(max),
     "min" : ee.Number(min)
   }).rename('VCI')).copyProperties(month,['system:time_start','system:time_end']);
});


/*
//@Jojo: modis list ist ein placeholder ([] bedeutet leeres array bzw list)
var modis_list = [];
// @Jojo: Outer loop: Monate  - inner Loop Jahre
for(var i = startmonth; i <=endmonth; i++) {
  for(var j = startyear; j <=endyear; j++){
    // @jojo hier wird dann jeweils der value für i, j hinten an die Liste (push) gesetzt. diese list enthält dann die werte
   var vci = img1.select('NDVI')
                          .subtract(MonthlyMIN).divide((MonthlyMAX.subtract(MonthlyMIN))).multiply(100)
                          .rename('vci');
    modis_list.push(monthlyNDVI
                     .first()); 
                     }
}
// @Jojo from Internet: .first() to get the first image of the collection. (This will usually be the oldest matching image.)

print(modis_list);

*/


/*
monthlyNDVI = monthlyNDVI.map(
        function(img1) {
        var vci = img1.select('NDVI')
                          .subtract(MonthlyMIN).divide((MonthlyMAX.subtract(MonthlyMIN))).multiply(100)
                          .rename('vci');
                          

    return img1.addBands([vci, Nratio]);
  }
);

var VCI = monthlyNDVI.select('vci')
var Nrat = monthlyNDVI.select('Nratio')

print (VCI);
*/



/*

*/
/*
var VCI = monthlyNDVI.map(function(image) {
  var equation = image.expression("100*(NDVI-min)/(max-min)", { 
      NDVI: monthlyNDVI.select("NDVI"),
      min: MonthlyMIN.select("min_NDVI"),
      max: MonthlyMAX.select("max_NDVI"),
    }).rename("VCI");
  return equation.copyProperties(image);
}); 

print (VCI2, 'VCI')
*/

/*
var vci = monthlyNDVI.map(function(img){
 var id = img.id();
 return img.expression(
   "100*(NDVI-min)/(max-min)",{
     "NDVI" : img,
     "max" : ee.Number(max_NDVI),
     "min" : ee.Number(min_NDVI)
   }).copyProperties(img,['system:time_start','system:time_end']);
});

print (vci);
*/


/*
var VCI = ee.ImageCollection.fromImages(months
  .map(function (m) { 
    var VCI2 = monthlyNDVI
      .filter(ee.Filter.eq("month", m))
      .select("NDVI")
      .reduce(ee.Reducer.stdDev())
      .rename("sdv_ND");
  return VCI2
    .set("month", m);})
  .flatten());




var VCI2 = monthlyNDVI.map(function(image) {
  var max = month.select('NDVI').reduce(ee.Reducer.percentile({percentiles: [90]}));
  var min = month.select('NDVI').reduce(ee.Reducer.percentile({percentiles: [10]}));
  return month.addBands(month.expression(
   "100*(NDVI-min)/(max-min)",{
     "NDVI" : month,
     "max" : ee.Number(max),
     "min" : ee.Number(min)
   }).rename('VCI')).copyProperties(month,['system:time_start','system:time_end']);
});


// Calculate standard deviation for each month over all years 
var sdvMonthlyND = ee.ImageCollection.fromImages(months
  .map(function (m) { 
    var sdvND = monthlyNDVI
      .filter(ee.Filter.eq("month", m))
      .select("NDVI")
      .reduce(ee.Reducer.stdDev())
      .rename("sdv_ND");
  return sdvND
    .set("month", m);})
  .flatten());

// Helper functions for joining 
var filtereq = ee.Filter.equals({
  leftField: "month",
  rightField: "month",
});

var join = ee.Join.saveFirst({
    matchKey: "match",
});

// Join monthly NDVI collection with mean
var ND_Join1 = ee.ImageCollection(join.apply(monthlyNDVI,meanMonthlyND,filtereq))
  .map(function(image) {
    return image.addBands(image.get("match"));
  });

// Join previously joined NDVI/mean with standard deviation
var ND_Join2 = ee.ImageCollection(join.apply(ND_Join1,sdvMonthlyND,filtereq))
  .map(function(image) {
    return image.addBands(image.get("match"));
  });

// Compute standardized anomaly

var standardizedAnomaly = function(image3) {
  var anomaly = image3.expression("(ndvi - average)/sd", { 
      ndvi: image3.select("NDVI"),
      average: image3.select("mean_ND"),
      sd: image3.select("sdv_ND"),
    }).rename("ND_anomaly");
  return anomaly.copyProperties(image3).set('system:time_start',image3.get("system:time_start"));
}; 
var ND_Anomaly = ND_Join2.map(standardizedAnomaly);
print("Standardized NDVI Anomalies", ND_Anomaly); 
// ---------------------------------------------------------------------------------------------

Map.addLayer(ND_Anomaly.first().select('ND_anomaly').clip(roi), {min:-3, max:3,  'palette': ['red','yellow', 'green']}, 'nd anomaly')


*/

Stack Exchange Network

How to calculate VCI index in GEE?

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
google-earth-engine
modis
vci
or ask your own question.

Hot Network Questions

How to calculate VCI index in GEE?

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

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

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
google-earth-engine
modis
vci
or ask your own question.