A16f Checkpoint

Import libraries¶

In [ ]:

Copied!

import ee
import geemap
import ee
import geemap

Create an interactive map¶

In [ ]:

Copied!

Map = geemap.Map(center=[40, -100], zoom=4)
Map = geemap.Map(center=[40, -100], zoom=4)

Add Earth Engine Python script¶

In [ ]:

Copied!





# Add Earth Engine dataset
image = ee.Image("USGS/SRTMGL1_003")

#  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#  Chapter:      A1.6 Health Applications
#  Checkpoint:   A16f
#  Author:       Dawn Nekorchuk
#  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Section 1: Data Import
woredas = ee.FeatureCollection(
    'projects/gee-book/assets/A1-6/amhara_woreda_20170207')
# Create region outer boundary to filter products on.
amhara = woredas.geometry().bounds()
gpm = ee.ImageCollection('NASA/GPM_L3/IMERG_V06')
LSTTerra8 = ee.ImageCollection('MODIS/061/MOD11A2') \
    .filterDate('2001-06-26', Date.now())
brdfReflect = ee.ImageCollection('MODIS/006/MCD43A4')
brdfQa = ee.ImageCollection('MODIS/006/MCD43A2')

# Visualize woredas with black borders and no fill.
# Create an empty image into which to paint the features, cast to byte.
empty = ee.Image().byte()
# Paint all the polygon edges with the same number and width.
outline = empty.paint({
    'featureCollection': woredas,
    'color': 1,
    'width': 1
})
# Add woreda boundaries to the map.
Map.setCenter(38, 11.5, 7)
Map.addLayer(outline, {
    'palette': '000000'
}, 'Woredas')

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 2: Handling of dates

# 2.1 Requested start and end dates.
reqStartDate = ee.Date('2021-10-01')
reqEndDate = ee.Date('2021-11-30')

# 2.2 LST Dates
# LST MODIS is every 8 days, and a user-requested date will likely not match.
# We want to get the latest previous image date,
#  i.e. the date the closest, but prior to, the requested date.
# We will filter later.
# Get date of first image.
LSTEarliestDate = LSTTerra8.first().date()
# Filter collection to dates from beginning to requested start date.
priorLstImgCol = LSTTerra8.filterDate(LSTEarliestDate,
    reqStartDate)
# Get the latest (max) date of this collection of earlier images.
LSTPrevMax = priorLstImgCol.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
LSTStartDate = ee.Date(LSTPrevMax.get('max'))
print('LSTStartDate', LSTStartDate)

# 2.3 Last available data dates
# Different variables have different data lags.
# Data may not be available in user range.
# To prevent errors from stopping script,
#  grab last available (if relevant) & filter at end.

# 2.3.1 Precipitation
# Calculate date of most recent measurement for gpm (of all time).
gpmAllMax = gpm.reduceColumns(ee.Reducer.max(), [
    'system:time_start'
])
gpmAllEndDateTime = ee.Date(gpmAllMax.get('max'))
# GPM every 30 minutes, so get just date part.
gpmAllEndDate = ee.Date.fromYMD({
    'year': gpmAllEndDateTime.get('year'),
    'month': gpmAllEndDateTime.get('month'),
    'day': gpmAllEndDateTime.get('day')
})

# If data ends before requested start, take last data date,
# otherwise use requested date.
precipStartDate = ee.Date(gpmAllEndDate.millis() \
    .min(reqStartDate.millis()))
print('precipStartDate', precipStartDate)

# 2.3.2 BRDF
# Calculate date of most recent measurement for brdf (of all time).
brdfAllMax = brdfReflect.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
brdfAllEndDate = ee.Date(brdfAllMax.get('max'))
# If data ends before requested start, take last data date,
# otherwise use the requested date.
brdfStartDate = ee.Date(brdfAllEndDate.millis() \
    .min(reqStartDate.millis()))
print('brdfStartDate', brdfStartDate)
print('brdfEndDate', brdfAllEndDate)

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 3: Precipitation

# Section 3.1: Precipitation filtering and dates

# Filter gpm by date, using modified start if necessary.
gpmFiltered = gpm \
    .filterDate(precipStartDate, reqEndDate.advance(1, 'day')) \
    .filterBounds(amhara) \
    .select('precipitationCal')

# Calculate date of most recent measurement for gpm
# (in the modified requested window).
gpmMax = gpmFiltered.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
gpmEndDate = ee.Date(gpmMax.get('max'))
precipEndDate = gpmEndDate
print('precipEndDate ', precipEndDate)

# Create a list of dates for the precipitation time series.
precipDays = precipEndDate.difference(precipStartDate, 'day')
precipDatesPrep = ee.List.sequence(0, precipDays, 1)

def makePrecipDates(n):
    return precipStartDate.advance(n, 'day')

precipDates = precipDatesPrep.map(makePrecipDates)

# Section 3.2: Calculate daily precipitation

