# Google Earth Engine: calculate percentile rank in ee.ImageCollection

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()
``````

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.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

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.

• Thanks for this comprehensive answer. I've been trying to adapt this to work on full images (ee.Arrays I guess?) rather than single points but I got a bit lost. Any idea on how to achieve this?
– e5k
Commented Mar 31, 2021 at 2:39
• MODIS/006/MYD13Q1 products contain global scale images (too big and too much masked pixels). So, I suggest to use your ROI for calculating the Vegetation Productivity Index placing a feature collection grid type on it. For each polygon of the grid, it would be determined its PR (mapping this feature collection; not by using an ee.Array object). Commented Apr 1, 2021 at 2:19

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
});
``````
• Percentile rank is not percentile function or percentile. You can read correct definition here: en.wikipedia.org/wiki/Percentile_rank . On the other hand, he needs to calculate PR of NDVI_img in relation to NDVI_ref values; not percentile 90. Commented Mar 31, 2021 at 0:02
• Ah yes, you're right, thanks! Commented Mar 31, 2021 at 18:12