// Google Earth Engine script for calculating EVI2 averages 
// over Sagehen & Independence catchments and subcatchments
// Version 1.0 
// Author: James Kirchner, ETH Zurich
// Copyright 2020 ETH Zurich and James Kirchner

// Public use permitted under GNU GPL3

// READ THIS CAREFULLY:
// ETH Zurich and James Kirchner make ABSOLUTELY NO WARRANTIES OF ANY KIND, including NO WARRANTIES, 
// expressed or implied, that this software is free of errors or is suitable for any particular 
// purpose.  Users are solely responsible for determining the suitability and
// reliability of this software for their own purposes.


// This example shows download of Terra satellite data.  
// For Aqua, unmask all references to Aqua and mask all references to Terra.


// list of polygons over which MODIS data will be extracted
var featureList = [
  ee.Feature(Independence01, {name: 'IND01'}),
  ee.Feature(LowerCulvert, {name: 'Lower'}),
  ee.Feature(MiddleCulvert, {name: 'Middle'}),
  ee.Feature(MainGauge, {name: 'Main'}),
  ee.Feature(KilnCreek, {name: 'Kiln'}),
  ee.Feature(SouthTrib1, {name: 'South1'}),
  ee.Feature(SouthTrib2, {name: 'South2'})];

var startDate = '2005-10-01'; var endDate = '2008-10-01'; 
//note: in google earth engine, start dates are INCLUSIVE but end dates are EXCLUSIVE


var patches = ee.FeatureCollection(featureList);  //cast polygons as a feature collection


//load MODIS daily surface reflectance
//MOD09GA is Terra (morning view)
//MYD09GA is Aqua (afternoon view)
var dataset = ee.ImageCollection('MODIS/006/MOD09GA')
//var dataset = ee.ImageCollection('MODIS/006/MYD09GA')
                  .filterBounds(domain1)
                  .filterDate(startDate, endDate);



// helper function to extract the QA bits
function getQABits(image, start, end, newName) {
    // Compute the bits we need to extract.
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    // Return a single band image of the extracted QA bits, giving the band
    // a new name.
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
}



// A function to mask out problem pixels.
function maskQuality(image) {
 // Select the QA band.

 var QA = image.select('QC_500m');
 var bit0 = getQABits(QA, 0, 1, 'internal_quality_flag'); // = 0 indicates ideal quality all bands


 var internalQuality = bit0;

 // Return an image masking out problem pixels.
 return image.updateMask(internalQuality.eq(0));
}

// create clean image collection
var clean = dataset.map(maskQuality);




// calculate EVI   NOTE THIS IS EVI2!
var addEVI = function(image){
  
  var red = image.select('sur_refl_b01');
  var nir = image.select('sur_refl_b02');
  
  var evi = image.expression(   //NOTE THIS IS EVI2
    '2.5 * ((NIR - RED) / (NIR + 2.4 * RED + 16000.0))', {  //note this is 16000 because we are working with the raw bands, which go from -100 to 16000
      'NIR':nir, 'RED':red
    }).rename('EVI');
  
  //mask anything that started with "fill values" of red or nir
  var good = red.gt(-100).and(nir.gt(-100));
  evi = evi.updateMask(good);
  
  return image.addBands(evi);
};


// add EVI to image collection
var cleanEVI = clean.map(addEVI);





//Function to extract values from image collection based on point file and export as a table 
var fill = function(img, ini) {
  // type cast
  var inift = ee.FeatureCollection(ini);
  // gets the values for the points in the current img
  // gets the date of the img
  var date = img.date();
  
  var first = date; first = first.advance(-3, "day");  //start date is inclusive
  var last = date; last = last.advance(4, "day");  //end date is exclusive
  
  var median = evi.filterDate(first, last).median();  //median composite of 7-day stack of evi
 
  var ft2 = median.reduceRegions({collection: patches, reducer: ee.Reducer.mean(), scale: 100});
  // writes the date in each feature
  var ft3 = ft2.map(function(f){return f.set("date", date.format())});
  // merges the FeatureCollections
  return inift.merge(ft3);
};  //end of function fill




//Function to calculate pixel completeness of each image in image collection and export as a table 
var tally = function(img, ini) {
  // type cast
  var inift = ee.FeatureCollection(ini);
  // create binary layer that will equal 1 for all non-missing values
  var nonmissing = img.lt(1).or(img.gt(0)); //note that this is always true....
  // now take sum of these 0 and 1 values
  var ft2 = nonmissing.reduceRegions({collection: patches, reducer: ee.Reducer.sum(), scale: 100});
  // gets the date of the img
  var date = img.date().format();
  // writes the date in each feature
  var ft3 = ft2.map(function(f){return f.set("date", date)});
  // merges the FeatureCollections
  return inift.merge(ft3);
};  //end of function tally





// Empty Collection to fill
var ft = ee.FeatureCollection(ee.List([]));

var evi = cleanEVI.select('EVI');
// Iterates over the ImageCollection
var evi_means = ee.FeatureCollection(evi.iterate(fill, ft));



// export patch means
Export.table.toDrive({
  collection: evi_means,
  description: 'evi2Means_TerraComposite_includeclouds',
//  description: 'evi2Means_AquaComposite_includeclouds',
  folder: 'Sagehen',
  fileFormat: 'CSV',
  selectors: 'name, date, mean'
});





// Empty Collection to fill
var ft2 = ee.FeatureCollection(ee.List([]));

var counts = ee.FeatureCollection(evi.iterate(tally, ft2));


// export patch counts
Export.table.toDrive({
  collection: counts,
 description: 'evi2_patch_counts_TerraComposite_includeclouds',
//  description: 'evi2_patch_counts_AquaComposite_includeclouds',
  folder: 'Sagehen',
  fileFormat: 'CSV',
  selectors: 'name, date, sum'
});