# Function to calculate daily precipitation:
def calcDailyPrecip(curdate):
    curdate = ee.Date(curdate)
    curyear = curdate.get('year')
    curdoy = curdate.getRelative('day', 'year').add(1)
    totprec = gpmFiltered \
        .filterDate(curdate, curdate.advance(1, 'day')) \
        .select('precipitationCal') \
        .sum() \
        .multiply(0.5) \
        .rename('totprec')

    return totprec \
        .set('doy', curdoy) \
        .set('year', curyear) \
        .set('system:time_start', curdate)

# Map function over list of dates.
dailyPrecipExtended =
    ee.ImageCollection.fromImages(precipDates.map(calcDailyPrecip))

# Filter back to the original user requested start date.
dailyPrecip = dailyPrecipExtended \
    .filterDate(reqStartDate, precipEndDate.advance(1, 'day'))

# Section 3.3: Summarize daily precipitation by woreda

# Filter precip data for zonal summaries.
precipSummary = dailyPrecip \
    .filterDate(reqStartDate, reqEndDate.advance(1, 'day'))

# Function to calculate zonal statistics for precipitation by woreda.
def sumZonalPrecip(image):
    # To get the doy and year,
    # convert the metadata to grids and then summarize.
    image2 = image.addBands([
        image.metadata('doy').int(),
        image.metadata('year').int()
    ])
    # Reduce by regions to get zonal means for each county.
    output = image2.select(['year', 'doy', 'totprec']) \
        .reduceRegions({
            'collection': woredas,
            'reducer': ee.Reducer.mean(),
            'scale': 1000
        })
    return output

# Map the zonal statistics function over the filtered precip data.
precipWoreda = precipSummary.map(sumZonalPrecip)
# Flatten the results for export.
precipFlat = precipWoreda.flatten()

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 4: Land surface temperature

# Section 4.1: Calculate LST variables

# Filter Terra LST by altered LST start date.
# Rarely, but at the end of the year if the last image is late in the year
#  with only a few days in its period, it will sometimes not grab
#  the next image. Add extra padding to reqEndDate and
#  it will be trimmed at the end.
LSTFiltered = LSTTerra8 \
    .filterDate(LSTStartDate, reqEndDate.advance(8, 'day')) \
    .filterBounds(amhara) \
    .select('LST_Day_1km', 'QC_Day', 'LST_Night_1km', 'QC_Night')

# Filter Terra LST by QA information.
def filterLstQa(image):
    qaday = image.select(['QC_Day'])
    qanight = image.select(['QC_Night'])
    dayshift = qaday.rightShift(6)
    nightshift = qanight.rightShift(6)
    daymask = dayshift.lte(2)
    nightmask = nightshift.lte(2)
    outimage = ee.Image(image.select(['LST_Day_1km',
        'LST_Night_1km'
    ]))
    outmask = ee.Image([daymask, nightmask])
    return outimage.updateMask(outmask)

LSTFilteredQa = LSTFiltered.map(filterLstQa)

# Rescale temperature data and convert to degrees Celsius (C).
def rescaleLst(image):
    LST_day = image.select('LST_Day_1km') \
        .multiply(0.02) \
        .subtract(273.15) \
        .rename('LST_day')
    LST_night = image.select('LST_Night_1km') \
        .multiply(0.02) \
        .subtract(273.15) \
        .rename('LST_night')
    LST_mean = image.expression(
        '(day + night) / 2', {
            'day': LST_day.select('LST_day'),
            'night': LST_night.select('LST_night')
        }
    ).rename('LST_mean')
    return image.addBands(LST_day) \
        .addBands(LST_night) \
        .addBands(LST_mean)

LSTVars = LSTFilteredQa.map(rescaleLst)

# Section 4.2: Calculate daily LST

# Create list of dates for time series.
LSTRange = LSTVars.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
LSTEndDate = ee.Date(LSTRange.get('max')).advance(7, 'day')
LSTDays = LSTEndDate.difference(LSTStartDate, 'day')
LSTDatesPrep = ee.List.sequence(0, LSTDays, 1)

def makeLstDates(n):
    return LSTStartDate.advance(n, 'day')

LSTDates = LSTDatesPrep.map(makeLstDates)

# Function to calculate daily LST by assigning the 8-day composite summary
# to each day in the composite period:
def calcDailyLst(curdate):
    curyear = ee.Date(curdate).get('year')
    curdoy = ee.Date(curdate).getRelative('day', 'year').add(1)
    moddoy = curdoy.divide(8).ceil().subtract(1).multiply(8).add(
        1)
    basedate = ee.Date.fromYMD(curyear, 1, 1)
    moddate = basedate.advance(moddoy.subtract(1), 'day')
    LST_day = LSTVars \
        .select('LST_day') \
        .filterDate(moddate, moddate.advance(1, 'day')) \
        .first() \
        .rename('LST_day')
    LST_night = LSTVars \
        .select('LST_night') \
        .filterDate(moddate, moddate.advance(1, 'day')) \
        .first() \
        .rename('LST_night')
    LST_mean = LSTVars \
        .select('LST_mean') \
        .filterDate(moddate, moddate.advance(1, 'day')) \
        .first() \
        .rename('LST_mean')
    return LST_day \
        .addBands(LST_night) \
        .addBands(LST_mean) \
        .set('doy', curdoy) \
        .set('year', curyear) \
        .set('system:time_start', curdate)

