BivariateTool.txt

//--------------------------------------------------------------------
//                        Description
//--------------------------------------------------------------------
// Master Code for NDTI/SST bivariate map and time series charts
// OSA Water Resources Project (III)
// NASA Develop, Fall 2018
// Bivariate related code adapted from SAVeTrEE and SST/NDTI Charts

//add comments in ui panel and readme for layer overlays

//--------------------------------------------------------------------
//                            Parameters
//--------------------------------------------------------------------
// Set and toggle parameters for analysis

// Gulf outline
var gulfOutline = ee.FeatureCollection('users/hmurayama/study_area');
// Study Area
var studyArea = ee.Geometry.Rectangle([-83.75, 8.35, -83.05, 8.95]);

// study years for bivariate map
var startYear = 2003;
var endYear = 2017;

// Bivariate color
var viz = ['eff5f5', 'b0d5df', '64acbe', 'e4acac', 
           'ad9ea5', '627f8c', 'c85a5a', '985356', '574249'];
           
////UI visualization
var BIV = {"opacity":1, "min":1, "max":9, "palette": viz};

// Set and enable panel
var panel = ui.Panel({style: {width:'200px'}});
ui.root.insert(0,panel);

// Stable variables
var bivariate,      // The bivariate layer to  be added to the map
    year,           // Years to loop through
    ndtiDiff_perc,  // NDTI percentiles for treshold 
    sstDiff_perc;   // SST percentiles for threshold

//--------------------------------------------------------------------
//                       Define functions                      
//--------------------------------------------------------------------
// Define bivariate function
// A function to specifically add the Bivariate layer to the map
function applyBivariate(year){
  // Classify NDTI & SST combinations
  var sstA = sstDiff.lt(ee.Number(sstDiff_perc.get("sst_p40")));
  var sstB = sstDiff.gte(ee.Number(sstDiff_perc.get("sst_p40"))).and(sstDiff.lt(ee.Number(sstDiff_perc.get("sst_p70"))));
  var sstC = sstDiff.gte(ee.Number(sstDiff_perc.get("sst_p70")));
  var ndti1 = ndtiDiff.lt(ee.Number(ndtiDiff_perc.get("nd_p40")));
  var ndti2 = ndtiDiff.gte(ee.Number(ndtiDiff_perc.get("nd_p40"))).and(ndtiDiff.lt(ee.Number(ndtiDiff_perc.get("nd_p70"))));
  var ndti3 = ndtiDiff.gte(ee.Number(ndtiDiff_perc.get("nd_p70")));
  // Create binary layers for each of the 9 classes
  var A1 = sstA.and(ndti1);var A2 = sstA.and(ndti2);var A3 = sstA.and(ndti3);
  var B1 = sstB.and(ndti1);var B2 = sstB.and(ndti2);var B3 = sstB.and(ndti3);
  var C1 = sstC.and(ndti1);var C2 = sstC.and(ndti2);var C3 = sstC.and(ndti3);
  // Create Class layers
  A1 = ee.Image.constant(1).updateMask(A1).unmask();
  A2 = ee.Image.constant(2).updateMask(A2).unmask();
  A3 = ee.Image.constant(3).updateMask(A3).unmask();
  B1 = ee.Image.constant(4).updateMask(B1).unmask();
  B2 = ee.Image.constant(5).updateMask(B2).unmask();
  B3 = ee.Image.constant(6).updateMask(B3).unmask();
  C1 = ee.Image.constant(7).updateMask(C1).unmask();
  C2 = ee.Image.constant(8).updateMask(C2).unmask();
  C3 = ee.Image.constant(9).updateMask(C3).unmask();
  // Combine into a single image
  var combo = A1.add(A2).add(A3).add(B1).add(B2).add(B3).add(C1).add(C2)
            .add(C3).clip(gulfOutline);
  bivariate = combo.select('constant').unmask(-9999);
  Map.centerObject(studyArea);
  Map.addLayer(bivariate.clip(gulfOutline), BIV,year.toString());
  print(bivariate.clip(gulfOutline));
}

// Apply Quality Assurance (QA) for the Landsat8 image collection
// (Adapted from Google developers tutorial)
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}


