2

I am having very basic skills in LandTrendr. I want to make a map showing year of detection and another magnitude of change, but I keep getting errors. I suspect the problem to be with the image collection I have prepared to use.

Here is my area of interest https://drive.google.com/open?id=1F7GNvYpnuEWzvft1JtJgINuagUecvejG

Here is the code:

var coefficients = {
  itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172]).multiply(10000),
  slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071]),
};

// Define function to get and rename bands of interest from OLI.
function renameOLI(img) {
  return img.select(
        ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
          ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']
    );
}

// Define function to get and rename bands of interest from ETM+.
function renameETM(img) {
  return img.select(
        ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
        ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']
  );
}


// Define function to apply harmonization transformation.
function etm2oli(img) {
  return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
    .multiply(coefficients.slopes)
    .add(coefficients.itcps)
    .round()
    .toShort()
    .addBands(img.select('pixel_qa')
    .copyProperties(img, ['system:time_start'])
  );
}

// Define function to mask out clouds and cloud shadows.
function fmask(img) {
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;
  var qa = img.select('pixel_qa');
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return img.updateMask(mask);
}


// Define function to prepare OLI images.
function prepOLI(img) {
  var orig = img;
  img = renameOLI(img);
  img = fmask(img);
    return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}

// Define function to prepare ETM+ images.
function prepETM(img) {
  var orig = img;
  img = renameETM(img);
  img = fmask(img);
  img = etm2oli(img);
  return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}



// Define AOI on the map.
Map.centerObject(aoi, 10);
Map.addLayer(aoi, {color: 'f8766d'}, 'AOI');
Map.setOptions('HYBRID');

// Get Landsat surface reflectance collections for OLI, ETM+ and TM sensors.
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol= ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');

// Define a collection filter.
var colFilter = ee.Filter.and(
  ee.Filter.bounds(aoi),
  ee.Filter.calendarRange(1, 365, 'day_of_year'),
  ee.Filter.lt('CLOUD_COVER', 50),
  ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
  ee.Filter.or(
    ee.Filter.eq('IMAGE_QUALITY', 9),
    ee.Filter.eq('IMAGE_QUALITY_OLI', 9)
  )
);

// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOLI);
etmCol= etmCol.filter(colFilter).map(prepETM);


// Merge the collections.
var col = oliCol
  .merge(etmCol);



// Define start and end years.
var startYear = 2000;
var endYear = 2018;

// Define a reducer.
var myReducer = ee.Reducer.mean();



// Make a list of years to generate composites for.
var yearList = ee.List.sequence(startYear, endYear);


// Map over the list of years to generate a composite for each year.
var yearCompList = yearList.map(function(year){
  var yearCol = col.filter(ee.Filter.calendarRange(year, year, 'year'));

   var yearComp = yearCol.reduce(myReducer);
  var imgList = yearCol.aggregate_array('constant');
  var systemStart = yearCol.reduceColumns(ee.Reducer.min(), ['system:time_start']).get('min');
  // Reduce (composite) the images for this year.
  var nBands = yearComp.bandNames().size();
 return yearComp.set({
    'year': year,
    'image_list': imgList,
    'n_bands': nBands,
    'system:time_start': systemStart
  });
});

//print("image",yearCompList)
// Convert the annual composite image list to an ImageCollection
var yearCompCol = ee.ImageCollection.fromImages(yearCompList);

// Filter out years with no bands.
//(can happen if there were no images to composite)
yearCompCol = yearCompCol.filter(ee.Filter.gt('n_bands', 0));

print("image",yearCompCol);

// define function to calculate a spectral index to segment with LT
var Ndvi = function(img) {
    var index = img.normalizedDifference(['NIR_mean', 'Red_mean'])
                   .select([0], ['NDVI'])
                   .multiply(1000)
                   .set('system:time_start', img.get('system:time_start'));
    return img.addBands(index) ;

};

var distDir = -1; // define the sign of spectral delta for vegetation loss for the segmentation index - 
                 