# Map the function over the image collection
dailyLstExtended =
    ee.ImageCollection.fromImages(LSTDates.map(calcDailyLst))

# Filter back to original user requested start date
dailyLst = dailyLstExtended \
    .filterDate(reqStartDate, LSTEndDate.advance(1, 'day'))

# Section 4.3: Summarize daily LST by woreda

# Filter LST data for zonal summaries.
LSTSummary = dailyLst \
    .filterDate(reqStartDate, reqEndDate.advance(1, 'day'))
# Function to calculate zonal statistics for LST by woreda:
def sumZonalLst(image):
    # To get the doy and year, we convert the metadata to grids
    #  and then summarize.
    image2 = image.addBands([
        image.metadata('doy').int(),
        image.metadata('year').int()
    ])
    # Reduce by regions to get zonal means for each county.
    output = image2 \
        .select(['doy', 'year', 'LST_day', 'LST_night', 'LST_mean']) \
        .reduceRegions({
            'collection': woredas,
            'reducer': ee.Reducer.mean(),
            'scale': 1000
        })
    return output

# Map the zonal statistics function over the filtered LST data.
LSTWoreda = LSTSummary.map(sumZonalLst)
# Flatten the results for export.
LSTFlat = LSTWoreda.flatten()

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 5: Spectral index NDWI

# Section 5.1: Calculate NDWI

# Filter BRDF-Adjusted Reflectance by date.
brdfReflectVars = brdfReflect \
    .filterDate(brdfStartDate, reqEndDate.advance(1, 'day')) \
    .filterBounds(amhara) \
    .select([
            'Nadir_Reflectance_Band1', 'Nadir_Reflectance_Band2',
            'Nadir_Reflectance_Band3', 'Nadir_Reflectance_Band4',
            'Nadir_Reflectance_Band5', 'Nadir_Reflectance_Band6',
            'Nadir_Reflectance_Band7'
        ],
        ['red', 'nir', 'blue', 'green', 'swir1', 'swir2', 'swir3'])

# Filter BRDF QA by date.
brdfReflectQa = brdfQa \
    .filterDate(brdfStartDate, reqEndDate.advance(1, 'day')) \
    .filterBounds(amhara) \
    .select([
            'BRDF_Albedo_Band_Quality_Band1',
            'BRDF_Albedo_Band_Quality_Band2',
            'BRDF_Albedo_Band_Quality_Band3',
            'BRDF_Albedo_Band_Quality_Band4',
            'BRDF_Albedo_Band_Quality_Band5',
            'BRDF_Albedo_Band_Quality_Band6',
            'BRDF_Albedo_Band_Quality_Band7',
            'BRDF_Albedo_LandWaterType'
        ],
        ['qa1', 'qa2', 'qa3', 'qa4', 'qa5', 'qa6', 'qa7', 'water'])

# Join the 2 collections.
idJoin = ee.Filter.equals({
    'leftField': 'system:time_end',
    'rightField': 'system:time_end'
})
# Define the join.
innerJoin = ee.Join.inner('NBAR', 'QA')
# Apply the join.
brdfJoined = innerJoin.apply(brdfReflectVars, brdfReflectQa,
    idJoin)

# Add QA bands to the NBAR collection.
def addQaBands(image):
    nbar = ee.Image(image.get('NBAR'))
    qa = ee.Image(image.get('QA')).select(['qa2'])
    water = ee.Image(image.get('QA')).select(['water'])
    return nbar.addBands([qa, water])

brdfMerged = ee.ImageCollection(brdfJoined.map(addQaBands))

# Function to mask out pixels based on QA and water/land flags.
def filterBrdf(image):
    # Using QA info for the NIR band.
    qaband = image.select(['qa2'])
    wband = image.select(['water'])
    qamask = qaband.lte(2).And(wband.eq(1))
    nir_r = image.select('nir').multiply(0.0001).rename('nir_r')
    swir2_r = image.select('swir2').multiply(0.0001).rename(
        'swir2_r')
    return image.addBands(nir_r) \
        .addBands(swir2_r) \
        .updateMask(qamask)

brdfFilteredVars = brdfMerged.map(filterBrdf)

# Function to calculate spectral indices:
def calcBrdfIndices(image):
    curyear = ee.Date(image.get('system:time_start')).get('year')
    curdoy = ee.Date(image.get('system:time_start')) \
        .getRelative('day', 'year').add(1)
    ndwi6 = image.normalizedDifference(['nir_r', 'swir2_r']) \
        .rename('ndwi6')
    return image.addBands(ndwi6) \
        .set('doy', curdoy) \
        .set('year', curyear)

# Map function over image collection.
brdfFilteredVars = brdfFilteredVars.map(calcBrdfIndices)

# Section 5.2: Calculate daily NDWI

# Create list of dates for full time series.
brdfRange = brdfFilteredVars.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
brdfEndDate = ee.Date(brdfRange.get('max'))
brdfDays = brdfEndDate.difference(brdfStartDate, 'day')
brdfDatesPrep = ee.List.sequence(0, brdfDays, 1)

