Originally published in RS

// JavaScript code to be implemented in Google Earth Engine(c) developed by H.A. Orengo
// to accompany the paper:

// Orengo, H.A. and Petrie, C.A. 2017. Large-scale, multi-temporal remote sensing of
// palaeo-river networks: a case study from northwest India and its implications for
// the Indus Civilisation. Remote Sensing, 9(7): 735; doi: 10.3390/rs9070735

//                ------------- ooo -------------

// TO EXECUTE THESE SCRIPTS PASTE THIS CODE INTO GOOGLE EARTH ENGINE CODE EDITOR AND PRESS RUN

//                ------------- ooo -------------

// The different algorithms are written here following the order adopted in the journal article:
// 1. NDVSI (Normalised Difference Vegetation Seasonality Index), line 78
// 2. SMTVI (Seasonal Multi-Temporal Vegetation Indices),         line 128
// 3. Wet Months Multi-Temporal Tasselled Cap Transformation,     line 176
// 4. Dry Months Multi-Temporal Tasselled Cap Transformation,     line 235
// 5. Wet Months Multi-Temporal Principal Component Analysis,     line 286
// 6. Dry Months Multi-Temporal Principal Component Analysis,     line 391

// For more information on how these scripts work and how to apply them refer to the text of the 
// article, which is freely available at: http://www.mdpi.com/2072-4292/9/7/735
// Suggestions and code improvementents are welcome! Please, contact Hector A. Orengo

// NOTES, READ CAREFULLY!
// 1. Visualisation parameters for each layer can be adjusted in the Layers dialogue in the map
//    area of Goggle Earth Engine. This might be convenient when visualising at different scales.
// 2. To apply these analyses in other areas simply change the 'geometry' variable and the
//    'Map.setCenter' coordinates. These two parameters can be deleted. Doing so will apply the
//    algorithms in whatever area is being displayed in the Google Earth Engine's Map area.
//    However, if the geometry variable is eliminated it will need to be also deleted from the
//    scripts incorporating it (for example when it is employed to clip the analysis area).
//    The user can also draw his/her own polygon using the tool provided in the top left of Google
//    Earth Engine’s map window. This will create a polygon named ‘geometry’ by default. The user
//    can then delete the variable called ‘geometry’ in the code lines 52-65 and the analysis will
//    be performed in the newly defined area. Make sure the central coordinates defined by
//    'Map.setCenter' (lines 67-69) are also changed so the map window will zoom to these
//    coordinates when executing the scripts. 
// 3. The execution of the algorithms can take some time (particularly the PCAs). To shorten their
//    visualisation, select only those you require in the layers panel of the Map window.
// 4. The images resulting from these scripts can be transferred to your Google Drive running the
//    analysis from the 'task' tab on the top right of the screen.
// 5. All these algorithms can be run separately, however, there are variables employed by several
//    of those that only appear the first time that they are required and are not repeated in
//    subsequent algorithms. The reader is advised to check carefully which variables are employed
//    by each algorithm and to what part of the code they refer before attempting to run a single
//    algorithm separately from the rest of the code written here.


//                ------------- ooo -------------


// In the section below a set of variables defining (1) the collection of all Landsat 5 images 
// available, (2) the Area of Interest (AOI) and (3) the central coordinate for visualisation
// will be defined for the use of the different algorithms

// Set all the available Landsat 5 calibrated top-of-atmosphere reflectance images as a variable
var L5_TOA = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
        .filter(ee.Filter.lte('CLOUD_COVER', 100)); //ideally less than a 2% cloud cover should be selected however, this woudl reduce very much the number of images available. THe previous version of L5 TOA incorporated Fmask, which covered clouds, check how to do this again.

// Define your AOI as a polygon with a set of WGS84 coordinates
var geometry = ee.Geometry.Polygon(
    [[[77.59,28.8],
      [77.59,30.29],
      [77.06,30.56],
      [76.83,30.74],
      [76.55,31.09],
      [74.45,31.09],
      [73.84,30.21],
      [73.84,28.8]]]);

