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);
aoiasset. Your script is successful using a different AOI: code.earthengine.google.com/59ff7de6808bfa2bd8cff3473d93908c