Google Earth Engine: calculate percentile rank in ee.ImageCollection

Question

I am trying to calculate the Vegetation Productivity Index sensu copernicus in Google Earth Engine, which requires estimating the percentile rank. I haven't found any lead so far on how to do so, does anyone know how to do it or recommend a way to explore?

Given the example below, how would you calculate the percentile rank of the NDVI values in img compared to those in ref? What I am expecting is a value from 0-1, where the value at each pixel is the rank compared to all historical values for the same pixel.

ref = ee.ImageCollection('MODIS/006/MYD13Q1').filterDate('2006-01-01', '2006-12-31').select('NDVI') 
img = ee.ImageCollection('MODIS/006/MYD13Q1').filterDate('2007-01-01', '2007-01-10').select('NDVI').first()

Answer 1

By definition, Percentile Rank (not to be confused with Percentile function or Percentile) is expressed as percentage so, if you want values from 0-1 you can only divide by 100. For calculating Percentile Rank (PR), the formula is as follows; where CF—the cumulative frequency—is the count of all scores less than or equal to the score of interest, F is the frequency for the score of interest, and N is the number of scores in the distribution.

In your case, the score of interest is img and all scores less than or equal to the score of interest are in ref. On the other hand, MODIS/006/MYD13Q1 products must be divided by 10000 for obtaining corresponding NDVI values. I'm going to assume an arbitrary point in USA (-106.59086720561201, 38.81137001051776), however, when it is used this point with 'reduceRegion' method, some nulls values are produced with your dates range for NDVI_ref and they need to be eliminated. Following code determines PR based in above assumptions and formula, calculating its respective CF, F and N values for your date ranges.

var pt = ee.Geometry
  .Point([-106.59086720561201, 38.81137001051776]);

function divideImages(image) {
  return image.divide(10000);
}

var ref = ee.ImageCollection('MODIS/006/MYD13Q1')
  .filterDate('2006-01-01', '2006-12-31')
  .select('NDVI')
  .map(divideImages); 

var img = ee.ImageCollection('MODIS/006/MYD13Q1')
  .filterDate('2007-01-01', '2007-01-10')
  .select('NDVI')
  .first();
  
print(ref);

img = ee.ImageCollection(img.divide(10000));

Map.addLayer(img);
Map.addLayer(pt);
Map.centerObject(pt);

var getNDVIvalues = function(image) {

  // Reducing region and getting value
  var ndvi_value = ee.Image(image)
    .reduceRegion(ee.Reducer.first(), pt)
    .get('NDVI');

  return ndvi_value;

};

var size1 = ref.size();

print("size list (including nulls)", size1);

var NDVI_ref = ref.toList(size1).map(getNDVIvalues);

var N = NDVI_ref.reduce(ee.Reducer.count());

print("size list (without nulls)", N);

print("NDVI ref", NDVI_ref);

var size2 = img.size();

var NDVI_img = img.toList(size2).map(getNDVIvalues).get(0);

print("NDVI img", NDVI_img);

var NDVI_ref = NDVI_ref.filter(ee.Filter.neq('item', null));

var CF_values = NDVI_ref.map(function count (ele){
  
  return ee.Number(ele).lte(NDVI_img).multiply(ele);
  
}).removeAll([0]);

print("CF_values", CF_values);

var CF = CF_values.size();

print("CF", CF);

var F = ee.Algorithms.If(NDVI_ref.frequency(NDVI_img).eq(0), 1, NDVI_ref.frequency(NDVI_img));

print("F", F);

//PR = percentile rank
var PR = ((CF.add(ee.Number(0.5).multiply(F))).divide(N)).multiply(100);

print("PR", PR);

After running above code in GEE code editor, I got result of following image:

where PR of NDVI_img (0.4865) is 40.476190476190474. Complete list of NDVI_ref values are:

0.4369, 0.475, 0.3211, 0.3146, 0.204, 0.4482, 0.2881, 0.4279, 0.5854, 0.5794, null, 0.6037, 0.0748, 0.6714, 0.655, 0.6442, 0.6564, 0.7149, 0.7389, 0.709, 0.4984, 0.5799, null

They can be used for corroborating manually above formula for PR.

Answer 2

I suggest doing a reduceRegion or reduceNeighborhood using ee.Reducer.percentile as your reducer.

Your code would look something like this:

var ndviPercentile = img.reduceNeighborhood({
  reducer: ee.Reducer.percentile([90]),
  geometry: inputGeometry,
  scale: /*depends on the size of your input geometry*/,
  maxPixels: 1e13
});