using landtrendr to determine duration 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
var viz = {
min: 1,
max: 2017-2001,
palette: ['#FF0000', '#FF7F00', '#FFFF00', '#00FF00', '#0000FF']
};
// run the dist extract function
var distImg = extractDisturbance(lt.select('LandTrendr'), distDir, distParams);
Map.addLayer(distImg.select(['dur']).clip(chl), viz, 'Duration'); // 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: 'Duration',
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);
Hellow, thanks for the post. It really helped me to understand the facilities of this incredible tool that the LT-GEE.
I’m wondering if it’s possible to use ENTINEL image collections, even though the LT-GEE was created, above all, to use Landsat collections.
If you know someone or a group that has done so, can you share names, codes, or pdfs for reading? I would be immensely grateful as I am working with SENTINEL image collection
LikeLike
you could use a similar approach.. though the timeseries is quite short.. you might wanna look into fusion products…
LikeLike
Thanks for sharing your insigths, mister. I’ve trying to use SENTINEL-2 in my studies and I can run a simple script for Disturbance Map, using the guide book from LT_GEE. Heres the code, if you interest: https://code.earthengine.google.com/517007befdbc5bef0146fe44df2e9978
I accept another contribuitions if you gave me the honor.
Thanks for all!
LikeLike