// define the segmentation parameters:
// reference: Kennedy, R. E., Yang, Z., & Cohen, W. B. (2010). Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sensing of Environment, 114(12), 2897-2910.
//            https://github.com/eMapR/LT-GEE
var run_params = { 
  maxSegments:            6,
  spikeThreshold:         0.9,
  vertexCountOvershoot:   3,
  preventOneYearRecovery: true,
  recoveryThreshold:      0.25,
  pvalThreshold:          0.05,
  bestModelProportion:    0.75,
  minObservationsNeeded:  6
};


// apply the function to calculate the segmentation index and adjust the values by the distDir parameter - flip index so that a vegetation loss is associated with a postive delta in spectral value
var ltCollection = yearCompCol.map(Ndvi)                                             // map the function over every image in the collection - returns a 1-band annual image collection of the spectral index

ltCollection = ltCollection.map(function(img) {return img.select("NDVI").multiply(distDir)           // ...multiply the segmentation index by the distDir to ensure that vegetation loss is associated with a positive spectral delta
                           .set('system:time_start', img.get('system:time_start'))});

//----- RUN LANDTRENDR -----
run_params.timeSeries = ltCollection;               // add LT collection to the segmentation run parameter object
var lt = ee.Algorithms.TemporalSegmentation.LandTrendr(run_params); // run LandTrendr spectral temporal segmentation algorithm

// define disturbance mapping filter parameters 
var treeLoss1  = 175;      // delta filter for 1 year duration disturbance, <= will not be included as disturbance - units are in units of segIndex defined in the following function definition
var treeLoss20 = 200;      // delta filter for 20 year duration disturbance, <= will not be included as disturbance - units are in units of segIndex defined in the following function definition
var preVal     = 400;      // pre-disturbance value threshold - values below the provided threshold will exclude disturbance for those pixels - units are in units of segIndex defined in the following function definition
var mmu        = 15;       // minimum mapping unit for disturbance patches - units of pixels

// assemble the disturbance extraction parameters
var distParams = {
  tree_loss1: treeLoss1,
  tree_loss20: treeLoss20,  
  pre_val: preVal           
};


// ----- function to extract greatest disturbance based on spectral delta between vertices 
var extractDisturbance = function(lt, distDir, params, mmu) {
  // select only the vertices that represents a change
  var vertexMask = lt.arraySlice(0, 3, 4); // get the vertex - yes(1)/no(0) dimension
  var vertices = lt.arrayMask(vertexMask); // convert the 0's to masked
  
  // construct segment start and end point years and index values
  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 

  // concatenate segment start year, delta, duration, and starting spectral index value to an array 
  var distImg = ee.Image.cat([startYear.add(1), mag, dur, startVal.multiply(distDir)]).toArray(0); // make an image of segment attributes - multiply by the distDir parameter to re-orient the spectral index if it was flipped for segmentation - do it here so that the subtraction to calculate segment delta in the above line is consistent - add 1 to the detection year, because the vertex year is not the first year that change is detected, it is the following year
 
  // sort the segments in the disturbance attribute image delta by spectral index change delta  
  var distImgSorted = distImg.arraySort(mag.multiply(-1));                                  // flip the delta around so that the greatest delta segment is first in order

  // slice out the first (greatest) delta
  var tempDistImg = distImgSorted.arraySlice(1, 0, 1).unmask(ee.Image(ee.Array([[0],[0],[0],[0]])));                                      // get the first segment in the sorted array

  // make an image from the array of attributes for the greatest disturbance
  var finalDistImg = ee.Image.cat(tempDistImg.arraySlice(0,0,1).arrayProject([1]).arrayFlatten([['yod']]),     // slice out year of disturbance detection and re-arrange to an image band 
                                  tempDistImg.arraySlice(0,1,2).arrayProject([1]).arrayFlatten([['mag']]),     // slice out the disturbance magnitude and re-arrange to an image band 
                                  tempDistImg.arraySlice(0,2,3).arrayProject([1]).arrayFlatten([['dur']]),     // slice out the disturbance duration and re-arrange to an image band
                                  tempDistImg.arraySlice(0,3,4).arrayProject([1]).arrayFlatten([['preval']])); // slice out the pre-disturbance spectral value and re-arrange to an image band
  
  // filter out disturbances based on user settings
  var threshold = ee.Image(finalDistImg.select(['dur']))                        // get the disturbance band out to apply duration dynamic disturbance magnitude threshold 
                    .multiply((params.tree_loss20 - params.tree_loss1) / 19.0)  // ...
                    .add(params.tree_loss1)                                     //    ...interpolate the magnitude threshold over years between a 1-year mag thresh and a 20-year mag thresh
                    .lte(finalDistImg.select(['mag']))                          // ...is disturbance less then equal to the interpolated, duration dynamic disturbance magnitude threshold 
                    .and(finalDistImg.select(['mag']).gt(0))                    // and is greater than 0  
                    .and(finalDistImg.select(['preval']).gt(params.pre_val));   // and is greater than pre-disturbance spectral index value threshold
  
  // apply the filter mask
  finalDistImg = finalDistImg.mask(threshold).int16(); 
  
   // patchify the remaining disturbance pixels using a minimum mapping unit
  if(mmu > 1){
    var mmuPatches = finalDistImg.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(mmu);                      // are the the number of pixels per candidate patch greater than user-defined minimum mapping unit?
    finalDistImg = finalDistImg.updateMask(mmuPatches);     // mask the pixels/patches that are less than minimum mapping unit
  } 
  
  return finalDistImg; // return the filtered greatest disturbance attribute image
};

