Hfrhyu

I want to execute Continuous Change Detection and Classification(CCDC, Zhu&Woodcock, 2014) to all pixels in a given region in GEE.

The first step is to build an image collection containing all landsat images according to the path and row number assigned with. Here is my code:

# Define the path and row
path = 120
row = 38

# Define the start and finish time
start = ee.Date.fromYMD(1983, 1, 1)
finish = ee.Date.fromYMD(2019, 1, 1)

# Select Landsat bands respectively as their different band configurations
l8_bandlist = ee.List(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'pixel_qa'])
l7_bandlist = ee.List(['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'B6', 'pixel_qa'])

# Unified the band names in the collection
rename_list = ee.List(['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'tbb', 'qa'])

# Function to rename the band names of images in the collection
def __renamelandsatbands(img):
return img.select(
 old_bandname,
 new_bandname
 )


# Function to calculate the evi
def __calculateEVI(img):
index = img.expression(
 '2.5 * ((NIR-RED) / (NIR +6 * RED -7.5* BLUE))', 
 'NIR': img.select('nir'),
 'RED': img.select('red'),
 'BLUE': img.select('blue')
 ).set('system:time_start', img.get('system:time_start'))
return img.addBands(index)

# Function to calculate the NBRT
def __calculateNBRT(img):
index = img.expression(
 '(NIR - 0.1 * SWIR * Temp) / (NIR + 0.1 * SWIR * Temp)',
 'NIR': img.select('nir'),
 'SWIR': img.select('swir1'),
 'Temp': img.select('tbb')
 ).set('system:time_start', img.get('system:time_start'))
return img.addBands(index)

#Landsat Collection
l8_sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filter(
ee.Filter.eq('WRS_PATH', path)).filter(ee.Filter.eq('WRS_ROW', row)).filterDate(
 start, finish).select(l8_bandlist, rename_list).sort('system:time_start')

l7_sr = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').filter(
ee.Filter.eq('WRS_PATH', path)).filter(ee.Filter.eq('WRS_ROW', row)).filterDate(
 start, finish).select(l7_bandlist, rename_list).sort('system:time_start')

l5_sr = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').filter(
ee.Filter.eq('WRS_PATH', path)).filter(ee.Filter.eq('WRS_ROW', row)).filterDate(
 start, finish).select(l7_bandlist, rename_list).sort('system:time_start')

# Add EVI and NBRT into Landsat SR dataset
#OLI
l8_srevi = l8_sr.map(__calculateEVI) 
l8_dataset = l8_srevi.map(__calculateNBRT)

#ETM+
l7_srevi = l7_sr.map(__calculateEVI)
l7_dataset = l7_srevi.map(__calculateNBRT)

# TM
l5_srevi = l5_sr.map(__calculateEVI)
l5_dataset = l5_srevi.map(__calculateNBRT)

# Rename images again
old_bandname = rename_list.add('constant').add('nir_1')
new_bandname = ee.List(['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'tbb', 'qa', 'evi', 'nbrt'])
l8_dataset = l8_dataset.map(__renamelandsatbands)
l7_dataset = l7_dataset.map(__renamelandsatbands)
l5_dataset = l5_dataset.map(__renamelandsatbands)

# Stack landsat series image collection
lcdataset = l8_dataset.merge(l7_dataset).merge(l5_dataset).sort('system:time_start')

According to these code above, i made a Landsat series image collection which contains all TM/ETM+/OLI images in path 120 and row 38. Each of the image contains ten bands, including seven surface reflectance bands('blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'tbb'), two derived feature bands('EVI', ‘NBRT') and one quality control band('qa').

After the data get prepared, What i want to do is to get the every pixel clear observation time series of all seven SR bands and the two derived feature bands which filtered by the qa band. Then, for every pixel, the first 24 observations are used to initialize a regression time series model, the break points detection start from the 25th clear observation based on a specific threshold determined by 3 folds of the RMSE of the regression model . The 25th observation will be identified as a break pointexceed if it exceeds the threshold, and the next 24 clear observations will be used to initialize a new model again. If the 25th observation is in the threshold, it will be joined into the first 24 observations to update the initial model, and the 26th observation will be assessed until all clear observations of the time series been checked.

AS i introduced above, the algorithm is based on pixel scale and use circulative iteration of nine time series trajectorys to find the break point in each pixel. In other words, using CCDC to an image collection facing two iterations, the time series break points iterative evaluation and loop this algorithm pixel by pixel to the whole image coverage. Is there a best way to do this in GEE?