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,
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