How to get the fitted data in GEE LandTrendr for the changes defined by vertices of initial veg. index (e.g. NBR)?

by Maksym Matsala   Last Updated October 19, 2019 15:22 PM

What I have: NBR layes of magnitude, year of disturbance, etc. for the AOI.

What I expect to have: fitted values of TCW, TCB, TCG defined by vertices of NBR time series I mentioned above.

I cannot: just create new layers of TCW or other index/band, since these disturbed pixels will be different from the ones created by NBR.

I want to: not parametirize filters for TCW, TCB or TCG, I want to have just fitted data for the pixels created by NBR.

Why I need that: to have a set of magnitude and disturbance rate of NBR and tasseled cap transformation for the pixels defined by NBR, and then use it as the predictor variables for my RF model (to classify disturbance causal agents).

The code of NBR layers below is typical and obtained from Robert Kennedy eMapR guide. Can I get from some of these variables fitted data for the TCT?

var boundary = ee.FeatureCollection(table);
var region = boundary.geometry();
var ltgee = require('users/emaprlab/public:Modules/LandTrendr.js');

var distDir = -1;
var mmu = 11;

var startYear = 1986;
var endYear = 2018;
var startDay= '06-01';
var endDay = '08-30';
var aoi = region;
var index = 'NBR';
var ftvList = ['TCB', 'TCG', 'TCW', 'NBR'];
var maskThese = ['cloud', 'shadow', 'snow', 'water'];

var runParams = {
  maxSegments: 6,
  spikeThreshold: 0.9,
  vertexCountOvershoot: 3,
  preventOneYearRecovery: true,
  recoveryThreshold: 0.25,
  pvalThreshold: 0.05,
  bestModelProportion: 0.75,
  minObservationsNeeded: 6
};

var lt = ltgee.runLT(startYear, endYear, startDay, endDay, aoi, index, ftvList, runParams, maskThese)
.select('LandTrendr');

// Getting segment information:
var vertexMask = lt.arraySlice(0, 3, 4); // slice out the 'Is Vertex' row - yes(1)/no(0)
var vertices = lt.arrayMask(vertexMask); // use the 'Is Vertex' row as a mask for all rows
var left = vertices.arraySlice(1, 0, -1);    // slice out the vertices as the start of segments
var right = vertices.arraySlice(1, 1, null); // slice out the vertices as the end of segments
var startYear = left.arraySlice(0, 0, 1);    // get year dimension of LT data from the segment start vertices
var startVal = left.arraySlice(0, 2, 3);     // get spectral index dimension of LT data from the segment start vertices
var endYear = right.arraySlice(0, 0, 1);     // get year dimension of LT data from the segment end vertices 
var endVal = right.arraySlice(0, 2, 3);      // get spectral index dimension of LT data from the segment end vertices

var dur = endYear.subtract(startYear);       // subtract the segment start year from the segment end year to calculate the duration of segments 
var mag = endVal.subtract(startVal);         // substract the segment start index value from the segment end index value to calculate the delta of segments
var rate = mag.divide(dur);                  // calculate the rate of spectral change

var segInfo = ee.Image.cat([startYear.add(1), endYear, startVal, endVal, mag, dur, rate])
                      .toArray(0)
                      .mask(vertexMask.mask());

// Isolate segment of interest:
var sortByThis = segInfo.arraySlice(0, 4, 5).toArray(0).multiply(-1); // need to flip the delta here, since arraySort is working by ascending order
var segInfoSorted = segInfo.arraySort(sortByThis); // sort the array by magnitude
var bigDelta = segInfoSorted.arraySlice(1, 0, 1); // get the first segment in the sorted array (greatest magnitude vegetation loss segment)

//
var bigDeltaImg = ee.Image.cat(bigDelta.arraySlice(0,0,1).arrayProject([1]).arrayFlatten([['yod']]),
                               bigDelta.arraySlice(0,1,2).arrayProject([1]).arrayFlatten([['endYr']]),
                               bigDelta.arraySlice(0,2,3).arrayProject([1]).arrayFlatten([['startVal']]),//.multiply(distDir),
                               bigDelta.arraySlice(0,3,4).arrayProject([1]).arrayFlatten([['endVal']]),//.multiply(distDir),
                               bigDelta.arraySlice(0,4,5).arrayProject([1]).arrayFlatten([['mag']]).multiply(distDir),
                               bigDelta.arraySlice(0,5,6).arrayProject([1]).arrayFlatten([['dur']]),
                               bigDelta.arraySlice(0,6,7).arrayProject([1]).arrayFlatten([['rate']]).multiply(distDir));

var distMask =  bigDeltaImg.select(['mag']).lt(-200)
                                           .and(bigDeltaImg.select(['dur']).lt(4))
                                           .and(bigDeltaImg.select(['startVal']).lt(-300));
var bigFastDist = bigDeltaImg.mask(distMask).int16(); // need to set as int16 bit to use connectedPixelCount for minimum mapping unit filter

var mmuPatches = bigFastDist.select(['yod'])                // patchify based on disturbances having the same year of detection
                            .connectedPixelCount(mmu, true) // count the number of pixel in a candidate patch
                            .gte(11);                       // are the the number of pixels per candidate patch greater than user-defined minimum mapping unit?
var bigFastDist = bigFastDist.updateMask(mmuPatches);       // mask the pixels/patches that are less than minimum mapping unit

// Mapping the greatest vegetation
var exportImg = bigFastDist.clip(region).unmask(0).short();
//Map.addLayer(exportImg);
Export.image.toDrive({
  image: exportImg, 
  description: 'NBR_new', 
  folder: 'NBR_new', 
  fileNamePrefix: 'NBR_new', 
  region: region, 
  scale: 30, 
  crs: 'EPSG:4326', 
  maxPixels: 1e13
});


Related Questions


Updated October 16, 2019 16:22 PM

Updated November 13, 2018 15:22 PM

Updated May 22, 2018 16:22 PM

Updated October 16, 2019 01:22 AM