GEE：LandTrendr时间序列曲线拟合
作者：互联网
毕竟是时间序列分析，LandTrendr的时间序列曲线是怎么拟合的？在一个像素点的位置上，它是怎么输出变化年份、扰动幅度、恢复幅度和时间的？
以下代码给出了回答，可以参考LT-GEE官方的说明文档。
用LandTrendr选取指定像素点的APP：https://emaprlab.users.earthengine.app/view/lt-gee-pixel-time-series
LT-GEE的结果包括：
每个像素时间序列的观测年份；二维频谱时间空间中的x轴值；（默认）
每个像素时间序列的观测值的源值；二维频谱时间空间中的y轴值；（默认）
拟合到每个像素时间序列的顶点之间的线段（FTV）的观测值的源值；二维频谱时间空间中的y轴值；（默认）
FTV值相对于源值的均方根误差（RMSE）；（默认）
集合中大于波段1的其他波段的完整时间序列FTV值；二维频谱时间空间中的y轴值；（可选的）
时间序列曲线拟合，结果如下图所示：
源代码：（复制粘贴到GEE里面就可以直接运行了）https://code.earthengine.google.com/12a8d012f2a566995a7c1de1c257b479?noload=true
//# LANDTRENDR SOURCE AND FITTING PIXEL TIME SERIES
// date: 2018-04-19
// author: Zhiqiang Yang | zhiqiang.yang@oregonstate.edu
// Justin Braaten | jstnbraaten@gmail.com
// Robert Kennedy | rkennedy@coas.oregonstate.edu
// website: https://github.com/eMapR/LT-GEE
// 首先输入一个点的经纬度
var long = -122.8848;
var lat = 43.7929;
// 定义数据集参数
var startYear = 1985; // 时间序列开始年份
var endYear = 2017; // 时间序列结束年份
var startDay = '06-01'; // 筛选每年的开始时间
var endDay = '09-30'; // 筛选每年的结束时间
// 计算光谱指数（NDVI等）
var segIndex = function(img) {
var index = img.normalizedDifference(['B4', 'B7'])
.multiply(1000)
.select([0], ['NBR'])
.set('system:time_start', img.get('system:time_start'));
return index ;
};
var distDir = -1; // define the sign of spectral delta for vegetation loss for the segmentation index -
// NBR delta is negetive for vegetation loss, so -1 for NBR, 1 for band 5, -1 for NDVI, etc
// 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
};
//################################################################################################
//########################################################################################################
//##### ANNUAL SR TIME SERIES COLLECTION BUILDING FUNCTIONS #####
//########################################################################################################
//----- MAKE A DUMMY COLLECTOIN FOR FILLTING MISSING YEARS -----
var dummyCollection = ee.ImageCollection([ee.Image([0,0,0,0,0,0]).mask(ee.Image(0))]); // make an image collection from an image with 6 bands all set to 0 and then make them masked values
//------ L8 to L7 HARMONIZATION FUNCTION -----
// slope and intercept citation: Roy, D.P., Kovalskyy, V., Zhang, H.K., Vermote, E.F., Yan, L., Kumar, S.S, Egorov, A., 2016, Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity, Remote Sensing of Environment, 185, 57-70.(http://dx.doi.org/10.1016/j.rse.2015.12.024); Table 2 - reduced major axis (RMA) regression coefficients
var harmonizationRoy = function(oli) {
var slopes = ee.Image.constant([0.9785, 0.9542, 0.9825, 1.0073, 1.0171, 0.9949]); // create an image of slopes per band for L8 TO L7 regression line - David Roy
var itcp = ee.Image.constant([-0.0095, -0.0016, -0.0022, -0.0021, -0.0030, 0.0029]); // create an image of y-intercepts per band for L8 TO L7 regression line - David Roy
var y = oli.select(['B2','B3','B4','B5','B6','B7'],['B1', 'B2', 'B3', 'B4', 'B5', 'B7']) // select OLI bands 2-7 and rename them to match L7 band names
.resample('bicubic') // ...resample the L8 bands using bicubic
.subtract(itcp.multiply(10000)).divide(slopes) // ...multiply the y-intercept bands by 10000 to match the scale of the L7 bands then apply the line equation - subtract the intercept and divide by the slope
.set('system:time_start', oli.get('system:time_start')); // ...set the output system:time_start metadata to the input image time_start otherwise it is null
return y.toShort(); // return the image as short to match the type of the other data
};
//------ RETRIEVE A SENSOR SR COLLECTION FUNCTION -----
var getSRcollection = function(year, startDay, endDay, sensor, aoi) {
// get a landsat collection for given year, day range, and sensor
var srCollection = ee.ImageCollection('LANDSAT/'+ sensor + '/C01/T1_SR') // get surface reflectance images
.filterBounds(aoi) // ...filter them by intersection with AOI
.filterDate(year+'-'+startDay, year+'-'+endDay); // ...filter them by year and day range
// apply the harmonization function to LC08 (if LC08), subset bands, unmask, and resample
srCollection = srCollection.map(function(img) {
var dat = ee.Image(
ee.Algorithms.If(
sensor == 'LC08', // condition - if image is OLI
harmonizationRoy(img.unmask()), // true - then apply the L8 TO L7 alignment function after unmasking pixels that were previosuly masked (why/when are pixels masked)
img.select(['B1', 'B2', 'B3', 'B4', 'B5', 'B7']) // false - else select out the reflectance bands from the non-OLI image
.unmask() // ...unmask any previously masked pixels
.resample('bicubic') // ...resample by bicubic
.set('system:time_start', img.get('system:time_start')) // ...set the output system:time_start metadata to the input image time_start otherwise it is null
)
);
// make a cloud, cloud shadow, and snow mask from fmask band
var qa = img.select('pixel_qa'); // select out the fmask band
var mask = qa.bitwiseAnd(8).eq(0).and( // include shadow
qa.bitwiseAnd(16).eq(0)).and( // include snow
qa.bitwiseAnd(32).eq(0)); // include clouds
// apply the mask to the image and return it
return dat.mask(mask); //apply the mask - 0's in mask will be excluded from computation and set to opacity=0 in display
});
return srCollection; // return the prepared collection
};
//------ FUNCTION TO COMBINE LT05, LE07, & LC08 COLLECTIONS -----
var getCombinedSRcollection = function(year, startDay, endDay, aoi) {
var lt5 = getSRcollection(year, startDay, endDay, 'LT05', aoi); // get TM collection for a given year, date range, and area
var le7 = getSRcollection(year, startDay, endDay, 'LE07', aoi); // get ETM+ collection for a given year, date range, and area
var lc8 = getSRcollection(year, startDay, endDay, 'LC08', aoi); // get OLI collection for a given year, date range, and area
var mergedCollection = ee.ImageCollection(lt5.merge(le7).merge(lc8)); // merge the individual sensor collections into one imageCollection object
return mergedCollection; // return the Imagecollection
};
//------ FUNCTION TO REDUCE COLLECTION TO SINGLE IMAGE PER YEAR BY MEDOID -----
/*
LT expects only a single image per year in a time series, there are lost of ways to
do best available pixel compositing - we have found that a mediod composite requires little logic
is robust, and fast
Medoids are representative objects of a data set or a cluster with a data set whose average
dissimilarity to all the objects in the cluster is minimal. Medoids are similar in concept to
means or centroids, but medoids are always members of the data set.
*/
// make a medoid composite with equal weight among indices
var medoidMosaic = function(inCollection, dummyCollection) {
// fill in missing years with the dummy collection
var imageCount = inCollection.toList(1).length(); // get the number of images
var finalCollection = ee.ImageCollection(ee.Algorithms.If(imageCount.gt(0), inCollection, dummyCollection)); // if the number of images in this year is 0, then use the dummy collection, otherwise use the SR collection
// calculate median across images in collection per band
var median = finalCollection.median(); // calculate the median of the annual image collection - returns a single 6 band image - the collection median per band
// calculate the different between the median and the observation per image per band
var difFromMedian = finalCollection.map(function(img) {
var diff = ee.Image(img).subtract(median).pow(ee.Image.constant(2)); // get the difference between each image/band and the corresponding band median and take to power of 2 to make negatives positive and make greater differences weight more
return diff.reduce('sum').addBands(img); // per image in collection, sum the powered difference across the bands - set this as the first band add the SR bands to it - now a 7 band image collection
});
// get the medoid by selecting the image pixel with the smallest difference between median and observation per band
return ee.ImageCollection(difFromMedian).reduce(ee.Reducer.min(7)).select([1,2,3,4,5,6], ['B1','B2','B3','B4','B5','B7']); // find the powered difference that is the least - what image object is the closest to the median of teh collection - and then subset the SR bands and name them - leave behind the powered difference band
};
//------ FUNCTION TO APPLY MEDOID COMPOSITING FUNCTION TO A COLLECTION -------------------------------------------
var buildMosaic = function(year, startDay, endDay, aoi, dummyCollection) { // create a temp variable to hold the upcoming annual mosiac
var collection = getCombinedSRcollection(year, startDay, endDay, aoi); // get the SR collection
var img = medoidMosaic(collection, dummyCollection) // apply the medoidMosaic function to reduce the collection to single image per year by medoid
.set('system:time_start', (new Date(year,8,1)).valueOf()); // add the year to each medoid image - the data is hard-coded Aug 1st
return ee.Image(img); // return as image object
};
//------ FUNCTION TO BUILD ANNUAL MOSAIC COLLECTION ------------------------------
var buildMosaicCollection = function(startYear, endYear, startDay, endDay, aoi, dummyCollection) {
var imgs = []; // create empty array to fill
for (var i = startYear; i <= endYear; i++) { // for each year from hard defined start to end build medoid composite and then add to empty img array
var tmp = buildMosaic(i, startDay, endDay, aoi, dummyCollection); // build the medoid mosaic for a given year
imgs = imgs.concat(tmp.set('system:time_start', (new Date(i,8,1)).valueOf())); // concatenate the annual image medoid to the collection (img) and set the date of the image - hard coded to the year that is being worked on for Aug 1st
}
return ee.ImageCollection(imgs); // return the array img array as an image collection
};
//########################################################################################################
//##### FUNCTIONS FOR EXTRACTING AND PLOTTING A PIXEL TIME SERIES #####
//########################################################################################################
// ----- FUNCTION TO GET LT DATA FOR A PIXEL -----
var getPoint = function(img, geom, z) {
return img.reduceRegion({
reducer: 'first',
geometry: geom,
scale: z
}).getInfo();
};
// ----- FUNCTION TO CHART THE SOURCE AND FITTED TIME SERIES FOR A POINT -----
var chartPoint = function(lt, pt, distDir) {
Map.centerObject(pt, 14);
Map.addLayer(pt, {color: "FF0000"});
var point = getPoint(lt, pt, 10);
var data = [['x', 'y-original', 'y-fitted']];
for (var i = 0; i <= (endYear-startYear); i++) {
data = data.concat([[point.LandTrendr[0][i], point.LandTrendr[1][i]*distDir, point.LandTrendr[2][i]*distDir]]);
}
print(ui.Chart(data, 'LineChart',
{
'hAxis':
{
'format':'####'
},
'vAxis':
{
'maxValue': 1000,
'minValue': -1000
}
},
{'columns': [0, 1, 2]}
)
);
};
//########################################################################################################
//##### BUILD COLLECTION AND RUN LANDTRENDR #####
//########################################################################################################
//----- BUILD LT COLLECTION -----
// build annual surface reflection collection
var aoi = ee.Geometry.Point(long, lat);
var annualSRcollection = buildMosaicCollection(startYear, endYear, startDay, endDay, aoi, dummyCollection); // put together the cloud-free medoid surface reflectance annual time series collection
// 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 = annualSRcollection.map(segIndex) // map the function over every image in the collection - returns a 1-band annual image collection of the spectral index
.map(function(img) {return img.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'))}); // ...set the output system:time_start metadata to the input image time_start otherwise it is null
//----- 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
//----- PLOT THE SOURCE AND FITTED TIME SERIES FOR THE GIVEN POINT -----
chartPoint(lt, aoi, distDir); // plot the x-y time series for the given point
标签：曲线拟合,img,image,time,collection,var,GEE,year,LandTrendr 来源： https://blog.csdn.net/qq_35591253/article/details/115333551