// Define a central point in your study area (as X,Y WGS84 decimal degrees) and a scale,
// just for visualisation purposes
Map.setCenter(75.564, 29.874, 8);


//                ------------- ooo -------------


// NDVSI (Normalised Difference Vegetation Seasonality Index)

// Load a Landsat 5 collection (1984-2013) and select the bands of interest: B3 (red) and B4 (NIR)
var L5_NDVSI = L5_TOA.select(['B3', 'B4']);

// Select the L5 images corresponding to the wet months and average their values
var wet1 = L5_NDVSI.filter(ee.Filter.dayOfYear(1,120))
    .mean();
var wet2 = L5_NDVSI.filter(ee.Filter.dayOfYear(181,240))
    .mean();
var wet = (wet1.add(wet2)).divide(2);

// Develop a simple vegetation ratio for the wet months
var dvi_wet = wet.expression(
    'b("B4") / (b("B4") + b("B3"))');

// Select the L5 images corresponding to the dry months and average their values
var dry1 = L5_NDVSI.filter(ee.Filter.dayOfYear(121,180))
    .mean();
var dry2 = L5_NDVSI.filter(ee.Filter.dayOfYear(241,360))
    .mean();
var dry = (dry1.add(dry2)).divide(2);

// Develop a simple vegetation ratio for the dry months
var dvi_dry = dry.expression(
    'b("B4") / (b("B4") + b("B3"))');

// Develop the NDVSI and clip it using the AOI
var ndvsi = ((dvi_dry.subtract(dvi_wet)))
    .divide((dvi_dry.add(dvi_wet)))
    .clip(geometry);

var ndvsi_palette =
    '011301, 011d01, 012e01, 023b01, 004c00, 056201, 207401, 3e8601, 529400,' +
    '66a000, 74a901, 99b718, fcd163, f1b555, df923d, ce7e45, 9d350e';

Map.addLayer(ndvsi, {min: -0.0705, max: 0.00408, palette: ndvsi_palette}, 'NDVSI');

Export.image.toDrive({
  image: ndvsi,
  description: 'NDVSI',
  scale: 30,
  maxPixels: 1e9,
  region: geometry
});


//                ------------- ooo -------------


// SMTVI (Seasonal Multi-Temporal Vegetation Indices)

// Import the image collection with all Landsat 5 8-days EVI composites  
var eviL5 = ee.ImageCollection('LANDSAT/LT5_L1T_8DAY_EVI');

// Filter them by two-moths periods and extract the average values 
var eviL5_JanFeb = eviL5.filter(ee.Filter.dayOfYear(1,60))
    .mean();
var eviL5_MarApr = eviL5.filter(ee.Filter.dayOfYear(61,120))
    .mean();
var eviL5_MayJun = eviL5.filter(ee.Filter.dayOfYear(121,180))
    .mean();
var eviL5_JulAug = eviL5.filter(ee.Filter.dayOfYear(181,240))
    .mean();
var eviL5_SepOct = eviL5.filter(ee.Filter.dayOfYear(241,300))
    .mean();
var eviL5_NovDec = eviL5.filter(ee.Filter.dayOfYear(301,360))
    .mean();

// Create a composite image with all the two-month average EVI composites and clip the area of analysis
// according to your AOI (the AOI is defined by the geometry variable, which can be created/modified using
// the map window below)
var composite = ee.Image([eviL5_JanFeb, eviL5_MarApr, eviL5_MayJun, eviL5_JulAug, eviL5_SepOct, eviL5_NovDec])
    .clip(geometry);

// Rename the composite bands so they can be easier to interpret
var SMTVI = composite.select(
    ['EVI', 'EVI_1', 'EVI_2', 'EVI_3', 'EVI_4', 'EVI_5'], // old names
    ['EVI_JanFeb', 'EVI_MarApr', 'EVI_MayJun', 'EVI_JulAug', 'EVI_SepOct', 'EVI_NovDec'] // new names
);

// Add the newly created layer to the map window. In this case we have created a RGB composite joining
// the bi-month EVI average values from July-August (R), March-April (G) and January-February (B) 
Map.addLayer(SMTVI, {bands: ["EVI_JulAug","EVI_MarApr","EVI_JanFeb"], min: 0.12340243216502006, max: 0.5330676407686823}, 'SMTVI');