def makeBrdfDates(n):
    return brdfStartDate.advance(n, 'day')

brdfDates = brdfDatesPrep.map(makeBrdfDates)

# List of dates that exist in BRDF data.
brdfDatesExist = brdfFilteredVars \
    .aggregate_array('system:time_start')

# Get daily brdf values.
def calcDailyBrdfExists(curdate):
    curdate = ee.Date(curdate)
    curyear = curdate.get('year')
    curdoy = curdate.getRelative('day', 'year').add(1)
    brdfTemp = brdfFilteredVars \
        .filterDate(curdate, curdate.advance(1, 'day'))
    outImg = brdfTemp.first()
    return outImg

dailyBrdfExtExists =
    ee.ImageCollection.fromImages(brdfDatesExist.map(
        calcDailyBrdfExists))

# Create empty results, to fill in dates when BRDF data does not exist.
def calcDailyBrdfFiller(curdate):
    curdate = ee.Date(curdate)
    curyear = curdate.get('year')
    curdoy = curdate.getRelative('day', 'year').add(1)
    brdfTemp = brdfFilteredVars \
        .filterDate(curdate, curdate.advance(1, 'day'))
    brdfSize = brdfTemp.size()
    outImg = ee.Image.constant(0).selfMask() \
        .addBands(ee.Image.constant(0).selfMask()) \
        .addBands(ee.Image.constant(0).selfMask()) \
        .addBands(ee.Image.constant(0).selfMask()) \
        .addBands(ee.Image.constant(0).selfMask()) \
        .rename(['ndvi', 'evi', 'savi', 'ndwi5', 'ndwi6']) \
        .set('doy', curdoy) \
        .set('year', curyear) \
        .set('system:time_start', curdate) \
        .set('brdfSize', brdfSize)
    return outImg

# Create filler for all dates.
dailyBrdfExtendedFiller =
    ee.ImageCollection.fromImages(brdfDates.map(calcDailyBrdfFiller))
# But only used if and when size was 0.
dailyBrdfExtFillFilt = dailyBrdfExtendedFiller \
    .filter(ee.Filter.eq('brdfSize', 0))
# Merge the two collections.
dailyBrdfExtended = dailyBrdfExtExists \
    .merge(dailyBrdfExtFillFilt)

# Filter back to original user requested start date.
dailyBrdf = dailyBrdfExtended \
    .filterDate(reqStartDate, brdfEndDate.advance(1, 'day'))

# Section 5.3: Summarize daily spectral indices by woreda

# Filter spectral indices for zonal summaries.
brdfSummary = dailyBrdf \
    .filterDate(reqStartDate, reqEndDate.advance(1, 'day'))

# Function to calculate zonal statistics for spectral indices by woreda:
def sumZonalBrdf(image):
    # To get the doy and year, we convert the metadata to grids
    #  and then summarize.
    image2 = image.addBands([
        image.metadata('doy').int(),
        image.metadata('year').int()
    ])
    # Reduce by regions to get zonal means for each woreda.
    output = image2.select(['doy', 'year', 'ndwi6']) \
        .reduceRegions({
            'collection': woredas,
            'reducer': ee.Reducer.mean(),
            'scale': 1000
        })
    return output


# Map the zonal statistics function over the filtered spectral index data.
brdfWoreda = brdfSummary.map(sumZonalBrdf)
# Flatten the results for export.
brdfFlat = brdfWoreda.flatten()

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 6: Map display of calculated environmental variables
displayDate = ee.Date('2021-10-01')

precipDisp = dailyPrecip \
    .filterDate(displayDate, displayDate.advance(1, 'day'))
brdfDisp = dailyBrdf \
    .filterDate(displayDate, displayDate.advance(1, 'day'))
LSTDisp = dailyLst \
    .filterDate(displayDate, displayDate.advance(1, 'day'))

# Select the image (should be only one) from each collection.
precipImage = precipDisp.first().select('totprec')
LSTmImage = LSTDisp.first().select('LST_mean')
ndwi6Image = brdfDisp.first().select('ndwi6')

# Palettes for environmental variable maps:
palettePrecip = ['f7fbff', '08306b']
paletteLst = ['fff5f0', '67000d']
paletteSpectral = ['ffffe5', '004529']

# Add layers to the map.
# Show precipitation by default,
#  others hidden until users picks them from layers drop down.
Map.addLayer({
    'eeObject': precipImage,
    'visParams': {
        'min': 0,
        'max': 20,
        'palette': palettePrecip
    },
    'name': 'Precipitation',
    'shown': True,
    'opacity': 0.75
})
Map.addLayer({
    'eeObject': LSTmImage,
    'visParams': {
        'min': 0,
        'max': 40,
        'palette': paletteLst
    },
    'name': 'LST Mean',
    'shown': False,
    'opacity': 0.75
})
Map.addLayer({
    'eeObject': ndwi6Image,
    'visParams': {
        'min': 0,
        'max': 1,
        'palette': paletteSpectral
    },
    'name': 'NDWI6',
    'shown': False,
    'opacity': 0.75
})

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------
# Add Earth Engine dataset
image = ee.Image("USGS/SRTMGL1_003")

