1

I'm relatively new to scripting in the Google Earth Engine and have an issue with calculating the Standard Vegetation Index (SVI) for a MODIS image collection. I want to write a function, which uses ee.Reducer.mean() and ee.Reducer.stdDev() of an image to calculate the SVI and map it over the image collection. However, GEE gives me the following Error-Code:

Error in map(ID=XXX_XX_XX): Image.constant: Invalid Image.constant type.

Where "ID=XXXX_XX_XX" seems to be a random date. Has anybody an idea how to solve that problem?

//Add ROI
var ROI = ee.Geometry.Polygon(
  [[-2.2535247523754443,5.881894109545312], 
  [-0.31993100237544425,5.881894109545312], 
  [-0.31993100237544425,9.669268383910472], 
  [-2.2535247523754443,5.881894109545312]]);

//Add MODIS ImageCollection
var MODIS = ee.ImageCollection("MODIS/006/MOD13Q1")
  .select(['EVI'])
  .filterDate('2000-01-01', '2020-01-01')
  .filterBounds(ROI)
  .map(function(image){return image.clip(ROI)});

//Function to Map SVI over image collection
var CalculateSVI = MODIS.map(function(SVI){
  var SelectImage = SVI.select(['EVI'], ['SVI']);

//Calculate Mean of Image  
  var Mean = ee.Number(SelectImage.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: ROI,
  bestEffort: true,
  }));

//Calculate Stdv of Image  
  var Stdv = ee.Number(SelectImage.reduceRegion({
  reducer: ee.Reducer.stdDev(),
  geometry: ROI,
  bestEffort: true,
  }));

//Calculate SVI
  var CalcSVI =  SelectImage.expression(
  ((SelectImage.subtract(Mean)).divide(Stdv))
  );

  return CalcSVI.addBands(SVI);
});
print(CalculateSVI);
1

reduceRegion() returns a ee.Dictionary(). Since you selected one band, you will always have one output value. You can then get that value from the dictionary using .values().get(0). Also, you don't need to use the expression() function.

//Calculate Mean of Image  
  var Mean = ee.Number(SelectImage.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: ROI,
  bestEffort: true,
  }).values().get(0));

//Calculate Stdv of Image  
  var Stdv = ee.Number(SelectImage.reduceRegion({
  reducer: ee.Reducer.stdDev(),
  geometry: ROI,
  bestEffort: true,
  }).values().get(0));

//Calculate SVI
  var CalcSVI = (SelectImage.subtract(Mean)).divide(Stdv)

see some other suggestion in the code. Note that this is just the solution of the code error you got. I have no idea if the actual SVI values are correctly calculated.

2
  • Thank you, works perfect! Unfortunately, the long timeline is needed. Is there a possibility to prevent a time out?
    – Torben
    May 13 '20 at 12:59
  • export the images or mean/stdv values as a table
    – Kuik
    May 13 '20 at 13:22
0

Apart from the problems @kuik points out, I'm wondering if you're calculating SVI correctly. I've never done it myself, but look at this article:

The SVI is a z-score deviation from the mean in units of the standard deviation, calculated from the NDVI or EVI values for each pixel location of a composite period for each year during a given reference period. The equation below shows the general calculation of the SVI

It seems like you should calculate your mean and standard deviation by pixel over time, not over your region. Here's my take on this, which might not be correct, but reflects my understanding of the above article:

var region = ee.Geometry.Polygon([
  [-2.2535247523754443, 5.881894109545312],
  [-0.31993100237544425, 5.881894109545312],
  [-0.31993100237544425, 9.669268383910472],
  [-2.2535247523754443, 5.881894109545312]
])

var startDate = ee.Date('2000-01-01')
var endDate = ee.Date('2020-01-01')

var collection = ee.ImageCollection("MODIS/006/MOD13Q1")
  .filterDate(startDate, endDate)
  .filterBounds(region)
  .map(maskClouds)
  .select(['EVI'])

