!pip install earthengine-api
import ee

ee.Authenticate()
ee.Initialize()

"""
Explanation: <a href="https://colab.research.google.com/github/jirikadlec2/FetchSnow/blob/master/Copy_of_time_series_visualization.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>
Time Series Visualization
Exploring drought and vegetation response with Google Earth Engine and Altair in Google Colab
Justin Braaten<br>
Geo for Good Summit, 2019
bit.ly/ts-vis
<br><br><br><br><br><br>

Introduction
What this is not:

an introduction to Earth Engine.
an introduction to Pandas.
an introduction to Altair.
an introduction to Colab.

What this is:
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.
<br><br><br><br><br><br>

Turning data into insights
...and ultimately into a narrative
<br><br><br>
<img src='https://raw.githubusercontent.com/jdbcode/G4G19/master/img/data-to-insight.png' width='700px' alt='data-to-insight.png'>
<br><br><br><br><br><br>

Data reduction
<img src='https://raw.githubusercontent.com/jdbcode/G4G19/master/img/data-reduction.png' width='700px' alt='data-to-insight.png'>

Earth Engine - geospatial data access and processing
Pandas - dataframe structure
Altair - data visualization
Colab - integrative Python environment 

<br><br><br><br><br><br>

Drought and vegetation response
Sierra Nevada ecoregion
<br>
Source: Google | Earth Engine
<br>
Source: Flickr | Matt McGrath
<br><br><br><br><br><br>

Datasets

Drought severity (PDSI)
Vegetation response (MODIS NDVI)
Vegetation response (Landsat NBR)
Historical climate (PRISM)
Projected Climate (NEX-DCP30)

<br><br><br><br><br><br>

General workflow

Filter the dataset
Reduce region by a statistic
Convert data format from list to dataframe
Alter the dataframe
Plot the dataframe

<br><br><br><br><br><br>

Python setup
Install the Earth Engine API

Use pip to install the latest version of the Earth Engine API.
Import the Earth Engine library.
Authenticate access (registration verification and Google account access).
Initialize the API.
End of explanation
"""


import numpy as np
import pandas as pd
import altair as alt

!pip install geopandas
import geopandas

!pip install eeconvert
import eeconvert

"""
Explanation: Import other libraries
End of explanation
"""


reReArgs = {
  'reducer': ee.Reducer.mean(),
  'geometry': ee.Geometry.Point([0, 0]),
  'scale': 1000,
  'crs': 'EPSG:5070',
  'bestEffort': True,
  'maxPixels': 1e14,
  'tileScale': 4}

"""
Explanation: Define region reduction function
Reduction arguments
Define a global dictionary that provides arguments to the Earth Engine reduceRegion image method.
End of explanation
"""


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}))

"""
Explanation: Reduction function
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.

Derive a series of date properties that will be used in charting.
Apply the reduceRegion method to the image using arguments provided by the global reReArgs dictionary.
Return an ee.Feature without geometry that contains the region reduction statistic result (an ee.dictionary), all image properties, and the derived date properties.
End of explanation
"""


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())

"""
Explanation: Reduction to list
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.

Apply the above defined regionReduce function to all images in the passed collection.
Remove any features that are null for the requested properties (ensure equal length property lists in the following step)
Aggregate requested properties into a list of lists
Extract the list of lists from the resulting dictionary and then send the data to the Colab VM client with getInfo()
End of explanation
"""


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

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'))

"""
Explanation: Reduce drought severity by region
Load PDSI and AOI

Load gridded Palmer Drought Severity Index (PDSI) data as an ee.ImageCollection
Load the EPA Level-3 ecoregion boundaries as an ee.FeatureCollection and filter it to include only the Sierra Nevada region, which defines the area of interest (AOI).
End of explanation
"""


# 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'])

"""
Explanation: Reduce PDSI by AOI

Set the reduce region arguments.
Reduce PDSI for the AOI and return the results as a list. (return all properties in the provided list)
End of explanation
"""


# Print the dimensions of the returned list.
print('n variables:', len(snPdsiList))
print('n observations:', len(snPdsiList[0]))

"""
Explanation: Format the table
Check the format
Print the list dimensions. These data are in memory on the Colab VM.
End of explanation
"""


snPdsiDf = pd.DataFrame(snPdsiList)
print(snPdsiDf.shape)

"""
Explanation: Convert list of lists to dataframe

Convert the list to a Pandas dataframe.
Print the shape of the new dataframe. shape returns (n rows, n columns)
End of explanation
"""


snPdsiDf = snPdsiDf.transpose()
print(snPdsiDf.shape)

"""
Explanation: Transpose the dataframe

Transpose the dataframe.
Print the shape again.
End of explanation
"""


snPdsiDf.head(5)

"""
Explanation: Preview the first 5 rows of the dataframe.
End of explanation
"""


snPdsiDf.columns = ['Year', 'Month', 'DOY', 'Date', 'PDSI']
snPdsiDf.head(5)

"""
Explanation: It's missing column names - add some and preview again.
End of explanation
"""


snPdsiDf.dtypes

"""
Explanation: Now check the data type of each column.
End of explanation
"""


