Earth Engine LandTrendr fitted RGB thumbnail time series

ee_landtrendr_fitted_rgb_thumbnails.py

"""
Copyright 2020 Justin Braaten

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

    https://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.
"""

import ee
import requests 
import os

################################################################################
# Inputs

aoi = ee.Geometry.Polygon(
  [[[-64.81494128992148, -15.5314882713746],
    [-64.81494128992148, -15.838227596066043],
    [-64.69958484460898, -15.838227596066043],
    [-64.69958484460898, -15.5314882713746]]], None, False);

outDir = 'C:/Users/user/proj/ee-lt-ts/mamore-river/'

# Setup vars to get dates.
startYear = 1984
endYear = 2019
startMonth = 6
startDay = 1
nDays = 150

################################################################################

# Make outdir if it does not exist.
if not os.path.exists(outDir):
  os.makedirs(outDir)

# Define a collection filter by date, bounds, and quality.
def colFilter(col, aoi):#, startDate, endDate):
  return(col.filterBounds(aoi))

# Landsat collection preprocessingEnabled
# Get Landsat surface reflectance collections for OLI, ETM+ and TM sensors.
LC08col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
LE07col = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
LT05col = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
LT04col = ee.ImageCollection('LANDSAT/LT04/C01/T1_SR')

# Define a collection filter by date, bounds, and quality.
def colFilter(col, aoi, startDate, endDate):
  return(col
    .filterBounds(aoi)
    .filterDate(startDate, endDate))
    #.filter('CLOUD_COVER < 5')
    #.filter('GEOMETRIC_RMSE_MODEL < 15')
    #.filter('IMAGE_QUALITY == 9 || IMAGE_QUALITY_OLI == 9'))

