change detection 1: year of detection

using landtrendr to determine year of change

LandTrendr (Landsat-based detection of Trends in Disturbance and Recovery) is an approach to extract spectral trajectories of land surface change from yearly Landsat time-series stacks (LTS) . The method brings together two themes in time-series analysis of LTS: capture of short-duration events and smoothing of long-term trends.

https://github.com/eMapR/LT-GEE

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.

Step 1: get the collection and area of interest.

// get annual collection
var annualSRcollection = ee.ImageCollection("projects/servir-mekong/yearlyComposites")

// set the area of interest
var chl = ee.FeatureCollection("users/servirmekong/Vietnam/CHL_Boundary").geometry();

Step 2: add the function to calculate the nbr

// define function to calculate a spectral index to segment with LT
var nbr = function(img) {

    var index = img.normalizedDifference(['nir', 'swir1'])
                   .select([0], ['NBR'])
                   .multiply(1000)
                   .set('system:time_start', img.get('system:time_start'));
    return img.addBands(index) ;
}

Step 3: apply the function to the image collection

var ltCollection = annualSRcollection.map(nbr)

Step 4: Define the segmentation parameters

// define the segmentation parameters:
var run_params = {
  maxSegments:            6,
  spikeThreshold:         0.9,
  vertexCountOvershoot:   3,
  preventOneYearRecovery: true,
  recoveryThreshold:      0.25,
  pvalThreshold:          0.05,
  bestModelProportion:    0.75,
  minObservationsNeeded:  6
};

Step 5: Multiply the segmentation index by the distDir to ensure that vegetation loss is associated with a positive spectral delta

ltCollection = ltCollection.map(function(img) {
return img.select("NBR").multiply(distDir)
                 .set('system:time_start', img.get('system:time_start'))

Step 6: Run the landtrendr algorithm

//----- RUN LANDTRENDR -----
run_params.timeSeries = ltCollection;
var lt = ee.Algorithms.TemporalSegmentation.LandTrendr(run_params);

Step 7: define disturbance mapping filter parameters

// define disturbance mapping filter parameters
 // 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 treeLoss1  = 175;
// 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 treeLoss20 = 200;
// 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 preVal     = 400;
 // minimum mapping unit for disturbance patches - units of pixels
var mmu        = 15;      

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

Step 7: add function to extract greatest disturbance based on spectral delta between vertices.

// ----- 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);
var vertices = lt.arrayMask(vertexMask);

// construct segment start and end point years and index values
var left = vertices.arraySlice(1, 0, -1);
var right = vertices.arraySlice(1, 1, null);
var startYear = left.arraySlice(0, 0, 1);
var startVal = left.arraySlice(0, 2, 3);
var endYear = right.arraySlice(0, 0, 1);
var endVal = right.arraySlice(0, 2, 3);

var dur = endYear.subtract(startYear);
var mag = endVal.subtract(startVal);

// 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);

// sort the segments in the disturbance attribute image delta by spectral index change delta
var distImgSorted = distImg.arraySort(mag.multiply(-1));

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

// 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']]),
tempDistImg.arraySlice(0,1,2).arrayProject([1]).arrayFlatten([['mag']]),
tempDistImg.arraySlice(0,2,3).arrayProject([1]).arrayFlatten([['dur']]),
tempDistImg.arraySlice(0,3,4).arrayProject([1]).arrayFlatten([['preval']]));

// filter out disturbances based on user settings
var threshold = ee.Image(finalDistImg.select(['dur']))
.multiply((params.tree_loss20 - params.tree_loss1) / 19.0)
.add(params.tree_loss1)
.lte(finalDistImg.select(['mag']))
.and(finalDistImg.select(['mag']).gt(0))
.and(finalDistImg.select(['preval']).gt(params.pre_val));

// 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'])
.connectedPixelCount(mmu, true)
.gte(mmu);
finalDistImg = finalDistImg.updateMask(mmuPatches);
}

return finalDistImg;
};

Step 8: Run the extraction function

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

Step 8: add the map and legend

Map.addLayer(distImg.select(['yod']).clip(chl), 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);

example