using landtrendr to determine magnitude 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: distParams.tree_loss1, max: 1000, palette: ['#0000FF', '#00FF00', '#FFFF00', '#FF7F00', '#FF0000'] }; // run the dist extract function var distImg = extractDisturbance(lt.select('LandTrendr'), distDir, distParams); Map.addLayer(distImg.select(['mag']).clip(chl), viz, 'Magnitude'); // 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: 'Magnitude', 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)
Thank you so much for your sharing and it does help to me, a rookie GEEER. When I try your code, there is a warning:
Magnitude: Tile error: Image must have a numeric start time property ‘system:time_start’. I don’t know what it means ..
LikeLike
How did you prepare the annual collection? was it harmoninized firstly?
LikeLike
Hello…
What would be the units for magnitude of change?
I got the values for gains between 300 to 700
What do they really indicate?
LikeLike
what is ‘distDir’? Really hard to understand the role of this object from the comment only.
LikeLike
this variable describes the direction of change. -1 for regressive and 1 for gain/recovery
LikeLike
that means there is no date attached to the image..
LikeLike