snPdsiDf = snPdsiDf.astype({'Year': int, 'Month': int, 'DOY': int, 'Date': str, 'PDSI': float})
snPdsiDf.dtypes

"""
Explanation: object is a string - set the data types and preview again.
End of explanation
"""


snPdsiDf['Month'] = snPdsiDf['Month'] + 1
snPdsiDf.head(5)

"""
Explanation: Get the month correct.
End of explanation
"""


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

"""
Explanation: Define a data structure translator
The above data structure transformation will be used again multiple times; make it a function.
End of explanation
"""


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)

"""
Explanation: Calendar heatmap
End of explanation
"""


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)

"""
Explanation: Bar chart
End of explanation
"""


def calculate_ndvi(image):
  ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI')
  return image.addBands(ndvi)

ndvi = ee.ImageCollection('MODIS/006/MOD13A2').filterDate('2017-01-01', '2020-03-25').select('NDVI')

reReArgs['scale'] = 1000
snNdviList = getReReList(ndvi, ['Year', 'DOY', 'NDVI'])

"""
Explanation: Reduce NDVI by region
A measure of photosynthetic capacity from MODIS.
Get data and reduce it

Load MODIS NDVI data as an Earth Engine ee.ImageCollection.
Set the region reduce scale argument to 1000 (meters).
Get the data from Earth Engine to the Colab VM as a list of lists.
End of explanation
"""


cols = {'Year': int, 'DOY': int, 'NDVI': float}
snNdviDf = eeList2Df(snNdviList, cols)
snNdviDf.head(5)

"""
Explanation: Convert the list to a tidy dataframe and preview it.
End of explanation
"""


snNdviDf['NDVI'] = snNdviDf['NDVI']/10000
snNdviDf.head(5)

"""
Explanation: Remove the NDVI scaling and preview again.
End of explanation
"""


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)

"""
Explanation: DOY line chart
End of explanation
"""


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

"""
Explanation: Reduce the plot to median and interquartile range.
End of explanation
"""


ndvi = ee.ImageCollection('COPERNICUS/S2').filterDate('2017-01-01', '2020-03-25').select('NDVI')

reReArgs['scale'] = 1000
snNdviList = getReReList(ndvi, ['Year', 'DOY', 'NDVI'])

ndviDoy = [224, 272] 

snNdviDfSub = snNdviDf[
  (snNdviDf['DOY'] >= ndviDoy[0]) & (snNdviDf['DOY'] <= ndviDoy[1])
]

snNdviDfSub = snNdviDfSub.groupby('Year').agg('min')

"""
Explanation: Relationship between PDSI and NDVI
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. 

Subset the data in the above plot where summer NDVI declines and then increases. Use the cursor tooltip to identify start and end day-of-year for this period.
Reduce the values within the year by the minimum observation.
End of explanation
"""


pdsiDoy = [1, 272]

snPdsiDfSub = snPdsiDf[
  (snPdsiDf['DOY'] >= pdsiDoy[0]) & (snPdsiDf['DOY'] <= pdsiDoy[1])
]

snPdsiDfSub = snPdsiDfSub.groupby('Year').agg('mean')

"""
Explanation: Subset the PDSI data to include all days prior to the end of the above defined NDVI DOY range.
Reduce the values within a year by mean of observations.
End of explanation
"""


ndviPdsiCombo = pd.merge(snNdviDfSub, snPdsiDfSub, how='left', on='Year') \
  .drop(columns=['DOY_x', 'DOY_y', 'Month']) \
  .reset_index()

ndviPdsiCombo.head(5)

"""
Explanation: Perform a join on 'Year' to combine the two dataframes.
End of explanation
"""


ndviPdsiCombo['Fit'] = np.poly1d(np.polyfit(ndviPdsiCombo['PDSI'], ndviPdsiCombo['NDVI'], 1))(ndviPdsiCombo['PDSI'])
ndviPdsiCombo.head(5)

"""
Explanation: Add a line of best fit between PDSI and NDVI. Including a line of best fit can be a helpful visual aid. Here, a 1D polynomial is fit throught xy point cloud defined by corresponding NDVI and PDSI observations. The resulting
fit is added to the dataframe as a new column Fit.
End of explanation
"""


# 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 

"""
Explanation: 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.
End of explanation
"""


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)

"""
Explanation: 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.
Patch-level vegetation change
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.
Find a point of interest
Zoom to a location and click the map to list the latitude and longitude.
End of explanation
"""


# Start and end dates.
startDay = 1
endDay = 365

"""
Explanation: Prepare a Landsat surface reflectance collection
Edit date window inputs based on NDVI curve plotted previously and set latitude and longitude variables from the map above.
End of explanation
"""


# 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)

"""
Explanation: Prepare a Landsat surface reflectance collection 2017-present.
End of explanation
"""


# 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'])

"""
Explanation: Reduce the region.
End of explanation
"""


cols = {'Date': str, 'MSAVI2': float}
pointNDVIDf = eeList2Df(pointNDVIList, cols)
pointNDVIDf.tail(5)

"""
Explanation: Transform and tidy the table.
End of explanation
"""


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)

"""
Explanation: Display intra-annual reduction a line plot.
End of explanation
"""


