WET2.0_Training.txt

//////// WET TRAINING SCRIPT //////// 
// Wetland Extent Tool (W.E.T) Training Script
// This script allows a user to create a custom image stack to use for training in the WET 2.0 tool.
// It utilizes Landsat 8 surface reflectance, Sentinel-2 surface reflectance, Sentinel-1 C-SAR, and elevation data.
// The script takes in a date range and area of interest as inputs and uploads an image stack as a GEE asset
// NDVI, MNDWI, and TCWGD are calculated for both Landsat 8 and Sentinel 2, DSWE is calculated for Landsat 8 only
// VV in natural values, VH in natural values, and the VV/VH ratio in natural values are calculated for Sentinel 1
// The user can edit what bands to include in the output image stack
// The script calculates all of these indices, applies a SNIC algorithm, and creates a stack of the clustered parameters.
// Once uploaded as a GEE asset, the user can import the new stack in WET 2.0 for a more customized classification

var training_data = ee.FeatureCollection('users/ericacarcelen/allBasin_fieldPoly_v2'),
    srtm = ee.Image("USGS/SRTMGL1_003"),
    studyarea = ee.FeatureCollection('users/ericacarcelen/2020spring_jpl_glwii_studyarea'),
    basins = ee.FeatureCollection('users/ericacarcelen/greatlakes_subbasins');



//////////////////////////////////////////////////////////////////////////////////
// Constants /////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////

// Define dates ranges to filter datasets.
var y2019 = ee.Filter.date('2019-01-01', '2019-12-31');


          var date = y2019; // CHANGE DATE HERE


// Define area of interest.

  // GL2 team created training stack for entire GL Basin
  // GL Basin was too large to process, following variables select individual lake basin polygons in the basins shapefile
  var erie = basins.filter(ee.Filter.eq('merge','lk_erie'))       //ee.Filter.eq tells GEE to select a feature where
  var superior = basins.filter(ee.Filter.eq('merge','lk_sup'))    //the field (first input) equals a specified attribute (second input)
  var huron = basins.filter(ee.Filter.eq('merge','lk_huron'))     //Ex:filter basins shapefile where feature equals lk_huron in the merge field
  var michigan = basins.filter(ee.Filter.eq('merge','lk_mich'))
  var ontario = basins.filter(ee.Filter.eq('merge','lk_ont'))


          var aoi = ontario // CHANGE AREA OF INTEREST HERE, CONSIDER DIVIDING AREA IF TOO LARGE OF AN AREA
          Map.centerObject(aoi)
          

// Define bands and indices to include in final stack. Dictionary needs to be defined because
// some functions result in the renaming of bands/indices.
// These are all the inputs/indices calculated by the GL1 and GL2 teams. If you are adding inputs other than
// what is already listed, follow the naming convention below and keep the order the same.
var post_snic_names = ['clusters', 'NDVI_mean', 'MNDWI_mean', 'ratio_mean', 'VV_mean','VH_mean',
    'TCWGD_mean','B4_mean','B3_mean','B2_mean','s2NDVI_mean','s2MNDWI_mean','s2TCWGD_mean',
    's2B4_mean', 's2B3_mean', 's2B2_mean', 'DSWE_mean']
var pre_snic_names =   ['clusters', 'NDVI', 'MNDWI', 'ratio', 'VV','VH','TCWGD','B4','B3','B2', 's2NDVI','s2MNDWI','s2TCWGD',
    's2_B4', 's2_B3', 's2_B2', 'DSWE']
var dictionary = ee.Dictionary.fromLists(pre_snic_names, post_snic_names)

// CHANGE BANDS HERE. Possible bands/indices include: 
//    Landsat 8: TCWGD, NDVI, MNDWI, B4 (red), B3 (green), B2 (blue), DSWE
//    Sentinel 2: TCWGD, NDVI, MNDWI, B4 (red), B3 (green), B2 (blue)
//    Sentinel 1: VV, VH, VV/VH ratio

            var l8_bands = ee.List(['DSWE'])
            var s2_bands = ee.List(['s2TCWGD','s2MNDWI','s2NDVI'])
            var s1_bands = ee.List(['VV','VH', 'ratio'])