#  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#  Chapter:      A1.6 Health Applications
#  Checkpoint:   A16f
#  Author:       Dawn Nekorchuk
#  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Section 1: Data Import
woredas = ee.FeatureCollection(
    'projects/gee-book/assets/A1-6/amhara_woreda_20170207')
# Create region outer boundary to filter products on.
amhara = woredas.geometry().bounds()
gpm = ee.ImageCollection('NASA/GPM_L3/IMERG_V06')
LSTTerra8 = ee.ImageCollection('MODIS/061/MOD11A2') \
    .filterDate('2001-06-26', Date.now())
brdfReflect = ee.ImageCollection('MODIS/006/MCD43A4')
brdfQa = ee.ImageCollection('MODIS/006/MCD43A2')

# Visualize woredas with black borders and no fill.
# Create an empty image into which to paint the features, cast to byte.
empty = ee.Image().byte()
# Paint all the polygon edges with the same number and width.
outline = empty.paint({
    'featureCollection': woredas,
    'color': 1,
    'width': 1
})
# Add woreda boundaries to the map.
Map.setCenter(38, 11.5, 7)
Map.addLayer(outline, {
    'palette': '000000'
}, 'Woredas')

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 2: Handling of dates

# 2.1 Requested start and end dates.
reqStartDate = ee.Date('2021-10-01')
reqEndDate = ee.Date('2021-11-30')

# 2.2 LST Dates
# LST MODIS is every 8 days, and a user-requested date will likely not match.
# We want to get the latest previous image date,
#  i.e. the date the closest, but prior to, the requested date.
# We will filter later.
# Get date of first image.
LSTEarliestDate = LSTTerra8.first().date()
# Filter collection to dates from beginning to requested start date.
priorLstImgCol = LSTTerra8.filterDate(LSTEarliestDate,
    reqStartDate)
# Get the latest (max) date of this collection of earlier images.
LSTPrevMax = priorLstImgCol.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
LSTStartDate = ee.Date(LSTPrevMax.get('max'))
print('LSTStartDate', LSTStartDate)

# 2.3 Last available data dates
# Different variables have different data lags.
# Data may not be available in user range.
# To prevent errors from stopping script,
#  grab last available (if relevant) & filter at end.

# 2.3.1 Precipitation
# Calculate date of most recent measurement for gpm (of all time).
gpmAllMax = gpm.reduceColumns(ee.Reducer.max(), [
    'system:time_start'
])
gpmAllEndDateTime = ee.Date(gpmAllMax.get('max'))
# GPM every 30 minutes, so get just date part.
gpmAllEndDate = ee.Date.fromYMD({
    'year': gpmAllEndDateTime.get('year'),
    'month': gpmAllEndDateTime.get('month'),
    'day': gpmAllEndDateTime.get('day')
})

# If data ends before requested start, take last data date,
# otherwise use requested date.
precipStartDate = ee.Date(gpmAllEndDate.millis() \
    .min(reqStartDate.millis()))
print('precipStartDate', precipStartDate)

# 2.3.2 BRDF
# Calculate date of most recent measurement for brdf (of all time).
brdfAllMax = brdfReflect.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
brdfAllEndDate = ee.Date(brdfAllMax.get('max'))
# If data ends before requested start, take last data date,
# otherwise use the requested date.
brdfStartDate = ee.Date(brdfAllEndDate.millis() \
    .min(reqStartDate.millis()))
print('brdfStartDate', brdfStartDate)
print('brdfEndDate', brdfAllEndDate)

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 3: Precipitation

# Section 3.1: Precipitation filtering and dates

# Filter gpm by date, using modified start if necessary.
gpmFiltered = gpm \
    .filterDate(precipStartDate, reqEndDate.advance(1, 'day')) \
    .filterBounds(amhara) \
    .select('precipitationCal')

# Calculate date of most recent measurement for gpm
# (in the modified requested window).
gpmMax = gpmFiltered.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
gpmEndDate = ee.Date(gpmMax.get('max'))
precipEndDate = gpmEndDate
print('precipEndDate ', precipEndDate)

# Create a list of dates for the precipitation time series.
precipDays = precipEndDate.difference(precipStartDate, 'day')
precipDatesPrep = ee.List.sequence(0, precipDays, 1)

def makePrecipDates(n):
    return precipStartDate.advance(n, 'day')

precipDates = precipDatesPrep.map(makePrecipDates)

# Section 3.2: Calculate daily precipitation

# Function to calculate daily precipitation:
def calcDailyPrecip(curdate):
    curdate = ee.Date(curdate)
    curyear = curdate.get('year')
    curdoy = curdate.getRelative('day', 'year').add(1)
    totprec = gpmFiltered \
        .filterDate(curdate, curdate.advance(1, 'day')) \
        .select('precipitationCal') \
        .sum() \
        .multiply(0.5) \
        .rename('totprec')

    return totprec \
        .set('doy', curdoy) \
        .set('year', curyear) \
        .set('system:time_start', curdate)