// Export the SMTVI image (with all bi-month EVI averages) to your Google drive
Export.image.toDrive({
  image: SMTVI,
  description: 'SMTVI',
  scale: 30,
  maxPixels: 1e9,
  region: geometry
});


//                ------------- ooo -------------


// Wet Months Multi-Temporal Tasselled Cap Transformation

// Define an Array of Tasselled Cap coefficients.
var coefficients = ee.Array([
  [0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863],
  [-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800],
  [0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572],
  [-0.8242, 0.0849, 0.4392, -0.0580, 0.2012, -0.2768],
  [-0.3280, 0.0549, 0.1075, 0.1855, -0.4357, 0.8085],
  [0.1084, -0.9022, 0.4120, 0.0573, -0.0251, 0.0238]
]);

// Load a Landsat 5 collection (1984-2013) and select the bands of interest
var L5_TOAbands = L5_TOA.select(['B1', 'B2', 'B3', 'B4', 'B5', 'B7']);

// Select the L5 images corresponding to the wet months, average their values and clip to AOI 
var L5_TOAbands_Wet1 = L5_TOAbands.filter(ee.Filter.dayOfYear(1,120))
    .mean();
var L5_TOAbands_Wet2 = L5_TOAbands.filter(ee.Filter.dayOfYear(181,240))
    .mean();
var L5_TOAbands_Wet = (L5_TOAbands_Wet1.add(L5_TOAbands_Wet2)).divide(2)
    .clip(geometry);

// Make an Array Image, with a 1-D Array per pixel.
var arrayImage1D_Wet = L5_TOAbands_Wet.toArray();

// Make an Array Image with a 2-D Array per pixel, 6x1.
var arrayImage2D_Wet = arrayImage1D_Wet.toArray(1);

// Do a matrix multiplication: 6x6 times 6x1.
var TCT_LS5_Wet = ee.Image(coefficients)
  .matrixMultiply(arrayImage2D_Wet)
  // Get rid of the extra dimensions.
  .arrayProject([0])
  .arrayFlatten(
    [['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth']]);

// Display the first three bands of the result and the input imagery.
var vizTCTwet = {
  bands: ['brightness', 'greenness', 'wetness'],
  min: [0.4299,0.003,-0.001655], max: [0.5304,0.0672975,0.0623]
};

// Add the resulting layer
Map.addLayer(TCT_LS5_Wet, vizTCTwet, 'TCT_LS5_Wet');

// Export the Multitemporal TCT of the wet months to your Google Drive
Export.image.toDrive({
  image: TCT_LS5_Wet,
  description: 'TCT_L5_Wet',
  scale: 30,
  maxPixels: 1e9,
  region: geometry
});


//                ------------- ooo -------------


// Dry Months Multi-Temporal Tasselled Cap Transformation

// NOTES!
// 1. The Array of Tasselled Cap coefficients has already been defined in the previous section
// 2. The bands of interest from the Landsat 5 collection have already been selected in the
//    previous section

// Select the L5 images corresponding to the dry months, average their values and clip to AOI
var L5_TOAbands_Dry1 = L5_TOAbands.filter(ee.Filter.dayOfYear(121,180))
    .mean();
var L5_TOAbands_Dry2 = L5_TOAbands.filter(ee.Filter.dayOfYear(241,360))
    .mean();
var L5_TOAbands_Dry = (L5_TOAbands_Dry1.add(L5_TOAbands_Dry2)).divide(2)
    .clip(geometry);

// Make an Array Image, with a 1-D Array per pixel.
var arrayImage1D_Dry = L5_TOAbands_Dry.toArray();

// Make an Array Image with a 2-D Array per pixel, 6x1.
var arrayImage2D_Dry = arrayImage1D_Dry.toArray(1);

