Run LandTrendr and Visualize Outputs

  • Creates a time series of Landsat composites, runs LandTrendr, and visualizes outputs

  • You can optionally export these outputs

Copyright 2025 Ian Housman

Licensed under the Apache License, Version 2.0 (the “License”); you may not use this file except in compliance with the License. You may obtain a copy of the License at

http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an “AS IS” BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

github github

#Example of how to get Landsat data using the getImagesLib, create median composites, run LandTrendr and then filter 
#LandTrendr output into usable data depicting where, when, and the magnitude of loss and gain
####################################################################################################
import os,sys
sys.path.append(os.getcwd())

#Module imports
try:
    import geeViz.getImagesLib as gil
except:
    !python -m pip install geeViz
    import geeViz.getImagesLib as gil


import geeViz.changeDetectionLib as cdl
ee = gil.ee

#Set up to mapper objects to use
#Can use the default one first
Map = gil.Map


print('done')

Initializing GEE
Successfully initialized
done

  • We will first take a look at the Global Forest Change output to see what LandTrendr can potentiall provide that it cannot

  • Take note how well it does with fast, abrupt changes. The Global Forest Change product does however struggle with long-term gradual changes common in the Western US

####################################################################################################
#Start function calls
####################################################################################################
#First, let's look at the Hansen Global Forest Change output
#This is a great product to get an idea of where loss has occurred 
Map.port = 1231
#First clear the map in case it has been populated with layers/commands earlier
Map.clearMap()

hansenStartYear = 2001
hansenEndYear = 2023

#Bring in Hansen data and add it to the map
hansen = ee.Image("UMD/hansen/global_forest_change_2023_v1_11").select(['lossyear']).add(2000).int16()
hansen = hansen.updateMask(hansen.neq(2000).And(hansen.gte(hansenStartYear)).And(hansen.lte(hansenEndYear)))
Map.addLayer(hansen,{'min':hansenStartYear,'max':hansenEndYear,'palette':cdl.lossYearPalette},'Hansen Loss Year',True)

#Bring in map
Map.turnOnInspector()
Map.view()

Adding layer: Hansen Loss Year
Starting webmap
Using default refresh token for geeView
Local web server at: http://localhost:1231/geeView/ already serving.
cwd a:\GEE\gee_py_modules_package\geeViz\examples
geeView URL: http://localhost:1231/geeView/?projectID=nlcd-tcc&accessToken=ya29.a0AeXRPp5QMy-DIsVjbb2NStjqAPKpsV_eQMJk6wpi7Ngvo-co2338OIP-IedaCoRPPeCIMZJfdAcoAYlkam_ydkkryoC1l-PuuogLeJ5_JfxfGsMNMxdMLVMkShQ-tZWPrC-gCkLq6BN1G-lSkcnpcmKtnuF_wgOU4idP4HwkBTUaCgYKAWwSARESFQHGX2Mi0b17bTJPlPc0L7VEu0Kn2g0178&accessTokenCreationTime=1741205464270

Create Annual Median Composites

  • LandTrendr is intended to be run with an annual time series

  • We will first create median composites for a range of years

Set up parameters

  • Below, we will provide a number of parameters for creating composites

  • We will just use the default Landsat image processing parameters. To learn how to improve the quality of composites, see the getLandsatWrapperNotebook.ipynb

# Define user parameters:

# Specify study area: Study area
# Can be a featureCollection, feature, or geometry
studyArea = gil.testAreas['CA']

# Update the startJulian and endJulian variables to indicate your seasonal 
# constraints. This supports wrapping for tropics and southern hemisphere.
# If using wrapping and the majority of the days occur in the second year, the system:time_start will default 
# to June 1 of that year.Otherwise, all system:time_starts will default to June 1 of the given year
# startJulian: Starting Julian date 
# endJulian: Ending Julian date
startJulian = 152
endJulian = 273