var viz = {
  min: 2001,
  max: 2017,
  palette: ['#9400D3', '#4B0082', '#0000FF', '#00FF00', '#FFFF00', '#FF7F00', '#FF0000']
};

// run the dist extract function
var distImg = extractDisturbance(lt.select('LandTrendr'), distDir, distParams);

Map.addLayer(distImg.select(['yod']).clip(aoi), viz, 'Year of Detection');    // add disturbance year of detection to map

// set position of panel
var legend = ui.Panel({
style: {
position: 'bottom-left',
padding: '8px 15px'
}
});
 
// Create legend title
var legendTitle = ui.Label({
value: 'Year',
style: {
fontWeight: 'bold',
fontSize: '18px',
margin: '0 0 4px 0',
padding: '0'
}
});
 
// Add the title to the panel
legend.add(legendTitle);
 
// create the legend image
var lon = ee.Image.pixelLonLat().select('latitude');
var gradient = lon.multiply((viz.max-viz.min)/100.0).add(viz.min);
var legendImage = gradient.visualize(viz);
 
// create text on top of legend
var panel = ui.Panel({
widgets: [
ui.Label(viz['max'])
],
});
 
legend.add(panel);
 
// create thumbnail from the image
var thumbnail = ui.Thumbnail({
image: legendImage,
params: {bbox:'0,0,10,100', dimensions:'10x200'},
style: {padding: '1px', position: 'bottom-center'}
});
 
// add the thumbnail to the legend
legend.add(thumbnail);
 
// create text on top of legend
var panel = ui.Panel({
widgets: [
ui.Label(viz['min'])
],
});
 
legend.add(panel);
 
Map.add(legend);
2

1 Answer 1

1
var aoi = 
    /* color: #d63000 */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[-47.05253256069382, -10.953258609139985],
          [-47.05253256069382, -11.025045428557927],
          [-46.98352468715866, -11.025045428557927],
          [-46.98352468715866, -10.953258609139985]]], null, false);


var coefficients = {
  itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172]).multiply(10000),
  slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071]),
};


// Define start and end years.
var startYear = 2009;
var endYear = 2019;

var viz = {
  min: 2009,
  max: 2019,
  palette: ['#9400D3', '#4B0082', '#0000FF', '#00FF00', '#FFFF00', '#FF7F00', '#FF0000']
};

// define disturbance mapping filter parameters 
var treeLoss1  = 175;      // delta filter for 1 year duration disturbance, <= will not be included as disturbance - units are in units of segIndex defined in the following function definition
var treeLoss20 = 80;//80;      // delta filter for 20 year duration disturbance, <= will not be included as disturbance - units are in units of segIndex defined in the following function definition
var preVal     = 400;      // pre-disturbance value threshold - values below the provided threshold will exclude disturbance for those pixels - units are in units of segIndex defined in the following function definition
var mmu        = 1;       // minimum mapping unit for disturbance patches - units of pixels