//A function to export Bivariate layer to drive
function exportBiv(year){ 
  applyBivariate(year);
  Export.image.toDrive({
    image: bivariate,
    description: year + '_Bivariate',
    fileNamePrefix: year +'_SST_NDTI_Bivariate', 
    folder: 'Bivariate SST NDTI',
    scale: 100,
    region: gulfOutline,
    skipEmptyTiles: true,
    crs: 'EPSG:3857'
  });
}

//--------------------------------------------------------------------
//                    Creat bivariate map by year
//--------------------------------------------------------------------
// Year loop
for (year=startYear; year<=endYear; year++) {
    print(year);
    
    if(year===2004 || year===2007) {
      print("Not enough Landsat data");
      continue;
    }
    
    var studyYear = year.toString();
    // Study start and ends will be updated according to set study year
    var studyRainStart = studyYear + '-02-01';
    var studyRainEnd = studyYear + '-05-31';
    var studyDryStart = studyYear + '-08-01';
    var studyDryEnd = studyYear + '-11-30';
    
  //--------------------------------------------------------------------
  //                                SST
  //--------------------------------------------------------------------
  // Load Aqua Modis Data and limit to study scope
  var aquaMODIS = ee.ImageCollection("NASA/OCEANDATA/MODIS-Aqua/L3SMI");
  var sstRain = aquaMODIS
    .select(['sst'])
    .filterDate(studyRainStart,studyRainEnd)
    .filterBounds(gulfOutline)
    .median();
  var sstDry = aquaMODIS
    .select(['sst'])
    .filterDate(studyDryStart,studyDryEnd)
    .filterBounds(gulfOutline)
    .median();
  // Calculate difference
  var sstDiff = sstDry.subtract(sstRain).abs().clip(gulfOutline);

  //Save yearly values for average later
  if (studyYear=="2015") {
    var sstDiff2015 = sstDiff;
  } else if (studyYear=="2016") {
    var sstDiff2016 = sstDiff;
  } else if (studyYear=="2017") {
    var sstDiff2017 = sstDiff;
  }

  //Calculate percentiles to obtain thresholds
  var sstDiff_perc = sstDiff.reduceRegion(ee.Reducer.percentile([10, 20, 30, 40, 50, 60, 70, 80, 90]), gulfOutline, 300);
  print(sstDiff_perc);

  // Get information about the MODIS projection.
  var modisProjection = ee.Image(aquaMODIS.filterDate(studyDryStart,studyDryEnd).first()).select('sst').projection();
  
  //--------------------------------------------------------------------
  //                                NDTI                               
  //--------------------------------------------------------------------
  // Load Landsat 8 image colection and filter date and region of images

  if (year < 2000) {
    var imc = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');
  } else if (year>=2000 & year<2014) {
    var imc = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
  } else {
    var imc = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
  }
  var landsat8 = imc;
  
  //RAIN
  var imcRain = landsat8
      .filterDate(studyRainStart, studyRainEnd)
      .filterBounds(gulfOutline)
      .map(maskL8sr);
      
  // Calculate median of all images    
  var imRain = imcRain.reduce(ee.Reducer.median());
  // Calculating Normalized Difference Turbidity Index (NDTI)
  var ndtiRain = imcRain.reduce(ee.Reducer.median()).normalizedDifference(['B4_median', 'B3_median']);
  
  //DRY
  var imcDry = landsat8
      .filterDate(studyDryStart, studyDryEnd)
      .filterBounds(gulfOutline)
      .map(maskL8sr);
    
  // Calculate median of all images    
  // Calculating Normalized Difference Turbidity Index (NDTI)
  var ndtiDry = imcDry.reduce(ee.Reducer.median()).normalizedDifference(['B4_median', 'B3_median']);
  
  //Difference
  var ndtiDiff = ee.Image(ndtiDry.subtract(ndtiRain).abs())
      .reproject(modisProjection)
      // Force the next reprojection to aggregate instead of resampling.
      .reduceResolution({
        reducer: ee.Reducer.median(),
        maxPixels: 1024
      })
      .clip(gulfOutline);
  
  //Save yearly values for average later
  if (studyYear=="2015") {
    var ndtiDiff2015 = ndtiDiff;
  } else if (studyYear=="2016") {
    var ndtiDiff2016 = ndtiDiff;
  } else if (studyYear=="2017") {
    var ndtiDiff2017 = ndtiDiff;
  }
  
  
  //Calculate percentiles to obtain thresholds
  ndtiDiff_perc = ndtiDiff.reduceRegion(ee.Reducer.percentile([10, 20, 30, 40, 50, 60, 70, 80, 90]), gulfOutline, 300);
  print(ndtiDiff_perc);

  
  //--------------------------------------------------------------------
  //                             Add to map                            
  //--------------------------------------------------------------------
  //Apply and Export NDTI and SST fluctuation maps
  exportBiv(studyYear);
  
  Export.image.toDrive({
  image: ndtiDiff.clip(studyArea).unmask(-9999),
  folder:'Bivariate SST NDTI',
  description: studyYear + '_NDTI',
  fileNamePrefix: studyYear +'_NDTI', 
  region: gulfOutline,
  scale: 100,
  });
  Export.image.toDrive({
    image: sstDiff.clip(studyArea).unmask(-9999),
    folder:'Bivariate SST NDTI',
    description: studyYear + '_SST',
    fileNamePrefix: studyYear +'_SST', 
    scale: 100,
    region: gulfOutline
  });
}