# Specify start and end years for all analyses
# More than a 3 year span should be provided for time series methods to work 
# well. If using Fmask as the cloud/cloud shadow masking method, or providing
# pre-computed stats for cloudScore and TDOM, this does not 
# matter
startYear = 1990  
endYear = 2024

Create Composites

  • Now we have some basic parameters defined, we will create simple annual median composites

  • This method will only use the Fmask cloud and cloud shadow masking methods. To learn how to improve the quality of composites, see the getLandsatWrapperNotebook.ipynb

  • We can visualize and inspect the composites in the Viewer.

  • Note: layer drawing and querying can be slow since all computation is being performed in GEE on-th-fly.

####################################################################################################
#Clear the map in case it has been populated with layers/commands earlier
Map.clearMap()

#Get images and then create median composites
allImages = gil.getProcessedLandsatScenes(studyArea,startYear,endYear,startJulian,endJulian)
dummyImage = allImages.first()
composites = ee.ImageCollection(ee.List.sequence(startYear,endYear)
                                .map(lambda yr: 
                                     gil.fillEmptyCollections(
                                         allImages.filter(ee.Filter.calendarRange(yr,yr,'year')),
                                         dummyImage)
                                     .median()
                                     .set('system:time_start',ee.Date.fromYMD(yr,6,1).millis())
                                    ))
Map.addTimeLapse(composites,gil.vizParamsFalse,'Composite Time Series')


#Bring in map
Map.centerObject(studyArea)
Map.turnOnInspector()
Map.view()

Get Processed Landsat: 
Start date: Jun 01 1990 , End date: Sep 29 2024
Applying scale factors for C2 L4 data
Applying scale factors for C2 L5 data
Applying scale factors for C2 L8 data
Only including SLC On Landsat 7
Applying scale factors for C2 L7 data
Applying scale factors for C2 L9 data
Applying Fmask Cloud Mask
Applying Fmask Shadow Mask
Adding layer: Composite Time Series
Starting webmap
Using default refresh token for geeView
Local web server at: http://localhost:8001/geeView/ already serving.
cwd a:\GEE\gee_py_modules_package\geeViz\examples
geeView URL: http://localhost:8001/geeView/?projectID=nlcd-tcc&accessToken=ya29.a0AeXRPp7XJBwco-XNr51PHM9PTl6ZFfJ-_dFI603HWNNZo3ivF9aWUhuw7lUAbS-Rro_YRnOxNt5y5a7I_cNpNWZwqFVJ4toS0OlC5brL56BXEa6K-7kD00njQtieIRBVDr8iOc-vTYl3UBXDHVgnTz-3X7K1CDebCB7zMGTUYX0aCgYKAaMSARESFQHGX2Mipura1U-rfBoTfqfO-gPzWw0178&accessTokenCreationTime=1741205796978

Run LandTrendr

Set up LandTrendr Parameters

  • LandTrendr has many parameters

  • We will set a few below. There are some additional parameters used by geeViz to help facilitate running and post-processing LandTrendr outputs to make them easier to use.

# Choose band or index
# NBR, NDMI, and NDVI tend to work best
# Other good options are wetness and tcAngleBG
indexName = 'NBR'

# How many significant loss and/or gain segments to include
# Do not make less than 1
# If you only want the first loss and/or gain, choose 1
# Generally any past 2 are noise
howManyToPull = 1

# Parameters to identify suitable LANDTRENDR segments

# Thresholds to identify loss in vegetation
# Any segment that has a change magnitude or slope less than both of these thresholds is omitted
lossMagThresh = -0.15
lossSlopeThresh = -0.05


# Thresholds to identify gain in vegetation
# Any segment that has a change magnitude or slope greater than both of these thresholds is omitted
gainMagThresh = 0.1
gainSlopeThresh = 0.05

# Number of years of duration to separate between slow and fast loss (>= this number will be called slow loss)
slowLossDurationThresh = 3

# Which segment to show change from
# Choose from: 'newest','oldest','largest','smallest','steepest','mostGradual','shortest','longest'
chooseWhichLoss = 'largest'
chooseWhichGain = 'largest'