var bands = l8_bands.cat(s2_bands).cat(s1_bands) // This variable is a list of all bands
var snic_bands = dictionary.values(bands) // Allows for bands to always be defined, even with name change

print('Bands and Indices Used in Classification:', bands)



///////////////////////////////////////////////////////////////////////////////////
// Sentinel 1 /////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

// Functions ----------------------------------------------------------------------

// Create function to convert to natural values.
function toNatural(image) {
  return ee.Image(10.0).pow(image.divide(10.0));
}

// Create functions to remove dark borders, threshold value depends on VV or VH.
// Thresholds were chosen based on trial an error and are different depending on bands.
// Spot check if area of interest changes significantly.

function removeBordersVV(image) {
  var vv = image.select('VV');
  var vv_edge = vv.lt(0.005); // Can change this threshold if needed
  var vv_maskedImage = vv.mask().and(vv_edge.not());
  return image.addBands(vv.updateMask(vv_maskedImage),['VV'],true);
}
function removeBordersVH(image) {
  var vh = image.select('VH');
  var vh_edge = vh.lt(0.0008); // Can change this threshold if needed
  var vh_maskedImage = vh.mask().and(vh_edge.not());
  return image.addBands(vh.updateMask(vh_maskedImage),['VH'],true);
}

// Function to add VV/VH band to each image in the collection.
function bandRatio(image){
  var vv = image.select('VV');
  var vh = image.select('VH');
  var ratio = vv.divide(vh).rename('ratio');
  return image.addBands(ratio);
}

// Data collection ----------------------------------------------------------------

// Filter, convert, and mask Sentinel-1 image collection.
var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  .filterBounds(aoi)
  .filter(date)
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
  .select('VV', 'VH')


// Create and print list of images in collection to the console.
print('Sentinel-1 C-SAR Images Used:', s1);

// Data processing ----------------------------------------------------------------

// Convert from dB to natural numbers, add VV/VH ratio as a band, remove borders
// (note that this may remove open water pixels).

var s1_processed = s1
  .map(toNatural)
  .map(removeBordersVV)
  .map(removeBordersVH)
  .map(bandRatio)

// Reduce VV/VH collection to mean values.
var s1_composite = s1_processed
  .reduce(ee.Reducer.mean())
  .select(['ratio_mean','VV_mean','VH_mean'],['ratio','VV','VH']);


// Map.addLayer(s1_composite.clip(aoi), {
//   bands: 'ratio',
//   min: 1, 
//   max: 31, 
//   gamma: 2}, 
//   'S1', 
//   false);

///////////////////////////////////////////////////////////////////////////////////
// Landsat 8 //////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

// Functions ----------------------------------------------------------------------

// Mask clouds using bits 3 and 5 of the pixel_qa band 
var cloudMask = function(image){
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;
  var qa = image.select('pixel_qa');
  var qa_mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(qa_mask);
};

// Function to add NDVI band to each image in the collection.
var addNDVI = function(image){
  var NDVI = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(NDVI);
};

// Function to add MNDWI band to each image in the collection.
var addMNDWI = function(image){
  var MNDWI = image.normalizedDifference(['B3', 'B6']).rename('MNDWI');
  return image.addBands(MNDWI);
};

// Function to add Tasseled Cap Wetness Greenness Difference to each image in collection
var addTCWGD = function(image){
  var TCW = image.expression(
    '(0.1511*BLUE)+(0.1973*GREEN)+(0.3283*RED)+(0.3407*NIR)+(-0.7117*SWIR1)+(-0.4559*SWIR2)',{
      'BLUE': image.select('B2'),
      'GREEN': image.select('B3'),
      'RED': image.select('B4'),
      'NIR': image.select('B5'),
      'SWIR1': image.select('B6'),
      'SWIR2': image.select('B7')
    })

  var TCG = image.expression(
      '(-0.2941*BLUE)+(-0.243*GREEN)+(-0.5424*RED)+(0.7276*NIR)+(0.0713*SWIR1)+(-0.1608*SWIR2)',{
        'BLUE': image.select('B2'),
        'GREEN': image.select('B3'),
        'RED': image.select('B4'),
        'NIR': image.select('B5'),
        'SWIR1': image.select('B6'),
        'SWIR2': image.select('B7')
      })
  var TCWGD = TCW.subtract(TCG).rename('TCWGD');
  return image.addBands(TCWGD);
};

