Run LandTrendr and Visualize Outputs#
Gets a time series of Landsat composites, runs LandTrendr, and visualizes the output
You can optionally export these outputs
Copyright 2024 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.
#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
Map1 = gil.Map
#Set up another
Map2 = gil.mapper()
print('done')
Initializing GEE
Cached project id file path: C:\Users\ihousman\.config\earthengine\credentials.proj_id
Cached project id: lcms-292214
Successfully initialized
geeViz package folder: c:\Users\ihousman\AppData\Local\Programs\Python\Python311\Lib\site-packages\geeViz
done
# 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
# 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
slowLossDurationThresh = 3
# 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
####################################################################################################
#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
Map1.port = 1231
#First clear the map in case it has been populated with layers/commands earlier
Map1.clearMap()
#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(startYear)).And(hansen.lte(endYear)))
Map1.addLayer(hansen,{'min':startYear,'max':endYear,'palette':cdl.lossYearPalette},'Hansen Loss Year',True)
#Bring in map
Map1.turnOnInspector()
Map1.centerObject(studyArea)
Map1.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=lcms-292214&accessToken=ya29.a0AeDClZDAyPF5nRhKZA4tLqBJuKamOPKUJlYPXvx_4IPBza3iLvVBpwhpLrSWfOt_OIBrm_MZDnYKUCiBCdOrjuuREK9QG5HxYeZ3Dz-EwcqD7vl6Dm1W2HCtAVz8b_5Ldhzh3Mo69gzUZ1NvWa4hgsq6eoGh_O5iETqEAyo_3sIaCgYKAZ0SARESFQHGX2MiTPzEKVCYqICaTZoILxhr4Q0178
####################################################################################################
#Clear the map in case it has been populated with layers/commands earlier
Map2.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())
))
Map2.addTimeLapse(composites,gil.vizParamsFalse,'Composite Time Series')
#Bring in map
Map2.centerObject(studyArea)
Map2.turnOnInspector()
Map2.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=lcms-292214&accessToken=ya29.a0AeDClZBeTiKfQyHMwd8gXMqd7OJSPWvEdkLXj21ZHHomYbeqVxXPnAYhldyu5siFLHUaDATMB1dQuBcF4WEAd0BdG97m9iFuLEe5Kq0CvE5ikXsvkdVVRKhSHka14FEDl2o4-PDlG4a0-cfLM-TU96aOkTbiWH9cEhbTf62AoWsaCgYKAX8SARESFQHGX2MiQosDow0Kr8IqVtQ6Q81KMg0178
#Clear the map in case it has been populated with layers/commands earlier
Map1.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
Map1.turnOnInspector()
Map1.centerObject(studyArea)
Map1.addLayer(studyArea, {'strokeColor': '0000FF'}, "Study Area", False)
Map1.view()
Converting LandTrendr from array output to Gain & Loss
Adding layer: Raw and Fitted Time Series
Adding layer: 1 NBR Loss Year
Adding layer: 1 NBR Loss Magnitude
Adding layer: 1 NBR Loss Duration
Adding layer: 1 NBR Gain Year
Adding layer: 1 NBR Gain Magnitude
Adding layer: 1 NBR Gain Duration
Adding layer: Study Area
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=lcms-292214&accessToken=ya29.a0AeDClZBbdJ7B_a_E-eZ8D8xiEyKbPNP9fin7igxNZDffmCiWDn0YQ-jm7jTDj7YdiKU-yLB-VqCwmr-efOQ-jCveOUE3siDXMM8mKgqKs-VesI-FbrNaAF5odepkdqaJCtC4rRCVlVAhcGeeGIa928g9FxgHVpqUQmK-R33eZQgaCgYKAfwSARESFQHGX2MipsP5wJR8PGEeQWtKmu85RA0178
# 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)