# LandTrendr Params
# run_params ={
#   'timeSeries': (ImageCollection) Yearly time-series from which to extract breakpoints. The first band is used to find breakpoints, and all subsequent bands are fitted using those breakpoints.
#   'maxSegments':            6,\ (Integer) Maximum number of segments to be fitted on the time series.
#   'spikeThreshold':         0.9,\ (Float, default: 0.9) Threshold for damping the spikes (1.0 means no damping).
#   'vertexCountOvershoot':   3,\(Integer, default: 3) The initial model can overshoot the maxSegments + 1 vertices by this amount. Later, it will be pruned down to maxSegments + 1.
#   'preventOneYearRecovery': False,\(Boolean, default: False): Prevent segments that represent one year recoveries.
#   'recoveryThreshold':      0.25,\(Float, default: 0.25) If a segment has a recovery rate faster than 1/recoveryThreshold (in years), then the segment is disallowed.
#   'pvalThreshold':          0.05,\(Float, default: 0.1) If the p-value of the fitted model exceeds this threshold, then the current model is discarded and another one is fitted using the Levenberg-Marquardt optimizer.
#   'bestModelProportion':    0.75,\(Float, default: 0.75) Allows models with more vertices to be chosen if their p-value is no more than (2 - bestModelProportion) times the p-value of the best model.
#   'minObservationsNeeded':  6\(Integer, default: 6) Min observations needed to perform output fitting.
# };

# Define landtrendr params
run_params = { \
  'maxSegments':            6,\
  'spikeThreshold':         0.9,\
  'vertexCountOvershoot':   3,\
  'preventOneYearRecovery': False,\
  'recoveryThreshold':      0.25,\
  'pvalThreshold':          0.05,\
  'bestModelProportion':    0.75,\
  'minObservationsNeeded':  6\
}

# Whether to add change outputs to map
addToMap = True

# Export params
# Whether to export LANDTRENDR change detection (loss and gain) outputs
exportLTLossGain = False

# Whether to export LandTrendr vertex array raw output
exportLTVertexArray = False

# Set up Names for the export
outputName = 'LT_Test'

# Provide location composites will be exported to
# This should be an asset imageCollection
exportPathRoot = 'users/username/someCollection'

# CRS- must be provided.  
# Common crs codes: Web mercator is EPSG:4326, USGS Albers is EPSG:5070, 
# WGS84 UTM N hemisphere is EPSG:326+ zone number (zone 12 N would be EPSG:32612) and S hemisphere is EPSG:327+ zone number
crs = 'EPSG:5070'

# Specify transform if scale is None and snapping to known grid is needed
transform = [30,0,-2361915.0,0,-30,3177735.0]

# Specify scale if transform is None
scale = None
####################################################################################################
#End user parameters
####################################################################################################
print('Done')

Done

Run LandTrendr

  • This will run LandTrendr using the median annual composites we created above.

  • It will then post-process the output to provide a year of the largest vegetation cover loss, gain, along with the magnitude and duration of any loss or gain.

  • You can inspect any of these layers. To see the classic raw vs LandTrendr fitted output, make sure the “Raw and Fitted Time Series” is turned on.

  • Note: layer drawing and querying can be slow since all computation is being performed in GEE on-th-fly.

#Clear the map in case it has been populated with layers/commands earlier
Map.clearMap()
# Important that the range of data values of the composites match the run_params spikeThreshold and recoveryThreshold
# e.g. Reflectance bands that have a scale of 0-1 should have a spikeThreshold around 0.9
# and a recoveryThreshold around 0.25
# If the reflectance values were scaled by 10000, the spikeThreshold would be around 9000 
# and a recoveryThreshold around 2500
#Run LANDTRENDR
ltOutputs = cdl.simpleLANDTRENDR(composites,startYear,endYear,indexName, run_params,lossMagThresh,lossSlopeThresh,\
                                                gainMagThresh,gainSlopeThresh,slowLossDurationThresh,chooseWhichLoss,\
                                                chooseWhichGain,addToMap,howManyToPull)