// Do a matrix multiplication: 6x6 times 6x1.
var TCT_LS5_Dry = ee.Image(coefficients)
  .matrixMultiply(arrayImage2D_Dry)
  // Get rid of the extra dimensions.
  .arrayProject([0])
  .arrayFlatten(
    [['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth']]);

// Display the first three bands of the result and the input imagery.
var vizTCTdry = {
  bands: ['brightness', 'greenness', 'wetness'],
  min: [0.3677,-0.0493,-0.1178], max: [0.6232,0.0086,0.0088]
};

// Add the resulting layer
Map.addLayer(TCT_LS5_Dry, vizTCTdry, 'TCT_LS5_Dry');

// Export the Multitemporal TCT of the dry months to your Google Drive
Export.image.toDrive({
  image: TCT_LS5_Dry,
  description: 'TCT_L5_Dry',
  scale: 30,
  maxPixels: 1e9,
  region: geometry
});


//                ------------- ooo -------------


// Wet Months Multi-Temporal Principal Component Analysis

// NOTES!
// 1. The bands of interest from the Landsat 5 collection have already been selected in the
//    previous section
// 2. The filter of seasonal images has already been done in previous sections


// Get some information about the input to be used later.
var scale = L5_TOAbands_Wet.projection().nominalScale();
var bandNames = L5_TOAbands_Wet.bandNames();

// Mean center the data to enable a faster covariance reducer
// and an SD stretch of the principal components.
var meanDict = L5_TOAbands_Wet.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: geometry,
    scale: scale,
    maxPixels: 1e9
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = L5_TOAbands_Wet.subtract(means);

// This helper function returns a list of new band names.
var getNewBandNames = function(prefix) {
  var seq = ee.List.sequence(1, bandNames.length());
  return seq.map(function(b) {
    return ee.String(prefix).cat(ee.Number(b).int());
  });
};

// This function accepts mean centered imagery, a scale and
// a region in which to perform the analysis.  It returns the
// Principal Components (PC) in the region as a new image.
var getPrincipalComponents = function(centered, scale, geometry) {
  // Collapse the bands of the image into a 1D array per pixel.
  var arrays = centered.toArray();

  // Compute the covariance of the bands within the region.
  var covar = arrays.reduceRegion({
    reducer: ee.Reducer.centeredCovariance(),
    geometry: geometry,
    scale: scale,
    maxPixels: 1e9
  });

  // Get the 'array' covariance result and cast to an array.
  // This represents the band-to-band covariance within the region.
  var covarArray = ee.Array(covar.get('array'));

  // Perform an eigen analysis and slice apart the values and vectors.
  var eigens = covarArray.eigen();

  // This is a P-length vector of Eigenvalues.
  var eigenValues = eigens.slice(1, 0, 1);
  // This is a PxP matrix with eigenvectors in rows.
  var eigenVectors = eigens.slice(1, 1);

  // Convert the array image to 2D arrays for matrix computations.
  var arrayImage = arrays.toArray(1);

  // Left multiply the image array by the matrix of eigenvectors.
  var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);

  // Turn the square roots of the Eigenvalues into a P-band image.
  var sdImage = ee.Image(eigenValues.sqrt())
    .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);

  // Turn the PCs into a P-band image, normalized by SD.
  return principalComponents
    // Throw out an an unneeded dimension, [[]] -> [].
    .arrayProject([0])
    // Make the one band array image a multi-band image, [] -> image.
    .arrayFlatten([getNewBandNames('pc')])
    // Normalize the PCs by their SDs.
    .divide(sdImage);
};

// Get the PCs at the specified scale and in the specified region
var pcaL5wet = getPrincipalComponents(centered, scale, geometry);

//Define the visualisation parameters
var vizPCAwet = {
  bands: ['pc2', 'pc4', 'pc1'],
  min: [-2.219236753194737], max: [1.6670496438431754]
};

// Display a composite with PCs 2,4,1 for visualisation
// Please change the visualization paramenters to combine different bands
// and/or reproduce the band combinations of the paper's figures
// NOTE! The exported image incorporates all Principal Components
Map.addLayer(pcaL5wet, vizPCAwet, 'PCA_composite_(2,4,1)_Wet');