var vizParams = {
  bands: ['b1', 'b2', 'b3'],
  min: 0,
  max: 256,
  //gamma: [0.95, 1.1, 1]
};


var MNN_1989 = ee.Image("users/maximoicmbio/MNN_19890609");
var MNN_1999 = ee.Image("users/maximoicmbio/MNN_19990605");
var MNN_2009 = ee.Image("users/maximoicmbio/MNN_20090702");
var MNN_2019 = ee.Image("users/maximoicmbio/MNN_20190612");
print(MNN_2019)
Map.addLayer(MNN_2019, vizParams, 'MNN_2019');


// Define function to get and rename bands of interest from OLI.
function renameOLI(img) {
  return img.select(
        ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
          ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']
    );
}

// Define function to get and rename bands of interest from ETM+.
function renameETM(img) {
  return img.select(
        ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
        ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']
  );
}

// Define function to get and rename bands of interest from TM.
function renameTM(img) {
  return img.select(
        ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
        ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']
  );
}

// Define function to apply harmonization transformation.
function etm2oli(img) {
  return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
    .multiply(coefficients.slopes)
    .add(coefficients.itcps)
    .round()
    .toShort()
    .addBands(img.select('pixel_qa')
    .copyProperties(img, ['system:time_start'])
  );
}

function tm2oli(img) {
  return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
    .multiply(coefficients.slopes)
    .add(coefficients.itcps)
    .round()
    .toShort()
    .addBands(img.select('pixel_qa')
    .copyProperties(img, ['system:time_start'])
  );
}
// Define function to mask out clouds and cloud shadows.
/*function fmask(img) {
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;
  var qa = img.select('pixel_qa');
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return img.updateMask(mask);
}
*/
// landsat 5 and landsat 7
var cloudMaskL457 = function(img) {
  var qa = img.select('pixel_qa');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));
  // Remove edge pixels that don't occur in all bands
  var mask2 = img.mask().reduce(ee.Reducer.min());
  return img.updateMask(cloud.not()).updateMask(mask2);
};

// landsat 8
function maskL8sr(img) {
  // 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 = img.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 img.updateMask(mask);
}

// Define function to prepare OLI images.
function prepOLI(img) {
  var orig = img;
  img = renameOLI(img);
  img = maskL8sr(img);
    return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}

// Define function to prepare ETM+ images.
function prepETM(img) {
  var orig = img;
  img = renameETM(img);
  img = cloudMaskL457(img);
  img = etm2oli(img);
  return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}

function prepTM(img) {
  var orig = img;
  img = renameTM(img);
  img = cloudMaskL457(img);
  img = tm2oli(img);
  return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}

// Define AOI on the map.
Map.centerObject(aoi, 12);
//Map.addLayer(aoi, {color: 'f8766d'}, 'AOI');
Map.setOptions('HYBRID');


// Get Landsat surface reflectance collections for OLI, ETM+ and TM sensors.
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol= ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
var tmCol= ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');

// Define a collection filter.
var colFilter = ee.Filter.and(
  ee.Filter.bounds(aoi),
  ee.Filter.calendarRange(1, 365, 'day_of_year'),
  ee.Filter.lt('CLOUD_COVER', 50),
  ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
  ee.Filter.or(
    ee.Filter.eq('IMAGE_QUALITY', 9),
    ee.Filter.eq('IMAGE_QUALITY_OLI', 9)
  )
);

// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOLI);
etmCol= etmCol.filter(colFilter).map(prepETM);
tmCol= tmCol.filter(colFilter).map(prepTM);

// Merge the collections.
var col = oliCol.merge(etmCol).merge(tmCol);


// Define a reducer.

var myReducer = ee.Reducer.mean();
//var myReducer = ee.Reducer.median();
//var myReducer = ee.Reducer.min();
//var myReducer = ee.Reducer.max();

// Make a list of years to generate composites for.
var yearList = ee.List.sequence(startYear, endYear);