# Map function over list of dates.
dailyPrecipExtended =
    ee.ImageCollection.fromImages(precipDates.map(calcDailyPrecip))

# Filter back to the original user requested start date.
dailyPrecip = dailyPrecipExtended \
    .filterDate(reqStartDate, precipEndDate.advance(1, 'day'))

# Section 3.3: Summarize daily precipitation by woreda

# Filter precip data for zonal summaries.
precipSummary = dailyPrecip \
    .filterDate(reqStartDate, reqEndDate.advance(1, 'day'))

# Function to calculate zonal statistics for precipitation by woreda.
def sumZonalPrecip(image):
    # To get the doy and year,
    # convert the metadata to grids and then summarize.
    image2 = image.addBands([
        image.metadata('doy').int(),
        image.metadata('year').int()
    ])
    # Reduce by regions to get zonal means for each county.
    output = image2.select(['year', 'doy', 'totprec']) \
        .reduceRegions({
            'collection': woredas,
            'reducer': ee.Reducer.mean(),
            'scale': 1000
        })
    return output

# Map the zonal statistics function over the filtered precip data.
precipWoreda = precipSummary.map(sumZonalPrecip)
# Flatten the results for export.
precipFlat = precipWoreda.flatten()

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 4: Land surface temperature

# Section 4.1: Calculate LST variables

# Filter Terra LST by altered LST start date.
# Rarely, but at the end of the year if the last image is late in the year
#  with only a few days in its period, it will sometimes not grab
#  the next image. Add extra padding to reqEndDate and
#  it will be trimmed at the end.
LSTFiltered = LSTTerra8 \
    .filterDate(LSTStartDate, reqEndDate.advance(8, 'day')) \
    .filterBounds(amhara) \
    .select('LST_Day_1km', 'QC_Day', 'LST_Night_1km', 'QC_Night')

# Filter Terra LST by QA information.
def filterLstQa(image):
    qaday = image.select(['QC_Day'])
    qanight = image.select(['QC_Night'])
    dayshift = qaday.rightShift(6)
    nightshift = qanight.rightShift(6)
    daymask = dayshift.lte(2)
    nightmask = nightshift.lte(2)
    outimage = ee.Image(image.select(['LST_Day_1km',
        'LST_Night_1km'
    ]))
    outmask = ee.Image([daymask, nightmask])
    return outimage.updateMask(outmask)

LSTFilteredQa = LSTFiltered.map(filterLstQa)

# Rescale temperature data and convert to degrees Celsius (C).
def rescaleLst(image):
    LST_day = image.select('LST_Day_1km') \
        .multiply(0.02) \
        .subtract(273.15) \
        .rename('LST_day')
    LST_night = image.select('LST_Night_1km') \
        .multiply(0.02) \
        .subtract(273.15) \
        .rename('LST_night')
    LST_mean = image.expression(
        '(day + night) / 2', {
            'day': LST_day.select('LST_day'),
            'night': LST_night.select('LST_night')
        }
    ).rename('LST_mean')
    return image.addBands(LST_day) \
        .addBands(LST_night) \
        .addBands(LST_mean)

LSTVars = LSTFilteredQa.map(rescaleLst)

# Section 4.2: Calculate daily LST

# Create list of dates for time series.
LSTRange = LSTVars.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
LSTEndDate = ee.Date(LSTRange.get('max')).advance(7, 'day')
LSTDays = LSTEndDate.difference(LSTStartDate, 'day')
LSTDatesPrep = ee.List.sequence(0, LSTDays, 1)

def makeLstDates(n):
    return LSTStartDate.advance(n, 'day')

LSTDates = LSTDatesPrep.map(makeLstDates)

# Function to calculate daily LST by assigning the 8-day composite summary
# to each day in the composite period:
def calcDailyLst(curdate):
    curyear = ee.Date(curdate).get('year')
    curdoy = ee.Date(curdate).getRelative('day', 'year').add(1)
    moddoy = curdoy.divide(8).ceil().subtract(1).multiply(8).add(
        1)
    basedate = ee.Date.fromYMD(curyear, 1, 1)
    moddate = basedate.advance(moddoy.subtract(1), 'day')
    LST_day = LSTVars \
        .select('LST_day') \
        .filterDate(moddate, moddate.advance(1, 'day')) \
        .first() \
        .rename('LST_day')
    LST_night = LSTVars \
        .select('LST_night') \
        .filterDate(moddate, moddate.advance(1, 'day')) \
        .first() \
        .rename('LST_night')
    LST_mean = LSTVars \
        .select('LST_mean') \
        .filterDate(moddate, moddate.advance(1, 'day')) \
        .first() \
        .rename('LST_mean')
    return LST_day \
        .addBands(LST_night) \
        .addBands(LST_mean) \
        .set('doy', curdoy) \
        .set('year', curyear) \
        .set('system:time_start', curdate)

# Map the function over the image collection
dailyLstExtended =
    ee.ImageCollection.fromImages(LSTDates.map(calcDailyLst))

# Filter back to original user requested start date
dailyLst = dailyLstExtended \
    .filterDate(reqStartDate, LSTEndDate.advance(1, 'day'))