Export.image.toDrive({
  image: pcaL5wet,
  description: 'PCA_L5_Wet',
  scale: 30,
  maxPixels: 1e9,
  region: geometry
});


//                ------------- ooo -------------


// Dry Months Multi-Temporal Principal Component Analysis

// NOTES!
// 1. The bands of interest from the Landsat 5 collection have already been selected in the
//    previous section
// 2. The filter of seasonal images has already been done in previous sections
// 3. region, scale, and bandNames variables have already been created in the previous section


// Mean center the data to enable a faster covariance reducer
// and an SD stretch of the principal components.
var meanDict_Dry = L5_TOAbands_Dry.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: geometry,
    scale: scale,
    maxPixels: 1e9
});
var means_Dry = ee.Image.constant(meanDict_Dry.values(bandNames));
var centered_Dry = L5_TOAbands_Dry.subtract(means_Dry);

// This helper function returns a list of new band names.
var getNewBandNames_Dry = function(prefix) {
  var seq_Dry = ee.List.sequence(1, bandNames.length());
  return seq_Dry.map(function(b) {
    return ee.String(prefix).cat(ee.Number(b).int());
  });
};

// This function accepts mean centered imagery, a scale and
// a region in which to perform the analysis.  It returns the
// Principal Components (PC) in the region as a new image.
var getPrincipalComponents_Dry = function(centered_Dry, scale, geometry) {
  // Collapse the bands of the image into a 1D array per pixel.
  var arrays_Dry = centered_Dry.toArray();

  // Compute the covariance of the bands within the region.
  var covar_Dry = arrays_Dry.reduceRegion({
    reducer: ee.Reducer.centeredCovariance(),
    geometry: geometry,
    scale: scale,
    maxPixels: 1e9
  });

  // Get the 'array' covariance result and cast to an array.
  // This represents the band-to-band covariance within the region.
  var covarArray_Dry = ee.Array(covar_Dry.get('array'));

  // Perform an eigen analysis and slice apart the values and vectors.
  var eigens_Dry = covarArray_Dry.eigen();

  // This is a P-length vector of Eigenvalues.
  var eigenValues_Dry = eigens_Dry.slice(1, 0, 1);
  // This is a PxP matrix with eigenvectors in rows.
  var eigenVectors_Dry = eigens_Dry.slice(1, 1);

  // Convert the array image to 2D arrays for matrix computations.
  var arrayImage_Dry = arrays_Dry.toArray(1);

  // Left multiply the image array by the matrix of eigenvectors.
  var principalComponents_Dry = ee.Image(eigenVectors_Dry).matrixMultiply(arrayImage_Dry);

  // Turn the square roots of the Eigenvalues into a P-band image.
  var sdImage_Dry = ee.Image(eigenValues_Dry.sqrt())
    .arrayProject([0]).arrayFlatten([getNewBandNames_Dry('sd')]);

  // Turn the PCs into a P-band image, normalized by SD.
  return principalComponents_Dry
    // Throw out an an unneeded dimension, [[]] -> [].
    .arrayProject([0])
    // Make the one band array image a multi-band image, [] -> image.
    .arrayFlatten([getNewBandNames_Dry('pc')])
    // Normalize the PCs by their SDs.
    .divide(sdImage_Dry);
};

// Get the PCs at the specified scale and in the specified region
var pcaL5dry = getPrincipalComponents_Dry(centered_Dry, scale, geometry);

//Define the visualisation parameters
var vizPCAdry = {
  bands: ['pc2', 'pc4', 'pc1'],
  min: [-2.219236753194737], max: [1.6670496438431754]
};

// Display a composite with PCs 2,4,1 for visualisation
// Please change the visualization paramenters to combine different bands
// and/or reproduce the band combinations of the paper's figures
// NOTE! The exported image incorporates all Principal Components
Map.addLayer(pcaL5dry, vizPCAdry, 'PCA_composite_(2,4,1)_Dry');

Export.image.toDrive({
  image: pcaL5dry,
  description: 'PCA_L5_Dry',
  scale: 30,
  maxPixels: 1e9,
  region: geometry
});