// Map over the list of years to generate a composite for each year.
var yearCompList = yearList.map(function(year){
  var yearCol = col.filter(ee.Filter.calendarRange(year, year, 'year'));

   var yearComp = yearCol.reduce(myReducer);
  var imgList = yearCol.aggregate_array('constant');
  var systemStart = yearCol.reduceColumns(ee.Reducer.min(), ['system:time_start']).get('min');
  // Reduce (composite) the images for this year.
  var nBands = yearComp.bandNames().size();
 return yearComp.set({
    'year': year,
    'image_list': imgList,
    'n_bands': nBands,
    'system:time_start': systemStart
  });
});

//print("image",yearCompList)
// Convert the annual composite image list to an ImageCollection
var yearCompCol = ee.ImageCollection.fromImages(yearCompList);

// Filter out years with no bands.
//(can happen if there were no images to composite)
yearCompCol = yearCompCol.filter(ee.Filter.gt('n_bands', 0));

print("image",yearCompCol);

// define function to calculate a spectral index to segment with LT
var Ndvi = function(img) {
    var index = img.normalizedDifference(['NIR_mean', 'Red_mean'])
                   .select([0], ['NDVI'])
                   .multiply(1000)
                   .set('system:time_start', img.get('system:time_start'));
    return img.addBands(index) ;

};


// NBR
/*var Ndvi = function(img) {

    var index = img.normalizedDifference(['NIR_mean', 'SWIR1_mean'])                
                   .select([0], ['NDVI'])//NBR
                   .multiply(1000)  
                   .set('system:time_start', img.get('system:time_start')); 
    return img.addBands(index) ;
}
*/

var distDir = -1; // define the sign of spectral delta for vegetation loss for the segmentation index - 



// define the segmentation parameters:
// reference: Kennedy, R. E., Yang, Z., & Cohen, W. B. (2010). Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sensing of Environment, 114(12), 2897-2910.
//            https://github.com/eMapR/LT-GEE
var run_params = { 
  maxSegments:            6,
  spikeThreshold:         0.9,
  vertexCountOvershoot:   3,
  preventOneYearRecovery: true,
  recoveryThreshold:      0.25,
  pvalThreshold:          0.05,
  bestModelProportion:    0.75,
  minObservationsNeeded:  6
};


// apply the function to calculate the segmentation index and adjust the values by the distDir parameter - flip index so that a vegetation loss is associated with a postive delta in spectral value
var ltCollection = yearCompCol.map(Ndvi)                                             // map the function over every image in the collection - returns a 1-band annual image collection of the spectral index

ltCollection = ltCollection.map(function(img) {return img.select("NDVI").multiply(distDir)           // ...multiply the segmentation index by the distDir to ensure that vegetation loss is associated with a positive spectral delta
                           .set('system:time_start', img.get('system:time_start'))});

//----- RUN LANDTRENDR -----
run_params.timeSeries = ltCollection;               // add LT collection to the segmentation run parameter object
var lt = ee.Algorithms.TemporalSegmentation.LandTrendr(run_params); // run LandTrendr spectral temporal segmentation algorithm

// assemble the disturbance extraction parameters
var distParams = {
  tree_loss1: treeLoss1,
  tree_loss20: treeLoss20,  
  pre_val: preVal           
};


