{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using AWS Lambda and PyWren for Landsat 8 Time Series\n",
    "This notebook is a simple demonstration of drilling a timeseries of NDVI values from the [Landsat 8 scenes held on AWS](https://landsatonaws.com/)\n",
    "\n",
    "### Credits\n",
    "- NDVI PyWren - [Peter Scarth](mailto:p.scarth@uq.edu.au?subject=AWS%20Lambda%20and%20PyWren) (Joint Remote Sensing Research Program)\n",
    "- [RemotePixel](https://github.com/RemotePixel/remotepixel-api) - Landsat 8 NDVI GeoTIFF parsing function\n",
    "- [PyWren](https://github.com/pywren/pywren) - Project by BCCI and riselab. Makes it easy to executive massive parallel map queries across [AWS Lambda](https://aws.amazon.com/lambda/)\n",
    "\n",
    "#### Additional notes\n",
    "The below remotely executed function will deliver results usually in under a minute for the full timeseries of more than 100 images, and we can simply plot the resulting timeseries or do further analysis. BUT, the points may well be cloud or cloud shadow contaminated. We haven’t done any cloud masking to the imagery, but we do have the scene metadata on the probable amount of cloud across the entire scene. We use this to weight a [smoothing spline](https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.interpolate.UnivariateSpline.html), such that an observation with no reported cloud over the scene has full weight, and an observation with a reported 100% of the scene with cloud has zero weight.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step by Step instructions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup Logging (optional)\n",
    "Only activate the below lines if you want to see all debug messages from PyWren. _Note: The output will be rather chatty and lengthy._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.INFO)\n",
    "%env PYWREN_LOGLEVEL=INFO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup all the necessary libraries\n",
    "This will setup all the necessary libraries to properly display our results and it also imports the library that allows us to query Landsat 8 data from the [AWS Public Dataset](https://aws.amazon.com/public-datasets/landsat/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import requests, json, numpy, datetime, os, boto3\n",
    "from IPython.display import HTML, display, Image\n",
    "import matplotlib.pyplot as plt\n",
    "import l8_ndvi\n",
    "from scipy.interpolate import UnivariateSpline\n",
    "import pywren\n",
    "\n",
    "# Function to return a Landsat 8 scene list given a Longitude,Latitude string\n",
    "# This uses the amazing developmentseed Satellite API\n",
    "# https://github.com/sat-utils/sat-api\n",
    "def getSceneList(lonLat):\n",
    "    scenes=[]\n",
    "    url = \"https://api.developmentseed.org/satellites/landsat\"\n",
    "    params = dict(\n",
    "        contains=lonLat,\n",
    "        satellite_name=\"landsat-8\",\n",
    "        limit=\"1000\") \n",
    "    # Call the API to grab the scene metadata\n",
    "    sceneMetaData = json.loads(requests.get(url=url, params=params).content)\n",
    "    # Parse the metadata\n",
    "    for record in sceneMetaData[\"results\"]:\n",
    "        scene = str(record['aws_index'].split('/')[-2]) \n",
    "        # This is a bit of a hack to get around some versioning problem on the API :(\n",
    "        # Related to this issue https://github.com/sat-utils/sat-api/issues/18 \n",
    "        if scene[-2:] == '01':\n",
    "            scene = scene[:-2] + '00'\n",
    "        if scene[-2:] == '02':\n",
    "            scene = scene[:-2] + '00'\n",
    "        if scene[-2:] == '03':\n",
    "            scene = scene[:-2] + '02'\n",
    "        scenes.append(scene)   \n",
    "    return scenes\n",
    "\n",
    "\n",
    "# Function to call a AWS Lambda function to drill a single pixel and compute the NDVI\n",
    "def getNDVI(scene):\n",
    "    return l8_ndvi.point(scene, eval(lonLat))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the code locally over a point of interest\n",
    "Let's have a look at Hong Kong, an urban area with some country parks surrounding the city: [114.1095,22.3964](https://goo.gl/maps/PhDLAdLbiQT2)\n",
    "\n",
    "First we need to retrieve the available Landsat 8 scenes from the point of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lonLat = '114.1095,22.3964'\n",
    "scenesHK = getSceneList('114.1095,22.3964')\n",
    "#print(scenesHK)\n",
    "display(HTML('Total scenes: <b>' + str(len(scenesHK)) + '</b>'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's find out the NDVI and the amount of clouds on a specific scene locally on our machine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lonLat = '114.1095,22.3964'\n",
    "thumbnail = l8_ndvi.thumb('LC08_L1TP_121045_20170829_20170914_01_T1', eval(lonLat))\n",
    "display(Image(url=thumbnail, format='jpg'))\n",
    "result = getNDVI('LC08_L1TP_121045_20170829_20170914_01_T1')\n",
    "#display(result)\n",
    "display(HTML('<b>Date:</b> '+result['date']))\n",
    "display(HTML('<b>Amount of clouds:</b> '+str(result['cloud'])+'%'))\n",
    "display(HTML('<b>NDVI:</b> '+str(result['ndvi'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, time to try this with an observation on a cloudier day. Please note that the NDVI drops too, as we are not able to actually receive much data fom the land surface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "lonLat = '114.1095,22.3964'\n",
    "thumbnail = l8_ndvi.thumb('LC08_L1GT_122044_20171108_20171108_01_RT', eval(lonLat))\n",
    "display(Image(url=thumbnail, format='jpg'))\n",
    "result = getNDVI('LC08_L1GT_122044_20171108_20171108_01_RT')\n",
    "#display(result)\n",
    "display(HTML('<b>Date:</b> '+result['date']))\n",
    "display(HTML('<b>Amount of clouds:</b> '+str(result['cloud'])+'%'))\n",
    "display(HTML('<b>NDVI:</b> '+str(result['ndvi'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Massively Parallel calculation with PyWren\n",
    "\n",
    "Now let's try this with multiple scenes and send it to PyWren, however to accomplish this we need to change our PyWren AWS Lambda function to include the necessary libraries such as rasterio and GDAL. Since those libraries are compiled C code, PyWren will not be able to pickle it up and send it to the Lambda function. Hence we will update the entire PyWren function to include the necessary binaries that have been compiled on an Amazon EC2 instance with Amazon Linux. We pre-packaged this and made it available via https://s3-us-west-2.amazonaws.com/pywren-workshop/lambda_function.zip\n",
    "\n",
    "You can simple push this code to your PyWren AWS Lambda function with below command, assuming you named the function with the default name pywren_1 and region us-west-2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambdaclient = boto3.client('lambda', 'us-west-2')\n",
    "\n",
    "response = lambdaclient.update_function_code(\n",
    "    FunctionName='pywren_1',\n",
    "    Publish=True,\n",
    "    S3Bucket='pywren-workshop',\n",
    "    S3Key='lambda_function.zip'\n",
    ")\n",
    "\n",
    "response = lambdaclient.update_function_configuration(\n",
    "    FunctionName='pywren_1',\n",
    "    Environment={\n",
    "        'Variables': {\n",
    "            'GDAL_DATA': '/var/task/lib/gdal'\n",
    "        }\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you look at the list of available scenes, we have a rather large amount. This is a good use-case for PyWren as it will allows us to have AWS Lambda perform the calculation of NDVI and clouds for us - furthermore it will have a faster connectivity to read and write from Amazon S3. If you want to know more details about the calculation, have a look at [l8_ndvi.py](/edit/Lab-4-Landsat-NDVI/l8_ndvi.py).\n",
    "\n",
    "Ok let's try this on the latest 200 collected Landsat 8 images GeoTIFFs of Hong Kong:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lonLat = '114.1095,22.3964'\n",
    "pwex = pywren.default_executor()\n",
    "resultsHK = pywren.get_all_results(pwex.map(getNDVI, scenesHK[:200]))\n",
    "display(resultsHK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display results\n",
    "Let's try to render our results in a nice HTML table first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Remove results where we couldn't retrieve data from the scene \n",
    "results = filter(None, resultsHK)\n",
    "\n",
    "#Render a nice HTML table to display result\n",
    "html = '<table><tr><td><b>Date</b></td><td><b>Clouds</b></td><td><b>NDVI</b></td></tr>'\n",
    "for x in results:    \n",
    "    html = html + '<tr>'\n",
    "    html = html + '<td>' + x['date'] + '</td>'\n",
    "    html = html + '<td>' + str(x['cloud']) + '%</td>'\n",
    "    html = html + '<td '\n",
    "    if (x['ndvi'] > 0.5):\n",
    "        html = html + ' bgcolor=\"#00FF00\">'\n",
    "    elif (x['ndvi'] > 0.1):\n",
    "        html = html + ' bgcolor=\"#FFFF00\">'\n",
    "    else:\n",
    "        html = html + ' bgcolor=\"#FF0000\">'\n",
    "    html = html + str(round(abs(x['ndvi']),2)) + '</td>'\n",
    "    html = html + '</tr>'\n",
    "html = html + '</table>'\n",
    "display(HTML(html))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This provides us a good overview but would quickly become difficult to read as the datapoints expand - let's use [Matplotlib](https://matplotlib.org/) instead to plot this out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeSeries = filter(None,resultsHK)\n",
    "\n",
    "# Extract the data trom the list of results\n",
    "timeStamps = [datetime.datetime.strptime(obs['date'],'%Y-%m-%d') for obs in timeSeries if 'date' in obs]\n",
    "ndviSeries = [obs['ndvi'] for obs in timeSeries if 'ndvi' in obs]\n",
    "cloudSeries = [obs['cloud']/100 for obs in timeSeries if 'cloud' in obs]\n",
    "\n",
    "# Create a time variable as the x axis to fit the observations\n",
    "# First we convert to seconds\n",
    "timeSecs = numpy.array([(obsTime-datetime.datetime(1970,1,1)).total_seconds() for obsTime in timeStamps])\n",
    "# And then normalise from 0 to 1 to avoid any numerical issues in the fitting\n",
    "fitTime = ((timeSecs-numpy.min(timeSecs))/(numpy.max(timeSecs)-numpy.min(timeSecs)))\n",
    "\n",
    "# Smooth the data by fitting a spline weighted by cloud amount\n",
    "smoothedNDVI=UnivariateSpline(\n",
    "    fitTime[numpy.argsort(fitTime)],\n",
    "    numpy.array(ndviSeries)[numpy.argsort(fitTime)],\n",
    "    w=(1.0-numpy.array(cloudSeries)[numpy.argsort(fitTime)])**2.0,\n",
    "    k=2,\n",
    "    s=0.1)(fitTime)\n",
    "\n",
    "fig = plt.figure(figsize=(16,10))\n",
    "plt.plot(timeStamps,ndviSeries, 'gx',label='Raw NDVI Data')\n",
    "plt.plot(timeStamps,ndviSeries, 'y:', linewidth=1)\n",
    "plt.plot(timeStamps,cloudSeries, 'b.', linewidth=1,label='Scene Cloud Percent')\n",
    "plt.plot(timeStamps,cloudSeries, 'b:', linewidth=1)\n",
    "#plt.plot(timeStamps,smoothedNDVI, 'r--', linewidth=3,label='Cloudfree Weighted Spline')\n",
    "plt.xlabel('Date', fontsize=16)\n",
    "plt.ylabel('NDVI', fontsize=16)\n",
    "plt.title('AWS Lambda Landsat 8 NDVI Drill (Hong Kong)', fontsize=20)\n",
    "plt.grid(True)\n",
    "plt.ylim([-.1,1.0])\n",
    "plt.legend(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the code over another location\n",
    "This test site is a cotton farming area in Queensland, Australia [147.870599,-28.744617](https://goo.gl/maps/GF5szf7vZo82)\n",
    "\n",
    "Let's first acquire some scenes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lonLat = '147.870599,-28.744617'\n",
    "scenesQLD = getSceneList(lonLat)\n",
    "#print(scenesQLD)\n",
    "display(HTML('Total scenes: <b>' + str(len(scenesQLD)) + '</b>'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first have a look at an individual observation first on our local machine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thumbnail = l8_ndvi.thumb('LC80920802017118LGN00', eval(lonLat))\n",
    "display(Image(url=thumbnail, format='jpg'))\n",
    "result = getNDVI('LC80920802017118LGN00')\n",
    "#display(result)\n",
    "display(HTML('<b>Date:</b> '+result['date']))\n",
    "display(HTML('<b>Amount of clouds:</b> '+str(result['cloud'])+'%'))\n",
    "display(HTML('<b>NDVI:</b> '+str(result['ndvi'])))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pywren Time\n",
    "Let's process this across all of the observations in parallel using AWS Lambda:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pwex = pywren.default_executor()\n",
    "resultsQLD = pywren.get_all_results(pwex.map(getNDVI, scenesQLD))\n",
    "display(resultsQLD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot this out again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeSeries = filter(None,resultsQLD)\n",
    "\n",
    "# Extract the data trom the list of results\n",
    "timeStamps = [datetime.datetime.strptime(obs['date'],'%Y-%m-%d') for obs in timeSeries if 'date' in obs]\n",
    "ndviSeries = [obs['ndvi'] for obs in timeSeries if 'ndvi' in obs]\n",
    "cloudSeries = [obs['cloud']/100 for obs in timeSeries if 'cloud' in obs]\n",
    "\n",
    "# Create a time variable as the x axis to fit the observations\n",
    "# First we convert to seconds\n",
    "timeSecs = numpy.array([(obsTime-datetime.datetime(1970,1,1)).total_seconds() for obsTime in timeStamps])\n",
    "# And then normalise from 0 to 1 to avoid any numerical issues in the fitting\n",
    "fitTime = ((timeSecs-numpy.min(timeSecs))/(numpy.max(timeSecs)-numpy.min(timeSecs)))\n",
    "\n",
    "# Smooth the data by fitting a spline weighted by cloud amount\n",
    "smoothedNDVI=UnivariateSpline(\n",
    "    fitTime[numpy.argsort(fitTime)],\n",
    "    numpy.array(ndviSeries)[numpy.argsort(fitTime)],\n",
    "    w=(1.0-numpy.array(cloudSeries)[numpy.argsort(fitTime)])**2.0,\n",
    "    k=2,\n",
    "    s=0.1)(fitTime)\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(16,10))\n",
    "plt.plot(timeStamps,ndviSeries, 'gx',label='Raw NDVI Data')\n",
    "plt.plot(timeStamps,ndviSeries, 'g:', linewidth=1)\n",
    "plt.plot(timeStamps,cloudSeries, 'b.', linewidth=1,label='Scene Cloud Percent')\n",
    "plt.plot(timeStamps,smoothedNDVI, 'r--', linewidth=3,label='Cloudfree Weighted Spline')\n",
    "plt.xlabel('Date', fontsize=16)\n",
    "plt.ylabel('NDVI', fontsize=16)\n",
    "plt.title('AWS Lambda Landsat 8 NDVI Drill (Cotton Farm QLD, Australia)', fontsize=20)\n",
    "plt.grid(True)\n",
    "plt.ylim([-.1,1.0])\n",
    "plt.legend(fontsize=14)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}