// Create average version
  print("Average");
  sstDiff = sstDiff2015.add(sstDiff2016).add(sstDiff2017).divide(3);
  var sstDiff_perc = sstDiff.reduceRegion(ee.Reducer.percentile([10, 20, 30, 40, 50, 60, 70, 80, 90]), gulfOutline, 300);
  print(sstDiff_perc);
  ndtiDiff = ndtiDiff2015.add(ndtiDiff2016).add(ndtiDiff2017).divide(3);
  var ndtiDiff_perc = ndtiDiff.reduceRegion(ee.Reducer.percentile([10, 20, 30, 40, 50, 60, 70, 80, 90]), gulfOutline, 300);
  print(ndtiDiff_perc);
  studyYear = "Average";
  exportBiv(studyYear);


// Create bivariate Legend
  // Text in legend
  var legendTitle = ui.Label('Bivariate Classification');
  var turbid = ui.Label('NDTI change',
    {margin: '0 0 0 10px',fontSize: '10px',color: 'gray', padding: '5'});
  var restoration = ui.Label('Poor restoration',
    {margin: '0 0 0 30px',fontSize: '10px',color: 'gray', padding: '5'});
  var textTop = ui.Panel([turbid, restoration], ui.Panel.Layout.flow('horizontal'));
  var temp = ui.Label('SST change',
    {margin: '5px 0 0 30px',fontSize: '10px',color: 'gray'});
  var good = ui.Label('Restoration',
    {margin: '5px 0 0 20px',fontSize: '10px',color: 'gray'});
  var textBottom = ui.Panel([good, temp], 
      ui.Panel.Layout.flow('horizontal'));
  panel.add(textTop);
  
  //Add 9 color boxes
  var makeRow = function(color1, color2, color3) {
    // Create the label that is actually the colored box.
    var colorBox1 = ui.Label({
      style: {
        backgroundColor: '#' + color1,
        padding: '10px',
        margin: '0 0 0 40px'
      }
    });
    var colorBox2 = ui.Label({
      style: {
        backgroundColor: '#' + color2,
        padding: '10px',
        margin: '0 0 0 0'
      }
    });
    var colorBox3 = ui.Label({
      style: {
        backgroundColor: '#' + color3,
        padding: '10px',
        margin: '0 0 0 0'
      }
    });
  
    return ui.Panel({
      widgets: [colorBox1, colorBox2, colorBox3],
      layout: ui.Panel.Layout.Flow('horizontal')
    });
  };
  
    for (var i = 2; i > -1; i = i-1) {
      panel.add(makeRow(viz[i], viz[i+3], viz[i+6]));
    }
  
  panel.add(textBottom);


  //Description text
  var description = ui.Label('The bivariate maps are overlapped with 100% opacity. Toggle through each layer using the Layers button.');
  panel.add(description);




//--------------------------------------------------------------------
//                          Add Land Mask
//--------------------------------------------------------------------
  // Add land mask (obtained and slightly modified from GEE tutorial)
  // Load or import the Hansen et al. forest change dataset.
  var hansenImage = ee.Image('UMD/hansen/global_forest_change_2015');
  // Select the land/water mask.
  var datamask = hansenImage.select('datamask');
  // Create a binary mask (1 is land and 2 is water)
  var maskl = datamask.eq(1);
  // Update the composite mask with the water mask.
  var maskedLand =landsat8.mean().updateMask(maskl).select('B3');
  var maskedNDTI =ndtiDiff.updateMask(maskl).unmask(-9999);
  // Visualize the defined mask
  var vismask = {bands: ['B3'], max: 0.5,"palette":["000000"], };
  
  Map.addLayer(maskedLand.clip(studyArea), vismask, 'masked');
  Export.image.toDrive({
    folder:'Bivariate SST NDTI',
    image: maskedLand,
    description: 'mask_data',
    region: studyArea,
    scale: 30,
  });
  Export.image.toDrive({
    folder:'Bivariate SST NDTI',
    image: maskedNDTI,
    description: 'mask_vis',
    region: studyArea,
    scale: 30,
  });



//--------------------------------------------------------------------
//                  Load Datasets / Define Parameters
//--------------------------------------------------------------------

// Load Landsat8 dataset
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
     .filterBounds(gulfOutline)
     .filterDate('2014-01-01', '2017-12-31')
     .map(maskL8sr);

var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
     .filterBounds(gulfOutline)
     .filterDate('2003-01-01', '2013-12-31')
     .map(maskL8sr);

// Load AQUA-MODIS dataset (2014-2017)
var sst1 = ee.ImageCollection("NASA/OCEANDATA/MODIS-Aqua/L3SMI")
     .select(['sst'])
     .filterBounds(gulfOutline)
     .filterDate('2014-01-01', '2017-12-31');

// Load AQUA-MODIS dataset (2003-2013)
var sst2 = ee.ImageCollection("NASA/OCEANDATA/MODIS-Aqua/L3SMI")
     .select(['sst'])
     .filterBounds(gulfOutline)
     .filterDate('2003-01-01', '2013-12-31');
// Map a function over the Landsat 8 TOA collection to add an NDTI band.
var l8NDTI = l8.map(function(image) {
  var ndti8 = image.normalizedDifference(['B4', 'B3']).rename('NDTI');
  return image.addBands(ndti8);
});

var l7NDTI = l7.map(function(image) {
  var ndti7 = image.normalizedDifference(['B3', 'B2']).rename('NDTI');
  return image.addBands(ndti7);
});

//--------------------------------------------------------------------
//                         Create charts
//--------------------------------------------------------------------

// SST Chart (2014 - 2017)
var chartSST1 = ui.Chart.image.series({
  imageCollection: sst1,
  region: gulfOutline,
  reducer: ee.Reducer.median(),
}).setOptions({
  title: 'SST Time Series (2014 - 2017)',
  vAxis: {title:'SST'},
  hAxis: {title: 'Date (MM-yy)', format: 'MM-yy', gridlines: {count: 12}}
});

// SST Chart (2003 - 2013)
var chartSST2 = ui.Chart.image.series({
  imageCollection: sst2,
  region: gulfOutline,
  reducer: ee.Reducer.median(),
}).setOptions({
  title: 'SST Time Series (2003 - 2013)',
  vAxis: {title:'SST'},
  hAxis: {title: 'Date (MM-yy)', format: 'MM-yy', gridlines: {count: 12}}
});

// NDTI Chart (2014 - 2017)
var chartNDTI8 = ui.Chart.image.series({
  imageCollection: l8NDTI.select('NDTI'),
  region: gulfOutline,
  reducer: ee.Reducer.median(),
}).setOptions({
  title: 'NDTI Time Series (2014 - 2017)',
  vAxis: {title:'NDTI'},
  hAxis: {title: 'Date (MM-yy)', format: 'MM-yy', gridlines: {count: 12}}
});

// NDTI Chart (2003 - 2013)
var chartNDTI7 = ui.Chart.image.series({
  imageCollection: l7NDTI.select('NDTI'),
  region: gulfOutline,
  reducer: ee.Reducer.median(),
}).setOptions({
  title: 'NDTI Time Series (2003 - 2013)',
  vAxis: {title:'NDTI'},
  hAxis: {title: 'Date (MM-yy)', format: 'MM-yy', gridlines: {count: 12}}
});


// Display the chart in the console.
print('please wait...'),
print('plotting charts might take a while'),
print(chartSST1);
print(chartSST2);
print(chartNDTI8);
print(chartNDTI7);