Create High Quality Landsat Composites#

Demonstrates various parameters and their impact for making good composites over cloud areas

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 and view outputs using the Python visualization tools
#Acquires Landsat data and then adds them to the viewer
####################################################################################################
import os,sys
sys.path.append(os.getcwd())

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


ee = getImagesLib.ee
Map = getImagesLib.Map

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

Setting up numerous parameters available for processing Landsat data#

#This example will use all default parameters to demonstrate how to use the basic composite functionality

#First clear the map in case it has been populated with layers/commands earlier
Map.clearMap()
#Define user parameters:

# Specify study area: Study area
# Can be a featureCollection, feature, or geometry
studyArea = ee.FeatureCollection('projects/lcms-292214/assets/R8/PR_USVI/Ancillary/prusvi_boundary_buff2mile').geometry().bounds()#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 = 153
endJulian = 152

# 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 providing pre-computed stats for cloudScore and TDOM, this does not 
# matter
startYear = 2009
endYear = 2009

#Call on master wrapper function to get Landat scenes and composites
lsAndTs = getImagesLib.getLandsatWrapper(studyArea,startYear,endYear,startJulian,endJulian)


#Separate into scenes and composites for subsequent analysis
processedScenes = lsAndTs['processedScenes']
processedComposites = lsAndTs['processedComposites']



Map.clearMap()
# Map.addLayer(processedComposites.select(['NDVI','NBR']),{'addToLegend':'false'},'Time Series (NBR and NDVI)',False)
for year in range(startYear,endYear + 1 ):
     t = processedComposites.filter(ee.Filter.calendarRange(year,year,'year')).mosaic()
     Map.addLayer(t.float(),getImagesLib.vizParamsFalse,'Default Params {} {}-{}'.format(year,startJulian,endJulian),'True')

Map.centerObject(ee.Geometry.Polygon(
        [[[-65.8337045819611, 18.329538797654042],
          [-65.8337045819611, 18.235653085671174],
          [-65.70461522649235, 18.235653085671174],
          [-65.70461522649235, 18.329538797654042]]], None, False))
Map.turnOnInspector()
Map.view()

Get Processed Landsat: 
Start date: Jun 02 2009 , End date: Jun 01 2010
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

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[2], line 28
endYear = 2009
#Call on master wrapper function to get Landat scenes and composites
---> 28 lsAndTs = getImagesLib.getLandsatWrapper(studyArea,startYear,endYear,startJulian,endJulian)
#Separate into scenes and composites for subsequent analysis
processedScenes = lsAndTs['processedScenes']

File c:\Users\ihousman\AppData\Local\Programs\Python\Python311\Lib\site-packages\geeViz\getImagesLib.py:5432, in getLandsatWrapper(studyArea, startYear, endYear, startJulian, endJulian, timebuffer, weights, compositingMethod, toaOrSR, includeSLCOffL7, defringeL5, applyCloudScore, applyFmaskCloudMask, applyTDOM, applyFmaskCloudShadowMask, applyFmaskSnowMask, cloudScoreThresh, performCloudScoreOffset, cloudScorePctl, zScoreThresh, shadowSumThresh, contractPixels, dilatePixels, correctIllumination, correctScale, exportComposites, outputName, exportPathRoot, crs, transform, scale, resampleMethod, preComputedCloudScoreOffset, preComputedTDOMIRMean, preComputedTDOMIRStdDev, compositingReducer, harmonizeOLI, landsatCollectionVersion, overwrite, verbose)
   ls = ls.map(lambda img: addZenithAzimuth(img, toaOrSR))
# Create composite time series
-> 5432 ts = compositeTimeSeries(
   ls=ls,
   startYear=startYear,
   endYear=endYear,
   startJulian=startJulian,
   endJulian=endJulian,
   timebuffer=timebuffer,
   weights=weights,
   compositingMethod=compositingMethod,
   compositingReducer=compositingReducer,
)
# Correct illumination
if correctIllumination:

File c:\Users\ihousman\AppData\Local\Programs\Python\Python311\Lib\site-packages\geeViz\getImagesLib.py:4315, in compositeTimeSeries(ls, startYear, endYear, startJulian, endJulian, timebuffer, weights, compositingMethod, compositingReducer)
   return composite.set(
       {
           "system:time_start": ee.Date.fromYMD(year + yearWithMajority, 6, 1).millis(),
   (...)
       }
   )