# Section 4.3: Summarize daily LST by woreda

# Filter LST data for zonal summaries.
LSTSummary = dailyLst \
    .filterDate(reqStartDate, reqEndDate.advance(1, 'day'))
# Function to calculate zonal statistics for LST by woreda:
def sumZonalLst(image):
    # To get the doy and year, we convert the metadata to grids
    #  and then summarize.
    image2 = image.addBands([
        image.metadata('doy').int(),
        image.metadata('year').int()
    ])
    # Reduce by regions to get zonal means for each county.
    output = image2 \
        .select(['doy', 'year', 'LST_day', 'LST_night', 'LST_mean']) \
        .reduceRegions({
            'collection': woredas,
            'reducer': ee.Reducer.mean(),
            'scale': 1000
        })
    return output

# Map the zonal statistics function over the filtered LST data.
LSTWoreda = LSTSummary.map(sumZonalLst)
# Flatten the results for export.
LSTFlat = LSTWoreda.flatten()

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 5: Spectral index NDWI

# Section 5.1: Calculate NDWI

# Filter BRDF-Adjusted Reflectance by date.
brdfReflectVars = brdfReflect \
    .filterDate(brdfStartDate, reqEndDate.advance(1, 'day')) \
    .filterBounds(amhara) \
    .select([
            'Nadir_Reflectance_Band1', 'Nadir_Reflectance_Band2',
            'Nadir_Reflectance_Band3', 'Nadir_Reflectance_Band4',
            'Nadir_Reflectance_Band5', 'Nadir_Reflectance_Band6',
            'Nadir_Reflectance_Band7'
        ],
        ['red', 'nir', 'blue', 'green', 'swir1', 'swir2', 'swir3'])

# Filter BRDF QA by date.
brdfReflectQa = brdfQa \
    .filterDate(brdfStartDate, reqEndDate.advance(1, 'day')) \
    .filterBounds(amhara) \
    .select([
            'BRDF_Albedo_Band_Quality_Band1',
            'BRDF_Albedo_Band_Quality_Band2',
            'BRDF_Albedo_Band_Quality_Band3',
            'BRDF_Albedo_Band_Quality_Band4',
            'BRDF_Albedo_Band_Quality_Band5',
            'BRDF_Albedo_Band_Quality_Band6',
            'BRDF_Albedo_Band_Quality_Band7',
            'BRDF_Albedo_LandWaterType'
        ],
        ['qa1', 'qa2', 'qa3', 'qa4', 'qa5', 'qa6', 'qa7', 'water'])

# Join the 2 collections.
idJoin = ee.Filter.equals({
    'leftField': 'system:time_end',
    'rightField': 'system:time_end'
})
# Define the join.
innerJoin = ee.Join.inner('NBAR', 'QA')
# Apply the join.
brdfJoined = innerJoin.apply(brdfReflectVars, brdfReflectQa,
    idJoin)

# Add QA bands to the NBAR collection.
def addQaBands(image):
    nbar = ee.Image(image.get('NBAR'))
    qa = ee.Image(image.get('QA')).select(['qa2'])
    water = ee.Image(image.get('QA')).select(['water'])
    return nbar.addBands([qa, water])

brdfMerged = ee.ImageCollection(brdfJoined.map(addQaBands))

# Function to mask out pixels based on QA and water/land flags.
def filterBrdf(image):
    # Using QA info for the NIR band.
    qaband = image.select(['qa2'])
    wband = image.select(['water'])
    qamask = qaband.lte(2).And(wband.eq(1))
    nir_r = image.select('nir').multiply(0.0001).rename('nir_r')
    swir2_r = image.select('swir2').multiply(0.0001).rename(
        'swir2_r')
    return image.addBands(nir_r) \
        .addBands(swir2_r) \
        .updateMask(qamask)

brdfFilteredVars = brdfMerged.map(filterBrdf)

# Function to calculate spectral indices:
def calcBrdfIndices(image):
    curyear = ee.Date(image.get('system:time_start')).get('year')
    curdoy = ee.Date(image.get('system:time_start')) \
        .getRelative('day', 'year').add(1)
    ndwi6 = image.normalizedDifference(['nir_r', 'swir2_r']) \
        .rename('ndwi6')
    return image.addBands(ndwi6) \
        .set('doy', curdoy) \
        .set('year', curyear)

# Map function over image collection.
brdfFilteredVars = brdfFilteredVars.map(calcBrdfIndices)

# Section 5.2: Calculate daily NDWI

# Create list of dates for full time series.
brdfRange = brdfFilteredVars.reduceColumns({
    'reducer': ee.Reducer.max(),
    'selectors': ['system:time_start']
})
brdfEndDate = ee.Date(brdfRange.get('max'))
brdfDays = brdfEndDate.difference(brdfStartDate, 'day')
brdfDatesPrep = ee.List.sequence(0, brdfDays, 1)

def makeBrdfDates(n):
    return brdfStartDate.advance(n, 'day')

brdfDates = brdfDatesPrep.map(makeBrdfDates)