var deltaDays = 7
var statsCollection = ee.ImageCollection(
  ee.List.sequence(0, 365, deltaDays).map(calculateStats)
)

var sviCollection = collection.map(toSVI)  
print(sviCollection)

Map.centerObject(region, 9)
var first = sviCollection
  .first()
  .clip(region)

var sviVis = {
    min: -2,
    max: 2,
    palette: ["ffffff", "ce7e45", "df923d", "f1b555", "fcd163", "99b718", "74a901", "66a000", "529400", "3e8601", "207401", "056201", "004c00", "023b01", "012e01", "011d01", "011301"]
  }  
Map.addLayer(first, sviVis, 'first')


function calculateStats(startDOY) {
  var endDOY = ee.Number(startDOY).add(deltaDays)
  return collection
    .filter(ee.Filter.calendarRange(startDOY, endDOY, 'day_of_year'))
    .reduce(ee.Reducer.stdDev().combine(ee.Reducer.mean(), null, true))
    .set('startDOY', startDOY)
    .set('endDOY', endDOY)  
}

function toSVI(image) {
  var doy = image.date().getRelative('day', 'year')
  var stats = statsCollection
    .filter(ee.Filter.lte('startDOY', doy))
    .filter(ee.Filter.gt('endDOY', doy))
    .first()
  return image
    .expression(
      '(evi - mean) / stdDev', {
        evi: image, 
        mean: stats.select('EVI_mean'), 
        stdDev: stats.select('EVI_stdDev')
    })
    .rename('svi')
    .copyProperties(image, image.propertyNames())
}

function maskClouds(image) {
  var qa = image.select('DetailedQA')

  // https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD13Q1
  var quality = bitwiseExtract(qa, 0, 1)
  var usefulness = bitwiseExtract(qa, 2, 5)
  var aerosolQuantity = bitwiseExtract(qa, 6, 7)
  var adjacentCloudDetected = bitwiseExtract(qa, 8)
  var atmosphereBrdfCorrection = bitwiseExtract(qa, 9)
  var mixedClouds = bitwiseExtract(qa, 10)
  var landWaterMask = bitwiseExtract(qa, 11, 13)
  var possibleSnowIce = bitwiseExtract(qa, 14)
  var possibleShadow = bitwiseExtract(qa, 15)
  var mask = quality.lte(1) // VI produced, but check other QA
    .and(usefulness.lte(2)) // Lower quality
    .and(aerosolQuantity.lte(2)) // Intermediate
    .and(adjacentCloudDetected.eq(0)) // No
    .and(mixedClouds.eq(0)) // No
    .and(possibleSnowIce.eq(0)) // No
    .and(possibleShadow.eq(0)) // No
  return image.updateMask(mask)
}

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

https://code.earthengine.google.com/1e552b88ad6f27f12acf58adcc01b8ba

8
  • Thank you. Yes, You have a point there. I was not sure whether to calculate it over time or region. I have a call with specialists tomorrow and will ask which way is the correct approach. I will let you know what they say about that. Thank you for giving a solution, in case I chose the wrong formula.
    – Torben
    May 13 '20 at 13:10
  • I do have a question regarding your masking. I had no idea how to handle the bit encoding of the "SummaryQA", so I found another approach. My mask affects mainly the southern part of Ghana, which seemed logical for me, so I did not investigate that further. But yours mainly masks out the northern regions. Here is my scrip so far: code.earthengine.google.com/… Can you tell me, whether I made a mistake with the mask? Thank you!
    – Torben
    May 13 '20 at 14:06
  • Not sure, whether I shared it correctly
    – Torben
    May 13 '20 at 14:19
  • My bad on the cloud masking. I've updated it properly now. Your masking code is correct, but it will only work for the first bits, like you have in the SummaryQA band. You include 1 - Marginal data, useful but look at detailed QA for more information. You might want to follow that advice and use the DetailedQA band. Then it's useful with a utility method like my bitwiseExtract(), where you just provide the bit range you want the value for. May 13 '20 at 14:41

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.