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