# Function to get and rename bands of interest from OLI.
def renameOli(img):
  return(img.select(
      ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
      ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']))

# Function to get and rename bands of interest from ETM+.
def renameEtm(img):
  return(img.select(
      ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
      ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']))

# Add NBR for LandTrendr segmentation.
def calcNbr(img):
  return(img.addBands(img.normalizedDifference(['NIR', 'SWIR2'])
    .multiply(-10000).rename('NBR')).int16())

# Define function to mask out clouds and cloud shadows in images.
# Use CFmask band included in USGS Landsat SR image product.
def fmask(img):
  cloudShadowBitMask = 1 << 3
  cloudsBitMask = 1 << 5
  qa = img.select('pixel_qa')
  mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) \
    .And(qa.bitwiseAnd(cloudsBitMask).eq(0))
  return(img.updateMask(mask))

# Define function to prepare OLI images.
def prepOli(img):
  orig = img
  img = renameOli(img)
  img = fmask(img)
  return (ee.Image(img.copyProperties(orig, orig.propertyNames()))
    .resample('bicubic'))

# Define function to prepare ETM+ images.
def prepEtm(img):
  orig = img
  img = renameEtm(img)
  img = fmask(img)
  return(ee.Image(img.copyProperties(orig, orig.propertyNames()))
    .resample('bicubic'))

# Get annual median collection.    
def getAnnualComp(y):
  startDate = ee.Date.fromYMD(
    ee.Number(y), ee.Number(startMonth), ee.Number(startDay))
  endDate = startDate.advance(ee.Number(nDays), 'day')
  
  # Filter collections and prepare them for merging.
  LC08coly = colFilter(LC08col, aoi, startDate, endDate).map(prepOli)
  LE07coly = colFilter(LE07col, aoi, startDate, endDate).map(prepEtm)
  LT05coly = colFilter(LT05col, aoi, startDate, endDate).map(prepEtm)
  LT04coly = colFilter(LT04col, aoi, startDate, endDate).map(prepEtm)
  
  # Merge the collections.
  col = LC08coly.merge(LE07coly).merge(LT05coly).merge(LT04coly)
  
  yearImg = col.median()
  nBands = yearImg.bandNames().size()
  yearImg = ee.Image(ee.Algorithms.If(
    nBands,
    yearImg,
    dummyImg))
  return(calcNbr(yearImg)
    .set({'year': y, 'system:time_start': startDate.millis(), 'nBands': nBands}))


################################################################################

# Make a dummy image for missing years.
bandNames = ee.List(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa'])
fillerValues = ee.List.repeat(0, bandNames.size())
dummyImg = ee.Image.constant(fillerValues).rename(bandNames) \
  .selfMask().int16()

################################################################################
# Get a list of years
years = ee.List.sequence(startYear, endYear)

################################################################################
# Make list of annual image composites.
imgList = years.map(getAnnualComp)

# Convert image composite list to collection
imgCol = ee.ImageCollection.fromImages(imgList)

################################################################################
# Run LandTrendr.
lt = ee.Algorithms.TemporalSegmentation.LandTrendr(
  timeSeries=imgCol.select(['NBR', 'SWIR1', 'NIR', 'Green']),
  maxSegments=10,
  spikeThreshold=0.7,
  vertexCountOvershoot=3,
  preventOneYearRecovery=True,
  recoveryThreshold=0.5,
  pvalThreshold=0.05,
  bestModelProportion=0.75,
  minObservationsNeeded=6)

################################################################################
# Get fitted imagery. This starts export tasks.
def getYearStr(year):
  return(ee.String('yr_').cat(ee.Algorithms.String(year).slice(0,4)))

yearsStr = years.map(getYearStr)

r = lt.select(['SWIR1_fit']).arrayFlatten([yearsStr]).toShort()
g = lt.select(['NIR_fit']).arrayFlatten([yearsStr]).toShort()
b = lt.select(['Green_fit']).arrayFlatten([yearsStr]).toShort()

for i, c in zip([r, g, b], ['r', 'g', 'b']):
  descr = 'mamore-river-'+c
  name = 'users/user/'+descr
  print(name)
  task = ee.batch.Export.image.toAsset(
   image=i,
   region=aoi.getInfo()['coordinates'],
   assetId=name,
   description=descr,
   scale=30,
   crs='EPSG:3857',
   maxPixels=1e13)
  task.start()


################################################################################
################################################################################

# Once export tasks are finished run this section.

"""
r = ee.Image('users/user/mamore-river-r')
g = ee.Image('users/user/mamore-river-g')
b = ee.Image('users/user/mamore-river-b')

def getLtRgb(year):
  return(ee.Image.cat(
  r.select([year]),
  g.select([year]),
  b.select([year]))
  .visualize(min=100, max=3500, gamma=0.8))

yearsStrPy = yearsStr.getInfo()
for yr in yearsStrPy:
  print(yr)
  img_url = getLtRgb(yr) \
    .getThumbUrl({
      'dimensions': 1200,
      'region': aoi.getInfo()['coordinates']})
  img_data = requests.get(img_url).content
  with open(outDir+'/'+yr+'.png', 'wb') as handler:
    handler.write(img_data)
"""

Dear Justin Braaten,

Thank you for providing the code. I would like to ask to check the possibility of removing the col.median() and instead get every available image such and used only col. I tried this but I get error ''ImageCollection' object has no attribute 'bandNames''. What changes would you suggest me to bring in order to analyse time series of every available image in the collection not only the median of all images.
Best,
Najib

Hi Najib, this script is using the LandTrendr algorithm to build an annual fitted image collection. The LandTrendr algorithm is expecting only a single image per year, so we composite multiple images per year as a representative. Maybe you can tell me what your end goal is and I can point you at a better script for your purpose.

Thank you for the prompt reply Justin, I really appreciate it. I just don't want the composite image. I have changed the code to do the analysis on monthly composite image (not yearly). As there is possibility of two image (acquisitions) during a month I get a single image to rund NDWI. I don't want to make a composite image and want to use individual Landsat image for every ~16 days. I am observing the area of water bodies (in dams) and would like to monitor the availability of water for Agriculture. The best would be to have the time series of ~16 days.
Please let me know if you want to know more .

You may just need this subset of the code then, it will return all TM, ETM+, and OLI images in the calendar date range that intersect your area of interest. It will ensure all of the band names are the same between the collections and mask each image for cloud and cloud shadows using the CFmask mask band. If you want to include other processing like adding NDWI, you can map another function over the merged image collection to calculate and add the band to each image. (the code is untested)

startDay = 175  # day of year
endDay = 245  # day of year
aoi = <ee.Geometry>  # area of interest


# Landsat collection preprocessingEnabled
# Get Landsat surface reflectance collections for OLI, ETM+ and TM sensors.
LC08col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
LE07col = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
LT05col = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
LT04col = ee.ImageCollection('LANDSAT/LT04/C01/T1_SR')

# Define a collection filter by calendar date and bounds.
def colFilter(col, aoi, startDay, endDay):
  return(col
    .filterBounds(aoi)
    .filter(ee.Filter.calendarRange(startDay, endDay, 'day_of_year'))

# Function to get and rename bands of interest from OLI.
def renameOli(img):
  return(img.select(
      ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
      ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']))

# Function to get and rename bands of interest from ETM+.
def renameEtm(img):
  return(img.select(
      ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
      ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']))

# Define function to mask out clouds and cloud shadows in images.
# Use CFmask band included in USGS Landsat SR image product.
def fmask(img):
  cloudShadowBitMask = 1 << 3
  cloudsBitMask = 1 << 5
  qa = img.select('pixel_qa')
  mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) \
    .And(qa.bitwiseAnd(cloudsBitMask).eq(0))
  return(img.updateMask(mask))

# Define function to prepare OLI images.
def prepOli(img):
  orig = img
  img = renameOli(img)
  img = fmask(img)
  return (ee.Image(img.copyProperties(orig, orig.propertyNames()))
    .resample('bicubic'))

# Define function to prepare ETM+ images.
def prepEtm(img):
  orig = img
  img = renameEtm(img)
  img = fmask(img)
  return(ee.Image(img.copyProperties(orig, orig.propertyNames()))
    .resample('bicubic'))
	
# Filter collections and prepare them for merging.
LC08coly = colFilter(LC08col, aoi, startDay, endDay).map(prepOli)
LE07coly = colFilter(LE07col, aoi, startDay, endDay).map(prepEtm)
LT05coly = colFilter(LT05col, aoi, startDay, endDay).map(prepEtm)
LT04coly = colFilter(LT04col, aoi, startDay, endDay).map(prepEtm)

# Merge the collections.
col = LC08coly.merge(LE07coly).merge(LT05coly).merge(LT04coly)

Dear Justin,
I would like to know if you could let me know why masking is introducing the effect as of Landsat7 images (removing pixel like parallel lines). The above images shows the NDWI values for masked and unmasked. It is reducing the quality of results. Do you have any suggestion to improve the mask?

Second Point: The LC08col = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR") is not anymore updated and ee.ImageCollection("LANDSAT/LC08/C02/T1_L2") is used as the new source based on Google Earth Engine message. The band combinations are different, do you have any plans to update your code?

Best regards,
Najib