// Google Earth Engine script for calculating NDSI 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-09-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
 
 QA = image.select('state_1km');
 var bit10 = getQABits(QA, 10, 10, 'internal_quality_flag'); //0 indicates no clouds
 var bit13 = getQABits(QA, 13, 13, 'internal_quality_flag'); //0 indicates no adjacent clouds


 var internalQuality = bit0.add(bit10.add(bit13));

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

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




// calculate ndsi  
var addNDSI = function(image){
  
  var green = image.select('sur_refl_b04');
  var swir = image.select('sur_refl_b06');
  
  var ndsi = image.normalizedDifference(['sur_refl_b04', 'sur_refl_b06']).rename('NDSI');

  //mask anything that started with "fill values" of green or swir
  var good = green.gt(-100).and(swir.gt(-100));
  ndsi = ndsi.updateMask(good);
  
  return image.addBands(ndsi);
};


// add NDSI to image collection
var cleanNDSI = clean.map(addNDSI);





//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");
  var last = date; last = last.advance(4, "day"); //because end dates are EXCLUSIVE
  
  var median = ndsi.filterDate(first, last).median();  //median composite of 7-day stack of ndsi
 
  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 ndsi = cleanNDSI.select('NDSI');
// Iterates over the ImageCollection
var ndsi_means = ee.FeatureCollection(ndsi.iterate(fill, ft));



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





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

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


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