// ----- function to extract greatest disturbance based on spectral delta between vertices 
var extractDisturbance = function(lt, distDir, params, mmu) {
  // select only the vertices that represents a change
  var vertexMask = lt.arraySlice(0, 3, 4); // get the vertex - yes(1)/no(0) dimension
  var vertices = lt.arrayMask(vertexMask); // convert the 0's to masked

  // construct segment start and end point years and index values
  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 

  // concatenate segment start year, delta, duration, and starting spectral index value to an array 
  var distImg = ee.Image.cat([startYear.add(1), mag, dur, startVal.multiply(distDir)]).toArray(0); // make an image of segment attributes - multiply by the distDir parameter to re-orient the spectral index if it was flipped for segmentation - do it here so that the subtraction to calculate segment delta in the above line is consistent - add 1 to the detection year, because the vertex year is not the first year that change is detected, it is the following year

  // sort the segments in the disturbance attribute image delta by spectral index change delta  
  var distImgSorted = distImg.arraySort(mag.multiply(-1));                                  // flip the delta around so that the greatest delta segment is first in order

  // slice out the first (greatest) delta
  var tempDistImg = distImgSorted.arraySlice(1, 0, 1).unmask(ee.Image(ee.Array([[0],[0],[0],[0]])));                                      // get the first segment in the sorted array

  // make an image from the array of attributes for the greatest disturbance
  var finalDistImg = ee.Image.cat(tempDistImg.arraySlice(0,0,1).arrayProject([1]).arrayFlatten([['yod']]),     // slice out year of disturbance detection and re-arrange to an image band 
                                  tempDistImg.arraySlice(0,1,2).arrayProject([1]).arrayFlatten([['mag']]),     // slice out the disturbance magnitude and re-arrange to an image band 
                                  tempDistImg.arraySlice(0,2,3).arrayProject([1]).arrayFlatten([['dur']]),     // slice out the disturbance duration and re-arrange to an image band
                                  tempDistImg.arraySlice(0,3,4).arrayProject([1]).arrayFlatten([['preval']])); // slice out the pre-disturbance spectral value and re-arrange to an image band

  // filter out disturbances based on user settings
  var threshold = ee.Image(finalDistImg.select(['dur']))                        // get the disturbance band out to apply duration dynamic disturbance magnitude threshold 
                    .multiply((params.tree_loss20 - params.tree_loss1) / 19.0)  // ...
                    .add(params.tree_loss1)                                     //    ...interpolate the magnitude threshold over years between a 1-year mag thresh and a 20-year mag thresh
                    .lte(finalDistImg.select(['mag']))                          // ...is disturbance less then equal to the interpolated, duration dynamic disturbance magnitude threshold 
                    .and(finalDistImg.select(['mag']).gt(0))                    // and is greater than 0  
                    .and(finalDistImg.select(['preval']).gt(params.pre_val));   // and is greater than pre-disturbance spectral index value threshold

  // apply the filter mask
  finalDistImg = finalDistImg.mask(threshold).int16(); 

   // patchify the remaining disturbance pixels using a minimum mapping unit
  if(mmu > 1){
    var mmuPatches = finalDistImg.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(mmu);                      // are the the number of pixels per candidate patch greater than user-defined minimum mapping unit?
    finalDistImg = finalDistImg.updateMask(mmuPatches);     // mask the pixels/patches that are less than minimum mapping unit
  } 

  return finalDistImg; // return the filtered greatest disturbance attribute image
};


// run the dist extract function

var exportImg = extractDisturbance(lt.select('LandTrendr'), distDir, distParams).select(['yod']).clip(aoi).unmask(0).short();
Export.image.toDrive({
  image: exportImg, 
  description: 'lt-gee_disturbance_map', 
  folder: 'lt-gee_disturbance_map', 
  fileNamePrefix: 'lt-gee_disturbance_map', 
  region: aoi, 
  scale: 30, 
  crs: 'EPSG:32723', 
  maxPixels: 1e13
});

var distImg = extractDisturbance(lt.select('LandTrendr'), distDir, distParams).select(['yod']).clip(aoi);

Map.addLayer(distImg, viz, 'Year_of_Detection');    // add disturbance year of detection to map

// set position of panel
var legend = ui.Panel({
style: {
position: 'bottom-left',
padding: '8px 15px'
}
});

// Create legend title
var legendTitle = ui.Label({
value: 'Year',
style: {
fontWeight: 'bold',
fontSize: '18px',
margin: '0 0 4px 0',
padding: '0'
}
});

// Add the title to the panel
legend.add(legendTitle);

// create the legend image
var lon = ee.Image.pixelLonLat().select('latitude');
var gradient = lon.multiply((viz.max-viz.min)/100.0).add(viz.min);
var legendImage = gradient.visualize(viz);

// create text on top of legend
var panel = ui.Panel({
widgets: [
ui.Label(viz['max'])
],
});

legend.add(panel);

// create thumbnail from the image
var thumbnail = ui.Thumbnail({
image: legendImage,
params: {bbox:'0,0,10,100', dimensions:'10x200'},
style: {padding: '1px', position: 'bottom-center'}
});

// add the thumbnail to the legend
legend.add(thumbnail);

// create text on top of legend
var panel = ui.Panel({
widgets: [
ui.Label(viz['min'])
],
});

legend.add(panel);

Map.add(legend);

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.