# List of dates that exist in BRDF data.
brdfDatesExist = brdfFilteredVars \
    .aggregate_array('system:time_start')

# Get daily brdf values.
def calcDailyBrdfExists(curdate):
    curdate = ee.Date(curdate)
    curyear = curdate.get('year')
    curdoy = curdate.getRelative('day', 'year').add(1)
    brdfTemp = brdfFilteredVars \
        .filterDate(curdate, curdate.advance(1, 'day'))
    outImg = brdfTemp.first()
    return outImg

dailyBrdfExtExists =
    ee.ImageCollection.fromImages(brdfDatesExist.map(
        calcDailyBrdfExists))

# Create empty results, to fill in dates when BRDF data does not exist.
def calcDailyBrdfFiller(curdate):
    curdate = ee.Date(curdate)
    curyear = curdate.get('year')
    curdoy = curdate.getRelative('day', 'year').add(1)
    brdfTemp = brdfFilteredVars \
        .filterDate(curdate, curdate.advance(1, 'day'))
    brdfSize = brdfTemp.size()
    outImg = ee.Image.constant(0).selfMask() \
        .addBands(ee.Image.constant(0).selfMask()) \
        .addBands(ee.Image.constant(0).selfMask()) \
        .addBands(ee.Image.constant(0).selfMask()) \
        .addBands(ee.Image.constant(0).selfMask()) \
        .rename(['ndvi', 'evi', 'savi', 'ndwi5', 'ndwi6']) \
        .set('doy', curdoy) \
        .set('year', curyear) \
        .set('system:time_start', curdate) \
        .set('brdfSize', brdfSize)
    return outImg

# Create filler for all dates.
dailyBrdfExtendedFiller =
    ee.ImageCollection.fromImages(brdfDates.map(calcDailyBrdfFiller))
# But only used if and when size was 0.
dailyBrdfExtFillFilt = dailyBrdfExtendedFiller \
    .filter(ee.Filter.eq('brdfSize', 0))
# Merge the two collections.
dailyBrdfExtended = dailyBrdfExtExists \
    .merge(dailyBrdfExtFillFilt)

# Filter back to original user requested start date.
dailyBrdf = dailyBrdfExtended \
    .filterDate(reqStartDate, brdfEndDate.advance(1, 'day'))

# Section 5.3: Summarize daily spectral indices by woreda

# Filter spectral indices for zonal summaries.
brdfSummary = dailyBrdf \
    .filterDate(reqStartDate, reqEndDate.advance(1, 'day'))

# Function to calculate zonal statistics for spectral indices by woreda:
def sumZonalBrdf(image):
    # To get the doy and year, we convert the metadata to grids
    #  and then summarize.
    image2 = image.addBands([
        image.metadata('doy').int(),
        image.metadata('year').int()
    ])
    # Reduce by regions to get zonal means for each woreda.
    output = image2.select(['doy', 'year', 'ndwi6']) \
        .reduceRegions({
            'collection': woredas,
            'reducer': ee.Reducer.mean(),
            'scale': 1000
        })
    return output


# Map the zonal statistics function over the filtered spectral index data.
brdfWoreda = brdfSummary.map(sumZonalBrdf)
# Flatten the results for export.
brdfFlat = brdfWoreda.flatten()

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

# Section 6: Map display of calculated environmental variables
displayDate = ee.Date('2021-10-01')

precipDisp = dailyPrecip \
    .filterDate(displayDate, displayDate.advance(1, 'day'))
brdfDisp = dailyBrdf \
    .filterDate(displayDate, displayDate.advance(1, 'day'))
LSTDisp = dailyLst \
    .filterDate(displayDate, displayDate.advance(1, 'day'))

# Select the image (should be only one) from each collection.
precipImage = precipDisp.first().select('totprec')
LSTmImage = LSTDisp.first().select('LST_mean')
ndwi6Image = brdfDisp.first().select('ndwi6')

# Palettes for environmental variable maps:
palettePrecip = ['f7fbff', '08306b']
paletteLst = ['fff5f0', '67000d']
paletteSpectral = ['ffffe5', '004529']

# Add layers to the map.
# Show precipitation by default,
#  others hidden until users picks them from layers drop down.
Map.addLayer({
    'eeObject': precipImage,
    'visParams': {
        'min': 0,
        'max': 20,
        'palette': palettePrecip
    },
    'name': 'Precipitation',
    'shown': True,
    'opacity': 0.75
})
Map.addLayer({
    'eeObject': LSTmImage,
    'visParams': {
        'min': 0,
        'max': 40,
        'palette': paletteLst
    },
    'name': 'LST Mean',
    'shown': False,
    'opacity': 0.75
})
Map.addLayer({
    'eeObject': ndwi6Image,
    'visParams': {
        'min': 0,
        'max': 1,
        'palette': paletteSpectral
    },
    'name': 'NDWI6',
    'shown': False,
    'opacity': 0.75
})

#  -----------------------------------------------------------------------
#  CHECKPOINT
#  -----------------------------------------------------------------------

Display the interactive map¶

In [ ]:

Copied!

Map
Map