// Function to calculate and add DSWE
var addDSWE = function(image, dem){
  // define function to add index as band to each image in collection
  var addNDVI = function(image) {
    var ndvi = image.normalizedDifference(['B5','B4']).rename('NDVI')
    return image.addBands(ndvi);
  };
  var addMNDWI = function(image){
    var mndwi = image. normalizedDifference(['B3','B6']).rename('MNDWI')
    return image.addBands(mndwi);
  };
  var addMBSRV = function(image){
    var mbsrv = image.select('B3').add(image.select('B4')).rename('MBSRV')
    return image.addBands(mbsrv);
  };
  var addMBSRN = function(image){
    var mbsrn= image.select('B5').add(image.select('B6')).rename('MBSRN')
    return image.addBands(mbsrn);
  };
  var addAWESH = function(image) {
  var awesh = (image
    .expression('Blue + (2.5 * Green) - (1.5 * mbsrn) - (0.25 * Swir2)', {
      'Blue': image.select(['B2']),
      'Green': image.select(['B3']),
      'mbsrn': addMBSRN(image).select(['MBSRN']),
      'Swir2': image.select(['B7']) })
      .rename('AWESH')
      );
    return image.addBands(awesh);
  };
  
  //apply index functions to image collection
  var dswe_indices = image.map(addNDVI)
                          .map(addMNDWI)
                          .map(addMBSRV)
                          .map(addMBSRN)
                          .map(addAWESH);
  
  //calculate mean composite for each band
  //output stack containing only the indices/bands needed for dswe
  var dswe_inputs = dswe_indices.mean()
                                .select(['MNDWI','MBSRV','MBSRN',
                                        'AWESH','NDVI','B2',
                                        'B5','B6','B7']);
                                         
  // define variables 
  var mndwi = dswe_inputs.select('MNDWI'),
      mbsrv = dswe_inputs.select('MBSRV'),
      mbsrn = dswe_inputs.select('MBSRN'),
      awesh = dswe_inputs.select('AWESH'),
      ndvi = dswe_inputs.select('NDVI'),
      swir1 = dswe_inputs.select('B6'),
      nir = dswe_inputs.select('B5'),
      blue = dswe_inputs.select('B2'),
      swir2 = dswe_inputs.select('B7'),
      slope = ee.Terrain.slope(dem),
      hillshade = ee.Terrain.hillshade(dem); 
  
  // define thresholds for each test
  var t1_thresh = mndwi.gt(0.124).rename('Test_1');
  var t2_thresh = mbsrv.gt(mbsrn).rename('Test_2');
  var t3_thresh = awesh.gt(0).rename('Test_3');
  var t4_thresh = mndwi.gt(-0.44)
                      .and(swir1.lt(900))
                      .and(nir.lt(1500))
                      .and(ndvi.lt(0.7))
                      .rename('Test_4');
  var t5_thresh = mndwi.gt(-0.5)
                      .and(blue.lt(1000))
                      .and(swir1.lt(3000))
                      .and(swir2.lt(1000))
                      .and(nir.lt(2500))
                      .rename('Test_5');
  
  //multiply booleans from test thresholds and add to create pixel values
  var tests = t1_thresh                        //test 1 true = 1
              .add(t2_thresh.multiply(10))     //test 2 true = 10
              .add(t3_thresh.multiply(100))    //test 3 true = 100
              .add(t4_thresh.multiply(1000))   //test 4 true = 1000
              .add(t5_thresh.multiply(10000))  //test 5 true = 10000
              .rename('DSWE_Tests');           //rename result 
  
  //apply decision rules to classify pixel in 5 classes
  var rules = tests.remap(
    /**first list - pixel values in test image, 
    * second list: new values for dswe output in following order:
      * 1: no water
      * 2: water-high
      * 3: water-moderate
      * 4: potential wetland
      * 5: water-low
    */
    [0, 1, 10, 100, 1000, 
    1111, 10111, 11011, 11101, 11110, 11111, 
    111, 1011, 1101, 1110, 10011, 10101, 10110, 11001, 11010, 11100, 
    11000,
    11, 101, 110, 1001, 1010, 1100, 10000, 10001, 10010, 10100], 
    [1, 1, 1, 1, 1, 
    2, 2, 2, 2, 2, 2, 
    3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
    4, 
    5, 5, 5, 5, 5, 5, 5, 5, 5, 5], 
    0)
    .rename('DSWE_Rules');
     
  //post-processing of image to apply hillshade/slope thresholds
  /**reclassifies to 1 (not water) if (in same order as coded):
      * water-high and slope >=30
      * water-moderate and slope >= 30
      * potential wetland and slope >= 20
      * water-low and slope >= 10
      * anywhere hillshade <= 110
  */
  var dswe =
  rules.where(slope.gte(30).and(rules.eq(2)), 1)
  rules.where(slope.gte(30).and(rules.eq(3)), 1)
  rules.where(slope.gte(20).and(rules.eq(4)), 1)
  rules.where(slope.gte(10).and(rules.eq(5)), 1)
  rules.where(hillshade.lte(110), 1);
  
  return(dswe.rename('DSWE'));
};