# Iterate across each year
-> 4315 ts = [yearCompositeGetter(yr) for yr in ee.List.sequence(startYear + timebuffer, endYear - timebuffer).getInfo()]
ts = ee.ImageCollection(ts).set(args)
return ts

File c:\Users\ihousman\AppData\Local\Programs\Python\Python311\Lib\site-packages\geeViz\getImagesLib.py:4315, in <listcomp>(.0)
   return composite.set(
       {
           "system:time_start": ee.Date.fromYMD(year + yearWithMajority, 6, 1).millis(),
   (...)
       }
   )
# Iterate across each year
-> 4315 ts = [yearCompositeGetter(yr) for yr in ee.List.sequence(startYear + timebuffer, endYear - timebuffer).getInfo()]
ts = ee.ImageCollection(ts).set(args)
return ts

File c:\Users\ihousman\AppData\Local\Programs\Python\Python311\Lib\site-packages\geeViz\getImagesLib.py:4290, in compositeTimeSeries.<locals>.yearCompositeGetter(year)
images = yearsTT.map(yrGetter)
lsT = ee.ImageCollection(ee.FeatureCollection(images).flatten())
-> 4290 count = lsT.select([0]).count().rename(["compositeObsCount"])
# Compute median or medoid or apply reducer
if compositingReducer != None:

TypeError: Collection.count() missing 1 required positional argument: 'property'

#It is clear the default parameters do not work very well in this area
#There are missing data and cloud artifacts
#You can access the parameters that were used through the properties of the returned collection
print(processedComposites.toDictionary().getInfo())
print(processedScenes.toDictionary().getInfo())

{'compositingMethod': 'medoid', 'compositingReducer': 'None', 'endJulian': 152, 'endYear': 2009, 'startJulian': 153, 'startYear': 2009, 'timebuffer': 0, 'weights': '[1]'}
{'addPixelQA': 'False', 'applyCloudScore': 'False', 'applyFmaskCloudMask': 'True', 'applyFmaskCloudShadowMask': 'True', 'applyFmaskSnowMask': 'False', 'applyTDOM': 'False', 'cloudScorePctl': 10, 'cloudScoreThresh': 10, 'defringeL5': 'False', 'endJulian': 152, 'endYear': 2009, 'harmonizeOLI': 'False', 'includeSLCOffL7': 'False', 'landsatCollectionVersion': 'C2', 'origin': 'Landsat', 'performCloudScoreOffset': 'True', 'preComputedCloudScoreOffset': 'None', 'preComputedTDOMIRMean': 'None', 'preComputedTDOMIRStdDev': 'None', 'resampleMethod': 'near', 'shadowSumBands': "['nir', 'swir1']", 'startJulian': 153, 'startYear': 2009, 'toaOrSR': 'SR', 'verbose': 'False', 'wrapOffset': 365, 'zScoreThresh': -1}

Improving composite outputs by including Landsat 7#

#Since there are not that many images available in this area for these years, let's try adding Landsat 7
includeSLCOffL7 = True
#Call on master wrapper function to get Landat scenes and composites
lsAndTs = getImagesLib.getLandsatWrapper(studyArea,startYear,endYear,startJulian,endJulian,includeSLCOffL7=includeSLCOffL7)


#Separate into scenes and composites for subsequent analysis
processedScenes = lsAndTs['processedScenes']
processedComposites = lsAndTs['processedComposites']

#Turn off layers from previous iteration
Map.turnOffAllLayers()

# Map.addLayer(processedComposites.select(['NDVI','NBR']),{'addToLegend':'false'},'Time Series (NBR and NDVI)',False)
for year in range(startYear,endYear + 1 ):
     t = processedComposites.filter(ee.Filter.calendarRange(year,year,'year')).mosaic()
     Map.addLayer(t.float(),getImagesLib.vizParamsFalse,'L7 added {} {}-{}'.format(year,startJulian,endJulian),'True')


Map.view()

#You'll notice this helps fill in the holes, but introduces many cloud-related artifacts

Get Processed Landsat: 
Start date: Jun 02 2009 , End date: Jun 01 2010
Applying scale factors for C2 L4 data
Applying scale factors for C2 L5 data
Applying scale factors for C2 L8 data
Including All 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: L7 added 2009 153-152
Starting webmap
Using default refresh token for geeView: C:\Users\ihousman/.config/earthengine/credentials
Local web server at: http://localhost:8001/geeView/ already serving.
cwd a:\GEE\gee_py_modules_package\geeViz\examples

