Calculate rainfall interception in GEE

Calculate how much rainfall is intercepted by the vegetation.

Interception is the amount of precipitation that does not reach the soil but is intercepted by the leaves and branches of plants. The formula below can be used to calculate the amount of interception.

a = emperical parameter (take a = 1 mm).
LAI = Leaf area Index.
p = Rainfall / number of rainy days.

The formula to calculate cc is shown below:

The steps below describe how to calculate the Interception

Step 1: Import these three Image Collections:
MCD15A3H.006 Terra+Aqua Leaf Area Index/FPAR 4-Day L4 Global 500 m
MYD13Q1.005 Vegetation Indices 16-Day Global 250m
CHIRPS: Climate Hazards Group InfraRed Precipitation with Station data (version 2.0 final)

Step 2: Import the fusion table for the Ca basin

// polygon for the Ca basin
var Ca = ee.FeatureCollection('ft:11uQewh4u7NFXneH99YPs6P9F1Au6wxBq959BZWt_','geometry');

Step 3: Set the date and time variables

// variables for start and end year
var startyear = 2003;
var endyear = 2004; 

// define start and end date
var startdate = ee.Date.fromYMD(startyear,1,1);
var enddate = ee.Date.fromYMD(endyear,12,1);

// create list for years
var years = ee.List.sequence(startyear,endyear);

// create list for months
var months = ee.List.sequence(1,12);

Step 4: Make a list with the average number of rainy days.

// raindays
var raindays = ee.List([7,10,13,9,14,9,13,16,15,11,7,8]);

Step 5: Filter the image collection based on location and time period.

// select P for time range
var P = chirps.filterDate(startdate, enddate)
  // Sort chronologically in descending order.
  .sort('system:time_start', false)
  .filterBounds(Ca);

// Leaf Area Index
var LAI = MCD15A.filterDate(startdate, enddate)
  // Sort chronologically in descending order.
  .sort('system:time_start', false)
  .filterBounds(Ca)
  .select('Lai');

// NDVI
var NDVI = MYD13Q1.filterDate(startdate, enddate)
  // Sort chronologically in descending order.
  .sort('system:time_start', false)
  .filterBounds(Ca)
  .select('NDVI');

Step 6: Calculate the monthly averages.

// calculate the LAI for each month
var MonthlyLAI =  ee.ImageCollection.fromImages(
  years.map(function (y) {
  return months.map(function(m){
  var w = LAI.filter(ee.Filter.calendarRange(y, y, 'year'))
           .filter(ee.Filter.calendarRange(m, m, 'month'))
           .mean();
  return w.set('year', y)
           .set('month', m)
           .set('date', ee.Date.fromYMD(y,m,1))
           .set('system:time_start',ee.Date.fromYMD(y,m,1))
           .set('index',m); 

});
}).flatten());

// calculate the p for each month
var MonthlyP =  ee.ImageCollection.fromImages(
  years.map(function (y) {
  return months.map(function(m){
  var w = P.filter(ee.Filter.calendarRange(y, y, 'year'))
           .filter(ee.Filter.calendarRange(m, m, 'month'))
           .sum();
  return  w.set('year', y)
           .set('month', m)
           .set('date', ee.Date.fromYMD(y,m,1))
           .set('system:time_start',ee.Date.fromYMD(y,m,1))
           .set('index',m);

});
}).flatten());

// calculate NDVI for each month
var MonthlyNDVI =  ee.ImageCollection.fromImages(
  years.map(function (y) {
  return months.map(function(m){
  var w = NDVI.filter(ee.Filter.calendarRange(y, y, 'year'))
           .filter(ee.Filter.calendarRange(m, m, 'month'))
           .mean();
  return  w.set('year', y)
           .set('month', m)
           .set('date', ee.Date.fromYMD(y,m,1))
           .set('system:time_start',ee.Date.fromYMD(y,m,1))
           .set('index',m) ;

});
}).flatten());

Step 7: Combine all data into one image collection

// Combine all data into one image collection
var Collection = MonthlyNDVI.combine(MonthlyLAI).combine(MonthlyP);

Step 8: Calculate the maximum and minimum NDVI

// calculate the maximum NDVI
var NDVImax = MonthlyNDVI.max();
// calculate the minimum NDVI
var NDVImin = MonthlyNDVI.min();

Step 8: Calculate the cc.

// Calculate cc
Collection = Collection.map(function(img){
  var myimg = img.expression("1 - ((ndvimax - ndvi)/(ndvimin - ndvi))",
  {
    ndvimax: NDVImax,
    ndvimin: NDVImin,
    ndvi: img.select("NDVI")
  });
  return img.addBands(myimg.rename('cc'));
});

Step 9: Calculate the interception.

// Calculate Interception
Collection = Collection.map(function(img){
  var nr = ee.Number(img.get("index")).subtract(1);
  var rainday = ee.Number(raindays.get(nr));
  var myimg = img.expression("1*lai*(1-(1/(1+(1+cc+(P/r))/(1*lai))))",
  {
    lai: img.select('Lai'),
    cc: img.select("cc"),
    P: img.select('precipitation'),
    r: rainday
  });
  return img.addBands(myimg.rename('I'));
});

Step 10: Add the layers to the map.

var viz = {min:0.0, max:10, palette:"000000,0000FF,f48042"};
Map.addLayer(Collection.select('I').mean().clip(Ca),viz,"iterception");
Map.centerObject(Ca,7);

See an example here.