// Data collection ---------------------------------------------------------------

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filter(date)
  .filterBounds(aoi);

print('Landsat 8 OLI Images Used:', l8);

// Data processing ---------------------------------------------------------------

// Apply Cloud Mask
var l8_masked = l8.map(cloudMask);

// Map the NDVI, MNDWI and TCWGD indices.
var l8_indices = l8_masked
  .map(addNDVI)
  .map(addMNDWI)
  .map(addTCWGD);
  
// Calculate DSWE
var l8_DSWE = addDSWE(l8_masked, srtm).float();

// Make a composite of the mean. Select only the bands you need at this point to
// reduce computing time when this goes into SNIC.
var l8_composite = l8_indices.mean().addBands(l8_DSWE).select(l8_bands);


// Map.addLayer(l8_composite.clip(aoi),[],'L8', false)

///////////////////////////////////////////////////////////////////////////////////
// Sentinel-2 /////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

// Functions ----------------------------------------------------------------------

// Mask clouds using bits 10 and 11 of the QA60 band 
var s2cloudMask = function(image) {
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  var qa = image.select('QA60');
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
};

// Function to add NDVI band to each image in the collection.
var s2addNDVI = function(image){
  var NDVI = image.normalizedDifference(['B8', 'B4']).rename('s2NDVI');
  return image.addBands(NDVI);
};

// Function to add MNDWI band to each image in the collection.
var s2addMNDWI = function(image){
  var MNDWI = image.normalizedDifference(['B3', 'B11']).rename('s2MNDWI');
  return image.addBands(MNDWI);
};

// Function to add Tasseled Cap Wetness Greenness Difference to each image in collection
var s2addTCWGD = function(image){
  var TCW = image.expression(
    '(0.2578*BLUE)+(0.2305*GREEN)+(0.0883*RED)+(0.1071*NIR)+(-0.7611*SWIR1)+(-0.5308*SWIR2)',{
      'BLUE': image.select('B2'),
      'GREEN': image.select('B3'),
      'RED': image.select('B4'),
      'NIR': image.select('B8'),
      'SWIR1': image.select('B11'),
      'SWIR2': image.select('B12')
    })

  var TCG = image.expression(
      '(-0.3599*BLUE)+(-0.3533*GREEN)+(-0.4734*RED)+(0.6633*NIR)+(0.0087*SWIR1)+(-0.2856*SWIR2)',{
        'BLUE': image.select('B2'),
        'GREEN': image.select('B3'),
        'RED': image.select('B4'),
        'NIR': image.select('B8'),
        'SWIR1': image.select('B11'),
        'SWIR2': image.select('B12')
      })
  var TCWGD = TCW.subtract(TCG).rename('s2TCWGD');
  return image.addBands(TCWGD);
};

      
// Data collection ----------------------------------------------------------------
var s2 = ee.ImageCollection("COPERNICUS/S2_SR")
  .filter(date)
  .filterBounds(aoi)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 10);

print('Sentinel 2  MSI Level-2 Images Used:', s2);

// Data processing ----------------------------------------------------------------