Improving cloud masking using the `cloudScore` method#

#Let's try to improve the cloud masking. Fmask is used by default, but misses some clouds
#We'll try the adding in the cloudScore method
applyCloudScore = True

#Call on master wrapper function to get Landat scenes and composites
lsAndTs = getImagesLib.getLandsatWrapper(studyArea,startYear,endYear,startJulian,endJulian,includeSLCOffL7=includeSLCOffL7,applyCloudScore=applyCloudScore)


#Separate into scenes and composites for subsequent analysis
processedScenes = lsAndTs['processedScenes']
processedComposites = lsAndTs['processedComposites']

#Turn off layers from previous iteration
Map.turnOffAllLayers()

# Map.addLayer(processedComposites.select(['NDVI','NBR']),{'addToLegend':'false'},'Time Series (NBR and NDVI)',False)
for year in range(startYear,endYear + 1 ):
     t = processedComposites.filter(ee.Filter.calendarRange(year,year,'year')).mosaic()
     Map.addLayer(t.float(),getImagesLib.vizParamsFalse,'L7 and CloudScore added {} {}-{}'.format(year,startJulian,endJulian),'True')


Map.view()

#You'll notice this cleans up the cloud masking a lot

Get Processed Landsat: 
Applying scale factors for C2 L4 data
Applying scale factors for C2 L5 data
Applying scale factors for C2 L8 data
Including All Landsat 7
Applying scale factors for C2 L7 data
Applying scale factors for C2 L9 data
Applying Cloud Score
Computing cloudScore offset
Applying Fmask Cloud Mask
Applying Fmask Shadow Mask
Adding layer: L7 and CloudScore added 2009 153-152
Starting webmap
Using default refresh token for geeView: C:\Users\ihousman/.config/earthengine/credentials
Local web server at: http://localhost:8001/geeView/ already serving.
cwd a:\GEE\gee_py_modules_package\geeViz\examples

Improving cloud shadow masking using TDOM#

#You'll still notice there are some dark areas likely due to cloud shadow masking omission
#Fmasks's cloud shadow mask misses a lot typically. A temporal outlier method called the 
#Temporal Dark Outlier Mask (TDOM) works well with masking cloud shadows

#We'll try the cloudScore method
applyTDOM = True

#Call on master wrapper function to get Landat scenes and composites
#In order to identify dark outliers, we will extend the dates by 6 years to get a larger sample
lsAndTs = getImagesLib.getLandsatWrapper(studyArea,startYear-3,endYear+3,startJulian,endJulian,includeSLCOffL7=includeSLCOffL7,applyCloudScore=applyCloudScore,applyTDOM=applyTDOM)


#Separate into scenes and composites for subsequent analysis
processedScenes = lsAndTs['processedScenes']
processedComposites = lsAndTs['processedComposites']

#Turn off layers from previous iteration
Map.turnOffAllLayers()

# Map.addLayer(processedComposites.select(['NDVI','NBR']),{'addToLegend':'false'},'Time Series (NBR and NDVI)',False)
for year in range(startYear,endYear + 1 ):
     t = processedComposites.filter(ee.Filter.calendarRange(year,year,'year')).mosaic()
     Map.addLayer(t.float(),getImagesLib.vizParamsFalse,'CloudScore and TDOM added {} {}-{}'.format(year,startJulian,endJulian),'True')


Map.view()

#You'll notice this cleans up the cloud masking a lot

Get Processed Landsat: 
Applying scale factors for C2 L4 data
Applying scale factors for C2 L5 data
Applying scale factors for C2 L8 data
Including All Landsat 7
Applying scale factors for C2 L7 data
Applying scale factors for C2 L9 data
Applying Cloud Score
Computing cloudScore offset
Applying Fmask Cloud Mask
Applying TDOM Shadow Mask
Computing irMean for TDOM
Computing irStdDev for TDOM
Applying Fmask Shadow Mask
Adding layer: CloudScore and TDOM added 2009 153-152
Starting webmap
Using default refresh token for geeView: C:\Users\ihousman/.config/earthengine/credentials
Local web server at: http://localhost:8001/geeView/ already serving.
cwd a:\GEE\gee_py_modules_package\geeViz\examples

There are many different parameters that can be changed in order to improve composites in different study areas This is just one example. Other parameters include changing date ranges, and reducers