#Bring in map
Map.turnOnInspector()
Map.centerObject(studyArea)
Map.addLayer(studyArea, {'strokeColor': '0000FF'}, "Study Area", False)
Map.view()

Converting LandTrendr from array output to Gain & Loss
Adding layer: Raw and Fitted Time Series
Adding layer: LandTrendr NBR Loss Year
Adding layer: LandTrendr NBR Loss Magnitude
Adding layer: LandTrendr NBR Loss Duration
Adding layer: LandTrendr NBR Gain Year
Adding layer: LandTrendr NBR Gain Magnitude
Adding layer: LandTrendr NBR Gain Duration
Adding layer: Study Area
Starting webmap
Using default refresh token for geeView
Local web server at: http://localhost:8001/geeView/ already serving.
cwd a:\GEE\gee_py_modules_package\geeViz\examples
geeView URL: http://localhost:8001/geeView/?projectID=nlcd-tcc&accessToken=ya29.a0AeXRPp7suyHEbMSG-0bO3LZtluvzk2xzW7J9maTpsDkhf-GYzSat6li0aWSCZGvthd86O42PNTM_3R-pfnQcqRgMscf99gyRI06OX-2dnQVFodYWz4BoJJqehLaSX7BzOuRAnpN85GQu98op07QhF5X5VkeIdjm-TmpOidQ32H8aCgYKAccSARESFQHGX2MiXSBld9eDW4-PCObpUmjrxA0178&accessTokenCreationTime=1741206123836

Exporting

  • Outputs can be exported to an asset for future use.

  • You can export the post-processed loss and gain change detection output for simple change detection needs and/or the raw LandTrendr output for use in more complex custom workflows such as use in machine learning.

# Export outputs if selected
if exportLTLossGain:
  lossGainStack = ltOutputs[1]
  #Export  stack
  exportName = f'{outputName}_LT_LossGain_Stack_{indexName}_{startYear}_{endYear}_{startJulian}_{endJulian}'
  exportPath = exportPathRoot + '/'+ exportName

  lossGainStack = lossGainStack.set({'startYear':startYear,
                                        'endYear':endYear,
                                        'startJulian':startJulian,
                                        'endJulian':endJulian,
                                        'band':indexName})
  lossGainStack =lossGainStack.set(run_params)
  
  #Set up proper resampling for each band
  #Be sure to change if the band names for the exported image change
  pyrObj = {'_yr_':'mode','_dur_':'mode','_mag_':'mean','_slope_':'mean'}
  possible = ['loss','gain']
  how_many_list = ee.List.sequence(1,howManyToPull).getInfo()
  outObj = {}
  for p in possible:
    for key in pyrObj.keys():
      for i in how_many_list:
        i = int(i)
        kt = indexName + '_LT_'+p + key+str(i)
        outObj[kt]= pyrObj[key]

  # print(outObj)
  # Export output
  gil.exportToAssetWrapper(lossGainStack,exportName,exportPath,outObj,studyArea,scale,crs,transform)


# Export raw LandTrendr array image
if exportLTVertexArray:
  rawLTForExport = ltOutputs[0]
  # Map.addLayer(rawLTForExport,{},'Raw LT For Export {}'.format(indexName),False)
  
  rawLTForExport = rawLTForExport.set({'startYear':startYear,
                                        'endYear':endYear,
                                        'startJulian':startJulian,
                                        'endJulian':endJulian,
                                        'band':indexName})
  rawLTForExport =rawLTForExport.set(run_params)
  exportName = '{}_LT_Raw_{}_{}_{}_{}_{}'.format(outputName,indexName,startYear,endYear,startJulian,endJulian)
  exportPath = exportPathRoot + '/'+ exportName
  gil.exportToAssetWrapper(rawLTForExport,exportName,exportPath,{'.default':'sample'},studyArea,scale,crs,transform)