// Apply Cloud Mask
var s2_masked = s2.map(s2cloudMask);

// Map the NDVI, MNDWI and TCWGD indices.
var s2_indices = s2_masked
  .map(s2addNDVI)
  .map(s2addMNDWI)
  .map(s2addTCWGD);

// Make a composite of the mean. Select only the bands you need at this point to
// reduce computing time when this goes into SNIC.
var s2_composite = s2_indices.mean().select(s2_bands);

// Map.addLayer(s2_composite.clip(aoi),[], 'S2', false)

//////////////////////////////////////////////////////////////////////////////////
// Assembling image stack ////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////

// This section below, is a image stack made without changing resolutions. It works fine
// when displaying data accross the whole state.
var stack = l8_composite
  .addBands(s2_composite)
  .addBands(s1_composite)
  .clip(aoi);

//////////////////////////////////////////////////////////////////////////////////
// Object Based Image Analysis ///////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////

//Image segmentation using Simple Non-Iterative Clustering (SNIC).
var seeds = ee.Algorithms.Image.Segmentation.seedGrid(10);

// Run SNIC on image stack. The parameters below are adjustable.
var SNIC = ee.Algorithms.Image.Segmentation.SNIC({
  image: stack.clip(aoi), 
  compactness: 0.01, // values closer to one produces "squarer" clusters
  connectivity: 8,
  neighborhoodSize: 250,
  seeds: seeds
})

//Reproject SNIC layer so that it doesn't change with different zoom levels
var l8_projection = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_025027_20170524').projection()
var SNIC_fixed = SNIC.select(snic_bands,bands).reproject(l8_projection,null,30)


// Outlines of clusters for evaluation
var minMax = SNIC_fixed.reduceNeighborhood(ee.Reducer.minMax(), ee.Kernel.square(1));
var perimeterPixels = minMax.select(0).neq(minMax.select(1)).rename('perimeter');
var perimeter = perimeterPixels.addBands(SNIC)
    .reduceConnectedComponents(ee.Reducer.sum(), 'clusters', 256);
var perimeter_outlines = perimeterPixels.remap([1],[1]);
// Map.addLayer(perimeter_outlines, [], 'Cluster outlines', false);

// Generate image stack excluding the clusters layer.
var slope = ee.Terrain.slope(srtm)
var stack = stack.addBands(slope).addBands(srtm)
var stack = stack.float();
print(stack)

// Map layers if you want to visualize
Map.addLayer(stack, {}, 'Training Multiband Image')


////////////////////////////////////////////////////////////////////////////////////
//// Image Export //////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////

// Export Training Image Stack 

// Export as Asset, this will store the image stack in your personal GEE repository
// In order to use the image stack in GEE, it must be stored as an Asset
Export.image.toAsset({                              // LEAVE IMAGE AS STACK, CHANGE REMAINING PARAMETERS AS NEEDED
  image: stack,                                     // DO NOT CHANGE
  description: 'Ontario_trainingStack_s2L2v2_ssn2019',   // CHANGE DESCRIPTION TO MATCH YOUR AOI AND DATE
  scale: 30,                                        // DEFINES PIXEL RESOLUTION
  region: aoi,                                      // DEFINES OUTPUT TIFF EXTENT
  assetId: 'ont_trainStack_L2v2_ssn2019',                    // "FILENAME" IN GEE, WILL NOT UPLOAD IF ASSET WITH SAME NAME EXISTS
  maxPixels:  1e13                            // DEFINES MAXIMUM NUMBER OF PIXELS ALLOWED, MAY NEED TO CHANGE FOR YOUR AOI**
});                                                 // **IF ERROR REGARDING BYTE LIMIT RETURNS, AOI TOO LARGE TO PROCESS, MUST DIVIDE INTO SMALLER AOIS

// Export to Google Drive, this will store the image stack in your personal Google Drive
// you can download the geotiff from your Google Drive to store on your computer
// however, the image stack must be stored as an asset for use in GEE (and the WET 2.0 tool)
Export.image.toDrive({
  image: stack, 
  description: 'Ontario_trainingStack_s2L2_3262020',
  scale: 30,
  region: aoi,
  maxPixels:1e13
});