This is an applied example of a workflow that combines the strengths of four technologies to interactively explore environmental time series data and to prototype more extensive analyses.
In [2]:
!pip install earthengine-api
import ee
ee.Authenticate()
ee.Initialize()
In [15]:
import numpy as np
import pandas as pd
import altair as alt
!pip install geopandas
import geopandas
!pip install eeconvert
import eeconvert
In [0]:
reReArgs = {
'reducer': ee.Reducer.mean(),
'geometry': ee.Geometry.Point([0, 0]),
'scale': 1000,
'crs': 'EPSG:5070',
'bestEffort': True,
'maxPixels': 1e14,
'tileScale': 4}
This function reduces the pixels intersecting a region to a statistic. It operates on a single ee.Image and returns an ee.Feature (after the reduction, the pixel values and geometry matter not, so just return an aspatial feature). The plan is to map this over an ee.ImageCollection.
reduceRegion method to the image using arguments provided by the global reReArgs dictionary.ee.Feature without geometry that contains the region reduction statistic result (an ee.dictionary), all image properties, and the derived date properties.
In [0]:
def regionReduce(img):
eeDate = img.date()
year = eeDate.get('year')
month = eeDate.getRelative('month', 'year')
doy = eeDate.getRelative('day', 'year')
date = eeDate.format('YYYY-MM-dd')
stat = img.reduceRegion(
reducer = reReArgs['reducer'],
geometry = reReArgs['geometry'],
scale = reReArgs['scale'],
crs = reReArgs['crs'],
bestEffort = reReArgs['bestEffort'],
maxPixels = reReArgs['maxPixels'],
tileScale = reReArgs['tileScale'])
return(ee.Feature(None, stat) \
.copyProperties(img, img.propertyNames())
.set({
'DOY': doy,
'Month': month,
'Year': year,
'Date': date}))
The result of the above region reduction function is an ee.FeatureCollection, but a list is easier to work with. The following takes a collection and a list of properties to aggregate into a list of lists.
regionReduce function to all images in the passed collection.null for the requested properties (ensure equal length property lists in the following step)getInfo()
In [0]:
def getReReList(col, props):
dict = col.map(regionReduce) \
.filter(ee.Filter.notNull(props)) \
.reduceColumns(
reducer = ee.Reducer.toList().repeat(len(props)),
selectors = props)
return(ee.List(dict.get('list')).getInfo())
In [0]:
def gdfToFc(gdf):
"""converts a geodataframe to a featurecollection
Args:
gdf (geoPandas.GeoDataFrame) : the input geodataframe
Returns:
fc (ee.FeatureCollection) : feature collection (server side)
"""
gdfCopy = gdf.copy()
gdfCopy["eeFeature"] = gdfCopy.apply(eeconvert.shapelyToEEFeature,1)
featureList = gdfCopy["eeFeature"].tolist()
fc = ee.FeatureCollection(featureList)
return fc
In [0]:
pdsi = ee.ImageCollection('IDAHO_EPSCOR/PDSI')
geojson_url = 'https://gist.githubusercontent.com/jirikadlec2/c503a718d16dfb2f65dbf4e5a6627d32/raw/a817c5629c10df73d83ebf4e30e46fc3f0f5a204/bhutan_test_sites.geojson'
aoi_geodataframe = geopandas.read_file(geojson_url)
aoi_geodataframe.head()
bhutan_sites = gdfToFc(aoi_geodataframe)
#sn = ee.FeatureCollection(geojson_url)
sn = bhutan_sites.filter(ee.Filter.eq('description', 'site3'))
#sn = ee.FeatureCollection('EPA/Ecoregions/2013/L3') \
# .filter(ee.Filter.eq('na_l3name', 'Sierra Nevada'))
In [0]:
# Reset global arguments.
reReArgs['scale'] = 5000
reReArgs['reducer'] = ee.Reducer.mean()
reReArgs['geometry'] = sn
reReArgs['crs'] = 'EPSG:3310'
# Get a list containing a series of lists, one for each property given.
snPdsiList = getReReList(pdsi, ['Year', 'Month', 'DOY', 'Date', 'pdsi'])
In [22]:
# Print the dimensions of the returned list.
print('n variables:', len(snPdsiList))
print('n observations:', len(snPdsiList[0]))
In [23]:
snPdsiDf = pd.DataFrame(snPdsiList)
print(snPdsiDf.shape)
In [24]:
snPdsiDf = snPdsiDf.transpose()
print(snPdsiDf.shape)
Preview the first 5 rows of the dataframe.
In [25]:
snPdsiDf.head(5)
Out[25]:
It's missing column names - add some and preview again.
In [26]:
snPdsiDf.columns = ['Year', 'Month', 'DOY', 'Date', 'PDSI']
snPdsiDf.head(5)
Out[26]:
Now check the data type of each column.
In [27]:
snPdsiDf.dtypes
Out[27]:
object is a string - set the data types and preview again.
In [28]:
snPdsiDf = snPdsiDf.astype({'Year': int, 'Month': int, 'DOY': int, 'Date': str, 'PDSI': float})
snPdsiDf.dtypes
Out[28]:
Get the month correct.
In [29]:
snPdsiDf['Month'] = snPdsiDf['Month'] + 1
snPdsiDf.head(5)
Out[29]:
In [0]:
cols = {'Year': int, 'Month': int, 'DOY': int, 'Date': str, 'PDSI': float}
def eeList2Df(list, cols):
df = pd.DataFrame(list).transpose()
df.columns = [k for k in cols.keys()]
return(df.astype(cols))
In [31]:
alt.Chart(snPdsiDf).mark_rect().encode(
x = 'Year:O',
y = 'Month:O',
color = alt.Color('PDSI:Q', scale=alt.Scale(scheme="redblue", domain=(-5, 5))),
tooltip = [
alt.Tooltip('Year:O', title='Year'),
alt.Tooltip('mean(Month):O', title='Month'),
alt.Tooltip('PDSI:Q', title='PDSI')
]).properties(width=700, height=300)
Out[31]:
In [32]:
alt.Chart(snPdsiDf).mark_bar(size=1).encode(
x = 'Date:T',
y = 'PDSI:Q',
color = alt.Color('PDSI:Q', scale=alt.Scale(scheme="redblue", domain=(-5, 5))),
tooltip = [
alt.Tooltip('Date:T', title='Date'),
alt.Tooltip('PDSI:Q', title='PDSI')
]
).properties(width=700, height=300)
Out[32]:
In [0]:
def calculate_ndvi(image):
ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI')
return image.addBands(ndvi)
In [0]:
ndvi = ee.ImageCollection('MODIS/006/MOD13A2').filterDate('2017-01-01', '2020-03-25').select('NDVI')
reReArgs['scale'] = 1000
snNdviList = getReReList(ndvi, ['Year', 'DOY', 'NDVI'])
Convert the list to a tidy dataframe and preview it.
In [47]:
cols = {'Year': int, 'DOY': int, 'NDVI': float}
snNdviDf = eeList2Df(snNdviList, cols)
snNdviDf.head(5)
Out[47]:
Remove the NDVI scaling and preview again.
In [48]:
snNdviDf['NDVI'] = snNdviDf['NDVI']/10000
snNdviDf.head(5)
Out[48]:
In [49]:
alt.Chart(snNdviDf).mark_line().encode(
alt.X('DOY:O'),
alt.Y('NDVI:Q', scale=alt.Scale(domain=(0.1, 0.88))),
alt.Color('Year:O', scale=alt.Scale(scheme="magma")),
tooltip=[
alt.Tooltip('Year:O', title='Year'),
alt.Tooltip('DOY:O', title='DOY'),
alt.Tooltip('NDVI:Q', title='NDVI')
]
).interactive().properties(width=700, height=300)
Out[49]:
Reduce the plot to median and interquartile range.
In [50]:
line = alt.Chart(snNdviDf).mark_line().encode(
x='DOY:O',
y='median(NDVI):Q',
).interactive()
band = alt.Chart(snNdviDf).mark_errorband(extent='iqr').encode(
x='DOY:O',
y=alt.Y('NDVI:Q', scale=alt.Scale(domain=(0.1, 0.7))),
).interactive().properties(width=700, height=300)
band + line
Out[50]:
A scatterplot is a good way to visualize the relationship between two variables. The x variable will be PDSI and the y variable will be NDVI. To create this plot both variables must exist in the same dataframe. Each row will be an observation in time and two columns for corresponding PDSI and NDVI observations. Right now PDSI and NDVI are in two different dataframes, the first step is to merge them.
In [0]:
In [0]:
ndvi = ee.ImageCollection('COPERNICUS/S2').filterDate('2017-01-01', '2020-03-25').select('NDVI')
reReArgs['scale'] = 1000
snNdviList = getReReList(ndvi, ['Year', 'DOY', 'NDVI'])
In [0]:
ndviDoy = [224, 272]
snNdviDfSub = snNdviDf[
(snNdviDf['DOY'] >= ndviDoy[0]) & (snNdviDf['DOY'] <= ndviDoy[1])
]
snNdviDfSub = snNdviDfSub.groupby('Year').agg('min')
In [0]:
pdsiDoy = [1, 272]
snPdsiDfSub = snPdsiDf[
(snPdsiDf['DOY'] >= pdsiDoy[0]) & (snPdsiDf['DOY'] <= pdsiDoy[1])
]
snPdsiDfSub = snPdsiDfSub.groupby('Year').agg('mean')
In [0]:
ndviPdsiCombo = pd.merge(snNdviDfSub, snPdsiDfSub, how='left', on='Year') \
.drop(columns=['DOY_x', 'DOY_y', 'Month']) \
.reset_index()
ndviPdsiCombo.head(5)
Out[0]:
Fit.
In [0]:
ndviPdsiCombo['Fit'] = np.poly1d(np.polyfit(ndviPdsiCombo['PDSI'], ndviPdsiCombo['NDVI'], 1))(ndviPdsiCombo['PDSI'])
ndviPdsiCombo.head(5)
Out[0]:
The dataframe is ready for plotting. Since this plot includes points and a line of best fit, two charts need to be created, one for each, and the combined.
In [0]:
# Point chart.
points = alt.Chart(ndviPdsiCombo).mark_circle(size=60).encode(
x=alt.X('PDSI:Q', scale=alt.Scale(domain=(-5, 5))),
y=alt.Y('NDVI:Q', scale=alt.Scale(domain=(0.4, 0.6))),
color=alt.Color('Year:O', scale=alt.Scale(scheme='magma')),
tooltip=['Year', 'PDSI', 'NDVI']
).interactive()
# Line chart.
fit = alt.Chart(ndviPdsiCombo).mark_line().encode(
x=alt.X('PDSI:Q', scale=alt.Scale(domain=(-5, 5))),
y=alt.Y('Fit:Q', scale=alt.Scale(domain=(0.4, 0.6))),
color=alt.value('#808080')
).properties(width=700, height=300)
# Combine charts.
fit + points
Out[0]:
As you can see there seems to be some degree of positive correlation between PDSI and NDVI. Some of the greatest outliers are 2016, 2017, 2018 - the three years following recovery from the long drought.
At a course regional scale there appears to be a relationship between drought and NDVI. This section will look more closely at effects of drought on vegetation at a patch level. Here, a Landsat time series collection is created for the period 1984-2019 to provide greater temporal context for change.
In [184]:
sn = bhutan_sites.filter(ee.Filter.eq('description', 'site3'))
# Import the Folium library.
import folium
# Set visualization parameters.
visParams = {'min':0, 'max':7000}
# Create a folium map object.
site_centroid = sn.geometry().centroid().getInfo()
site_lon = site_centroid['coordinates'][0]
site_lat = site_centroid['coordinates'][1]
fMap = folium.Map(
location=[site_lat, site_lon],
tiles='https://mt.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
attr='Map Data: Google',
zoom_start=17,
height=500,
width=700
)
# Add AOI to map.
folium.GeoJson(
sn.geometry().getInfo(),
name='geojson',
style_function=lambda x: {'fillColor': '#00000000', 'color': '#FFFFFF'},
).add_to(fMap)
# Add a lat lon popup.
folium.LatLngPopup().add_to(fMap)
# Display the map.
display(fMap)
Edit date window inputs based on NDVI curve plotted previously and set latitude and longitude variables from the map above.
In [0]:
# Start and end dates.
startDay = 1
endDay = 365
Prepare a Landsat surface reflectance collection 2017-present.
In [0]:
# Make lat and lon an `ee.Geometry.Point` - AOI.
point = sn.geometry().centroid()
polygon = sn.geometry()
#point = ee.Geometry.Point([Longitude, Latitude])
# Define function to get and rename bands of interest from OLI.
def renameOLI(img):
return(img.select(
ee.List(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa']),
ee.List(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa'])))
# Define function to get and rename bands of interest from ETM+.
def renameETM(img):
return(img.select(
ee.List(['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa']),
ee.List(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa'])))
# Define function to get and rename bands of interest from Sentinnel2
def renameS2(img):
return(img.select(
ee.List(['B2', 'B3', 'B4', 'B8A', 'B11', 'B12', 'SCL']),
ee.List(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa'])))
# Define function to mask out clouds and cloud shadows.
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 mask out clouds and cloud shadows.
def sclMask(img):
qa = img.select('pixel_qa')
mask = qa.eq(4)\
.Or(qa.eq(5))\
.Or(qa.eq(6))\
.Or(qa.eq(7))\
.Or(qa.eq(11))
return(img.updateMask(mask))
# Define function to add year.
def setYear(img):
year = ee.Image(img).date().get('year')
return(img.set('Year', year))
def setMonth(img):
month = ee.Image(img).date().get('month')
return img.set('Month', month)
# Define function to calculate NBR.
def calcNBR(img):
return(img.normalizedDifference(ee.List(['NIR', 'SWIR2'])).rename('NBR'))
# Define function to calculate NDVI.
def calcNDVI(img):
return(img.normalizedDifference(ee.List(['NIR', 'Red'])).rename('NDVI'))
def calcMSAVI2(img):
nir = img.select('NIR')
red = img.select('Red')
return img.expression(
'(2 * NIR + 1 - sqrt(pow((2 * NIR + 1), 2) - 8 * (NIR - RED)) ) / 2',
{
'NIR': img.select('NIR'),
'RED': img.select('Red')
}
).rename('MSAVI2')
#return ((2 * nir + 1 - math.sqrt( (2 * nir + 1)**2 - 8 * (nir - red) )) / 2).rename('MSAVI2')
# Define function to prepare OLI images.
def prepOLI(img):
orig = img
img = renameOLI(img)
img = fmask(img)
#img = calcNBR(img)
#img = calcNDVI(img)
img = calcMSAVI2(img)
img = img.copyProperties(orig, orig.propertyNames())
return(setYear(img))
# Define function to prepare ETM+ images.
def prepETM(img):
orig = img
img = renameETM(img)
img = fmask(img)
#img = calcNBR(img)
#img = calcNDVI(img)
img = calcMSAVI2(img)
img = img.copyProperties(orig, orig.propertyNames())
return(setYear(img))
# Define function to prepare S2 images.
def prepS2(img):
orig = img
img = renameS2(img)
img = sclMask(img)
#img = calcNBR(img)
#img = calcNDVI(img)
img = calcMSAVI2(img)
img = img.copyProperties(orig, orig.propertyNames())
return(setYear(img))
# get data
tmCol = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR")
etmCol = ee.ImageCollection("LANDSAT/LE07/C01/T1_SR")
oliCol = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
s2Col = ee.ImageCollection("COPERNICUS/S2_SR")
# Filter collections and prepare them for merging.
oliCol = oliCol.filterBounds(point).filter(ee.Filter.calendarRange(startDay, endDay, 'day_of_year')).map(prepOLI)
etmCol = etmCol.filterBounds(point).filter(ee.Filter.calendarRange(startDay, endDay, 'day_of_year')).map(prepETM)
tmCol = tmCol.filterBounds(point).filter(ee.Filter.calendarRange(startDay, endDay, 'day_of_year')).map(prepETM)
s2Col = s2Col.filterBounds(point).filter(ee.Filter.calendarRange(startDay, endDay, 'day_of_year')).map(prepS2)
# Merge the collections.
landsatCol = oliCol \
.merge(etmCol) \
.merge(tmCol) \
.merge(s2Col) \
.filterDate('2017-01-01', '2020-03-30')
# Get distinct year collection.
distinctYearMonthCol = landsatCol.distinct('YearMonth')
# Define a filter that identifies which images from the complete collection
# match the DOY from the distinct DOY collection.
joinFilter = ee.Filter.equals(leftField='YearMonth', rightField='YearMonth')
# Define a join.
join = ee.Join.saveAll('year_month_matches')
# Apply the join and convert the resulting FeatureCollection to an
# ImageCollection.
joinCol = ee.ImageCollection(join.apply(distinctYearMonthCol, landsatCol, joinFilter))
# Apply mean reduction among matching year collections.
def reduceByJoin(img):
yearCol = ee.ImageCollection.fromImages(ee.Image(img).get('year_month_matches'));
return(yearCol.reduce(ee.Reducer.mean()) \
.rename('MSAVI2') \
.set('system:time_start', ee.Image(img).date().update(day=1).millis()))
#landsatCol = joinCol.map(reduceByJoin)
Reduce the region.
In [0]:
# reset global args
reReArgs['reducer'] = ee.Reducer.mean()
reReArgs['scale'] = 30
reReArgs['geometry'] = polygon
# Get a list containing a series of lists, one for each property given.
pointNDVIList = getReReList(landsatCol, ['Date', 'MSAVI2'])
#pointNbrList = getReReList(landsatCol, ['Date', 'NBR'])
Transform and tidy the table.
In [182]:
cols = {'Date': str, 'MSAVI2': float}
pointNDVIDf = eeList2Df(pointNDVIList, cols)
pointNDVIDf.tail(5)
Out[182]:
Display intra-annual reduction a line plot.
In [183]:
alt.Chart(pointNDVIDf).mark_line().encode(
x='Date:T',
y='MSAVI2:Q',
tooltip=[
alt.Tooltip('Date:T', title='Date'),
alt.Tooltip('NDVI:Q', title='MSAVI2')
]
).interactive().properties(width=700, height=300)
Out[183]: