{ "cells": [ { "attachments": { "ndvi-formula-300x123.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAASwAAAB7CAAAAADr6uJSAAAVyklEQVR42uU9CXgUVdKVkwg5CBDCIVcA4edYQFEUkiXyE388QGX9/SPKgsEouIZ/lYgKri7rgforakDArAa8uGRZBUXYcAniciggVyDcIeQgZJLJOed7/3uv75npnqbTMyRame+l503N6341VfWq6lV1A0KY/JGWNeQfAb5FWq1BQEFtEWKTMqsFxM1epA/iiSZg+X6xj6++DS7wV2peC+IbnkAY8TTjD3wzl9EfO7i8hUx+ARbIg0WWwULrf+bNWxLNHhFkOosjl0wO1an82wTgCcT/6AoIkCIJ2suYXtVoQZQ45FtB+n79RjlLXAd58skVPOIl32dr7BXc5VA+GxNakKwEr1VPc86ahNRojX/T6PnMawUx9EkDdQvNoOAHmbXka5UZLSC5eYV8WOqqq7IxIQymgg+EUephyelpW4ZCNt0oFaWQ02NyudESKWMaEgWVs3T+8PpbkJmiGF+VK2fgjNfAPTRVwYuuMc8wMk2uscoYX8yDCabrLCwKiLdK96M+Da7kQY3SmCqHoAzG+PDbVXxo41EHmace+NZUqxREL0pUU3oCC00IaQWvNZ1VQUPiAhEpxcFeEM0N0ch1lYpBpcJxhkJMZkcCNOOy5psOWDKqfK2Hvlvjv3MwqWUya4kxeCkoEFDOwiiIrdkOD2DRVZTipWIQBgXA1GrhRqn8J9C7Hhr3P4L3Ml/By0wSD+NU0zZugnho8Lqj+HQVdlWXlzv4Hoel3EL6XZbyK6THeaWchwqrQ9NAr7ssYFbWuzRoQM+FXRWny1y6dBaW8RRCHnygfjkGf20/V3Rkauo+XJObmbmH7yl+J3MR6be8mTmXvL/8XCYPM+Z8dFRrtE1ZAub/vramRP2c9FzY/s/xUy/q4S7wNlWVc1PhDmN85ccfqp4Tm2zFVyYBTKziegrSIJ18WpIMfQjO+e4gQnTaXo3RcmIlzISHK1QR6bkwPpwKC/Q4bCBuFkoWvI7VsClCqD7utkGwlptA2w+5HkosLCdW+7EM+hEijD2nPlpObMQwhpjSiWBmO9QQOWLZFrTvVqLHgpdt2Hg6PAHYv9S0i6qei5hQwk0gZNRR1s8RCzFiIUSIlVrEoHBNf4h+XX20nNgOKxni+QPzEiHpuJq7xxELFd5CWMu/eyj6hjKFLu1gqOosoy+tS0H7Bkd84uAmAK1fqMe+iDWOx3fv7gxja1RHy4lN3Mr/Pg3PtY5fxi/hTodLWtDdDidPLOx+sXXvyzo3WRG66oBeAEyHxhxIIfzEJtATBm50IkFneRKLvuy9IblEdTSeWOxStybG5tBTuy59v2LZPw5buMtwFOd/8tne0w8Bc/hOXt/ucx3uDvZhFEg5D1obJwZazWhFxUR42kKOKbE+bQ2PFJNjRizEiIVEYjH8ql4wukJ1NImzMF7bPm4J+efc8kgCYdnfzT3uJu8aN0+KB+icdRvjLOxIjZjaqGOTVR70Rd6xP3XjwVhUSH2z6WiH2FwXz1mlD0DCpzYPzsKMWBy+c3X78D/aVUdjxOJOenlqZMI35HjXEIhNeTAZWv/xJOGyncMgZtx9nVq1Ao4C78OtZ/3rLIVoIN0Og+HYnzqXfwMDt2Kes1z5/eHW8wJnsdUQsdUw9QKFc/veGBLW4R/qo1EFzzALN2S0hZvLECofC0lLD184/Nf4mNes2PIY9HrvxJkvbgaeWLuhX75f5QNSlMF7Vz4wjrRatO59GMOWQEosVD0Lwl/1YWe17k0hqVtsCPytXn20nNiwLgyzV6dW0HkPEbzV7ePziKChimy47Qg6EA3PVSNs/7IPJ4a4DDos1hF1QNhnjBRhf8G/JjCXT8iEB0tF2wcduyMktMDbzhIg5sa/uzTGkozS8A6TD1GxvRP6W+0UtibCWvcHMPxHyhvFE4E3niIhy2/YAZRJfnqzaJoQ6VVt02ESdQA5zsLORQmQUeOls9qnpaWl9gV4+LRb9t1TmxmcEntyYiNvIpijEqHL21WIGyEll8HLvWGJ7Vm4+wzFdD4pEGs0ZDb4c9cAyxU5uoptCeNZUyqtkli4+IHw+OXH5GLIVkOCWbd5DHR9p1r23Rc5LnpR7KEKnlxm+eL+MHiTkxyduVXyfwix2LkoZPFiiMUeTTvLaymRhzjUjVJjzKW12NDLpcc8sdCmjqG3b/BBLITde26FLqtl35URi+vhiYUbcrvCyOOk4+QtMPBhAf7lEEiDGLGQjFiajrRsP9Jrez0A8SwN04ERS+QshBuyIfrBFJ+mg+2zHtC3WPru7sUMdnuZDqhsOsBsct7iUTCzWIAalA33XmDXlMWvhkQMH63VZTog1RCK+ell6mlMHsQiE+0CUVHeFjzFqXyE8pH4XYeNgUPsEe0s9+5bIPpnjKy3w1Ab6zm5cKcVLYTkA3SKtY8KxGoFf3L7J5ZiQ1qvvxyIUO8CGHsMS2JIYGkkkS2l6TCOxz98a0jH79SvkBdDcuj8MBFGVWL81zbt1xDlhSyzoU0+2t065p0GhF3bhvE66zLEzNdjOiBJrSDdRDG2p6up01bD0J1IjATQnvJ7KLGwp85i+Iviwu4tUR2NEYvTsWV/CI/LsaMDA2DANyXWokXtYHIxrpwK/7GytGrPfcCvhj9D0no9poOYlSfmtYlJFWbvBGhm1v3cJo7YTjIxxM4N3cCHGDL8yvGQkGNTG40XQ3bG/I4w9AByfBADPac++1A89F7nJI7iEOj16DPJIFjwH8JNR/xb8J4RLM/Igm8aG62l0frpSsfAbKtkOtD+ilnX+TAdOPwf2sHtR9RGE1ZDCjayUDxvxbX/dyMlTUTKylrqSG9IiSBvUofzjvS4sPutuvPgZRmfeu1wQ/lZGsth/RswrpAcN2zNzeX73Sc+zd1Mjuu+zl1F3ltX5H4n4jtW5644pzbaoWWfFYne/sXc3E01GNX99NFLs19adZyxI7YdXz5nzvITO3JzKU5ht7jFSJ8jjZROjK48bEOBZW3J3tmz1ZcuGqKz26UAnd1BPxL/OSV8h92u6vG47A63FOWzc++Qva62ThrAUVfnwNy50NK47oV68uDVp6BlXwYiC74iK3RyxTXZiyVW2Fxd5SjyxF6l7eAsKypu5LV9Y0lRJfm8vpiLgV+80sjNv7a4qFIiR0NRkQVXXyyuN0Qs97qe8duvxXa0a1mXtgV6U44EciGlO1L5cvq0HTzmyWfTPyd9hx5PZzBpxiv5VRTp2Kz0NyW5/j592rd4w9THD6lIvR/LonRK+ENOUzPQ9GUxFN0Hf9KXJimFxj0JSdagiEmlHDX3/Y56XiSeLTijkQPm0hD/pWQYWCr+lo9Bv61YvhJdnX3m3nRD9H7zMvX0gv3jjoMO6a3dUZbQSaqJEAu6LHOy95RY5OOtiZEjxhMYQ8gVm0XWFVf2dYlbBLfvVB+4vUo0CH3G4LU5vfHn7RUmZyjoANe57b/YdRFLkfqlzHWgxAq5qwApOCvh6yoClSff6wutyHKOv+903RzBS14WH7cEyVwN3/Fok1OItaxBs5PZFDmiyo8pscLhOWooIkYsRDhL2AhwruoBPciWaM0QGFXMDVI7BpJOYD9iGFxamZsgp1U0QIjVZgT02+KSc5ZIiJqMyPZrCObr0Gcjd217BsMU0YkVaP6TDOp/xUUDhFjdPoqHKZUiZzFiCTPenhg5nfw7AtHzHMzImx8HW7AnZ90sg3MtvWhA2rn32gqjHlnRQ9B2rQ/OopuiPYk6J37DIHighKm4B6BXHZIRi0u0l0Fh8ERQcuBMbEEqK/RaPiixXOuSYITVlxhidx+ygU7+vwqDt9H32wbD89iLs/4ig7JgKy2TAWSZpJ45WWxr05oZGvq2mzcdZMQiODewbAP8Y1joGzbix78eGr4be4RHSLcM3F71nlXNFhwqkVKkshXGBUaOJod03udLDFEy5SyEy++EacQuPXU3PMiMJF7Bq9a0KcJpzRY2ahilski8SEuOWM7FcRGPXOGIxUwH0ebvwxGr4W0YtBe5Nye2Wmqj/VxShjDMWRnYvPdLmi+xfOaUeri4MkbgMwwq7oHrP98jF0NOxmqTYDTZu8Luf3WCXHv9azB0vyKkyw12vwwuenHvy80WCn3YDqBYaJVWKb8BRZImwu5f6WE60D2DrmH/w/IPz90PWdWXboHMMuwthorV8BqYDma+QBZv91Rb/J4KdmZC9IReXjorJ5YliWFc/xL0L9oYkbDcKSOW8GvIOetCsG9VYHY5ChZ+BGXIFEkbUETH94KwEA/TwX02GXpwvjra2Dfkyz9D6mEhKUNuOtTLwN3SKyxkRU6elgq/p0KCAUtiABS+odP6y+RW8Cx/SUX3w7TEiKdtvoh1DX3DANQbqqf1i5yFzz/AZRLQ/MzZywgsfmH0dTDwAn9d9heiAHpv4K9R0FnmlRo05WV2jbRa0YDAWcTtW9ld4CxJXY/4RpTYTSQJ6PflYtYdx1lXXzRgusYy3d1RLxoQOQvhy09EKiOlkcNf3msXMUtToM08t39iBUI2zE7lNFg0UL997SbhvAVfraWB9ZINaxl8tf9MrVsaxrVz7frzwldPfLWhxFDRgOktNrnVKhpAbrdbEFnh2M2Dhx+j6PH6VF/RgOmt+Y60LC/rqtZeozdYUQfnxRMVblRZWt4g9JSWVpFcb3d1qYUUEDgqSgWoseuwSBX4WikpjaWl9dh24myt31Q8kHJHvctRNDM8DNXhamEUpt9LqipyMjK/dfA9GRlLybZ7wxcZSysRKvpbhgBzFv7U6Pd8nviqszmYkZGPLZPuWuHUUUKHkbL6PpDZbFozfAXSCmi6Itx9SshaYamL1myaJISPjJQW4qiRaxr8/TJK/HU2VcxvAd4Xzu4vD95rq01ffpbBfSn1hflS1w7v2mkaLrRZ1Mgv1SwllxEL0cl3H3MHhSFk+sP2+VvslfjD96tqB0YsdPLGmAUN/jP/xJuFIK9z+qttv/qUI/WP34YRp1kiLMCAU7x7IeMsRCY//eQlCuc2pEDY/HrljuP5Uxa3oscDf0GD2u/OcVbDW5B2QkdiiDwnSGm2mn43SXUT+XzPuKX0OB1CQuAJh5cY0slnW3n8vQPgvysUIxz9ffiHNkUPh894wrWT4bNjJ8mokQwMp93NOIsc7h8UstzhZytMdrNbaWdaupek+Tc7UWHzxXHdaAiJiOH1SaGwhbs4lpIriWG2lccvngx3litGIJ8usXn2UOKy6TB8elxzcFXuJ1vPOrn+qr2f5X51ej0TQ1SdzVLL/bg7sow/pFc1mZ7aXTcxdIKDS/AeM/d6uLMUeRBL5CyKf24im7xsBEYsz55sKy/8ZybABJKZ4T4zn+qvhAlfsmBS6TMDyLu7/syIhd0fx3Y56sfdaSZFA8eHwwq+dCDtxxlRMW/T0kxGLA8xpDj2Nd3g8SrFCDyxlD2UWPS09lWd4alqEg94KhYGPDi+E/T4mvTXTIeIIek3RcQwnUX2iEfAqpZRNLA+CX7hSwfSCg52CRl52AdnPX70PIV9r/eD2JUO/5ylxMcN77eFJ/ddOP1VStiIMwh92T5qyt6igzOB4ywaaMpuGUUD82PAzpcOpBU4X4fImbXexIrt1p1Cl2iIebpU/O6ZqX379u0eBR37kP+fiGMq8f9CttVP3wMPF9E0TBJoehfVTQjtfo4gV9/NE6s2A1JaRNGA60mIRnzBETEOG4dC369dXquhAHEjV9dK3yWVYxK8LzcdBIin+Ah9DRHrGlkpxjgY23BwGLzHEFfxxMLPhIK7JRQNkHTudO6YEguhdfEw8YzXatjrjrvuumNwFNx7TP7d88+kpKTcFAN9ksn/tWI/j582MAomHaO9tiXQ+XmuyGcc2cTLv4FkZtD+ozyJ0ZIEsLSEogFKLO6YEQtbpkGbJY2+FHzjj5OiWs+2yL5rv1JSUrJ9OLxxgfyvE/t5BV+3i8PHdH9TguSSVd3YbhPGhbyCR2KPetGAx9349OfLGr0xmy9gxEJI5CzX9kEw6KQP04HAqf+CjjkNHgNwCt6jh1s9C+4g+MTxrp9PjIYpPLxayUhDKVvI2VkYiz1adpYyVCajpY5bBl89b/m8Fo5Y9JjjLGSd1wZm++AsirOLuHq7PUbgiaXs4ews906AkfsxspFa320lPFQ48/vBLhmxCCYRw1I/dpZ6ObSWx2MsfqmqzOxPQF/umCl4WiORCm08HWlmN5GSnrkQnnVZOQJnwXv0ZHMVJs5nIGJWJXJ+Aa3Wcxr8rX+exIeGQS6b6g5BDJ8OBZf/PHg5sTwCm5rTNnpXKB/ts2HAEiW41ZAlqHYGzxANTyx8LgXi1zoVIxRM6PqJ3TNEI7g7HD76aTA8WkozPPMTw9Jx1fjQ22iqRuNUXgzrMqC/39RuGQ/pNjxNhxXXwyHJdKBgeTjCF7FYLvYXCXDnecX3rXu3FynZQubu2D9vR/GvZEfBi6ct5VtTwxK3Y7SmfeunCiwXPxCMUpKMN93v7g5WhkX1LmsGM9NVroXko36BZTqLHO+6wae7wwrIZrSC9/wEtGTEpY43vIvde/8zHCbMfGxwBDxFAjYWYtulZU1qA7wY7h8JuXpqd5B0exWdN8Fowo17fLZXRsN4uc6iW7evqOgsYsNuHQpJp7V/GSmkQxbX7wZAUhG275jamtoNXWeV0f6yJ+PIm3ZZnBi6SWL6fn+JIYbdHQPLoUZxxRvR3S/Q4y1531n4/uK8vB2kdsj+77zttaQsen3ev+0Cfv22vLyL2ozM4fPvrfl5eZeIo3Np41tzXlq6y8L1V2x+a86C78vy8oiRi6tmsciE360wlaKBgGRdqOnPI93il9Orcjqc4uPIHNyxi/yjYTuHS4q3uRwOt58dd4YvqAyng6upczfW19vFaZJ3Nnoamr95JCVkYaPxooFAPHdHfXYzIaXsGhYN2BbCqCNNKBpowk6NgScNnGjX+WPntUs5upAaNa+m5TxpYAa5McW1KxpYBiP2NKVoIABZYloXc2BwwjLntSoaqBzY7s2GphQNBCSepT47x7EfLrqvVdGA7Yf9/iuOtYoG/FmlhlozH/vprzU9i0ajaCAQMfhgFu8EohxFtWhAOzBlyN0J7vNrA5CtrFo04O8xdAZify07XVmraMBPCqIxMfy1Fg0EIgb/qy0aCNBW2K+zaMCUmyb8VooGDKUOBeL5BC2haCDoTxpoyUUDONhPGmj2RQNSaMY7jTEguztBlkEzX6AdDtB6ZoDBrcMm8eXVCr3pRQPqTxrQXmmMPpYhqA8pMtWbBmPR0EDfpal5Wlqg0FNIb2jB4K8dfD/a3Ecje9NI5/3gm+DwBG01NNs3bDZFAwFpzVbwXnEyXc9kbTYTMN/NUG9BcX9SbwNSsxiuubcBcHfUzd7AbFgE7RWQ+2eJCf3YW9NrPKHZqCcdNKs0AKYDxh63C1YkbJkb/gh2FN7kR1OLiSG+NlkNJQ9pPs6pZa+G/w+5ThNd/V7OLQAAAABJRU5ErkJggg==" } }, "cell_type": "markdown", "metadata": {}, "source": [ "# Reinvent 2018 - ARC329 - Massively Parallel Data Processing at Scale\n", "This notebook will walk you through the session Massively Parallel Data Processing at Scale. We will use Landsat-8 satellite imagery to calculate a Normalized Differnece Vegetation Index (NDVI) for a few interesting locations around the world evaluating the GeoTIFF data across multiple spectrum bands. [Landsat 8 scenes held on AWS](https://landsatonaws.com/)\n", "\n", "### What is NDVI (Normalized Difference Vegetation Index)?\n", "\n", "Normalized Difference Vegetation Index (NDVI) quantifies vegetation by measuring the difference between near-infrared (which vegetation strongly reflects) and red light (which vegetation absorbs). NDVI always has a range from -1 to +1 with no distinct boundary for types of land cover. For example, when you have negative values, it is likely that the majority of the scene is water. In contrast, if the NVDI for a given scene calculates to a value close to +1, there is a higher possibility that the scens contains dense green vegetation. When NDVI is close to zero, there are no green leaves and it could even be an urbanized area.\n", "\n", "The formula to calculate NDVI is as follows. For more details pertaining to how we are calculating NDVI in our examples, refer to file **l8_nvdi.py** in the source provided. This file contains the functions we will be using in this session to calulate NDVI.\n", "\n", "\n", "![ndvi-formula-300x123.png](attachment:ndvi-formula-300x123.png)\n", "\n", "### Credits\n", "- PyWren-WorkShops - This Session builds upon work done in Lab-4 of the PyWren workshop\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. This project makes it easy to execute massive parallel map queries across [AWS Lambda](https://aws.amazon.com/lambda/)\n", "\n", "#### Additional notes\n", "The functions we use in this session will deliver results normally 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 have not done any cloud masking to the imagery, but we do have the scene metadata that indicates the probable amount of clouds 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 clouds over the scene has full weight, and an observation with 100% of the scene covered with clouds has zero weight. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Step by Step Instructions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Prepare the Environment\n", "This section will setup all the necessary libraries we will be using to properly display our results, and also defines some variables and functions used throughout the session. We will import a library that allows us to query Landsat 8 data from the [AWS Public Dataset](https://aws.amazon.com/public-datasets/landsat/):\n", "\n", "Let's take a couple of minutes to walk through this initial code. We will need to update a few variables with data we gathered in the PyWren setup." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "import requests, json, numpy, random, time, datetime, os, boto3, folium\n", "from IPython.display import HTML, display, Image\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sb\n", "import l8_ndvi\n", "from scipy.interpolate import UnivariateSpline\n", "import pywren\n", "\n", "#AWS Region and Pywren Lambda Function Name - Update wth your values from the PyWren Setup\n", "AWSRegion = 'PyWren Setup AWS Region'\n", "LambdaName = 'PyWren Setup Function Name'\n", "LambdaS3Bucket='PyWren Setup S3 Bucket Name'\n", "\n", "\n", "\n", "\n", "LambdaS3Key='lambda_function.zip'\n", "\n", "\n", "#The following define the locations we will be pulling Landsat-8 data for and processing\n", "ColumbiaCoffeeLabel = 'Columbia Coffee Triangle'\n", "ColumbiaCoffee = '-75.6816,5.4716'\n", "ColumbiaCoffeeEstimatedNonMPP = 0\n", "\n", "CaNapaValleyLabel = 'California Napa Valley'\n", "CaNapaValley = '-122.2868,38.2975'\n", "CaNapaValleyEstimatedNonMPP = 0\n", "\n", "DeathValleyLabel = 'Death Valley USA'\n", "DeathValley = '-116.9325,36.5322'\n", "DeathValleyEstimatedNonMPP = 0\n", "\n", "SulaVineyardLabel = 'Sula Vineyard, Maharashtra India'\n", "SulaVineyard = '73.6885,20.0068'\n", "SulaVineyardEstimatedNonMPP = 0\n", "\n", "lonLat = ''\n", "\n", "\n", "\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(locallonLat):\n", " lonLat = locallonLat\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", " #print(sceneMetaData)\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", "\n", "# Function to Display Single Scene Data and plot location on map using folium\n", "def displaySingleScene(scene,locallonLat,locationName,locationTotalSceneCount):\n", " lonLat = locallonLat\n", " localTotalTimeEstimate = 0\n", " thumbnail = l8_ndvi.thumb(scene, eval(lonLat))\n", " display(Image(url=thumbnail, format='jpg'))\n", " SceneMap = None\n", " SceneMap = folium.Map(location=[float(lonLat.split(',')[1]),float(lonLat.split(',')[0])],zoom_start=4,width=400,height=300)\n", " folium.Marker(location=[float(lonLat.split(',')[1]),float(lonLat.split(',')[0])]).add_to(SceneMap)\n", " starttime=time.time()\n", " endtime = time.time()\n", " localTotalEstimate = (endtime-starttime)*locationTotalSceneCount\n", " result = getNDVI(scene)\n", " if result != {}:\n", " display(HTML('Date: '+result['date']))\n", " display(HTML('SceneID: '+ scene))\n", " display(HTML('Location: '+ locationName))\n", " display(HTML('Amount of clouds: '+str(result['cloud'])+'%'))\n", " display(HTML('NDVI: '+str(result['ndvi'])))\n", " display(HTML('Process Time in Sec: '+str(endtime-starttime)))\n", " display(HTML('Estimated Time in Sec to process all images in for this location ' +str(localTotalEstimate) ))\n", " display(SceneMap)\n", " else:\n", " display(HTML('Image Metadata Issue Re-run Last Cell'))\n", " return localTotalEstimate\n", " \n", "\n", "#Function to process All Landsat-8 images in a given Scene using PyWren in a MPP manner. \n", "def processScenesMPP(scenes,locallonLat,locationName,localTotalEstimate):\n", " lonLat = locallonLat\n", " index = len(scenes)\n", " pwex = pywren.default_executor()\n", " starttime=time.time()\n", " results = pywren.get_all_results(pwex.map(getNDVI, scenes[:index]))\n", " endtime = time.time()\n", " EstMPPTimeSaving = localTotalEstimate - (endtime-starttime)\n", " results = filter(None, results)\n", " display(HTML('Location Name: '+locationName))\n", " display(HTML('Scenes Processed: '+str(len(scenes))))\n", " display(HTML('Scenes With Data: '+str(len(results))))\n", " display(HTML('Process Time in Sec: '+str(endtime-starttime)))\n", " display(HTML('Estimated Time Saved using MPP: '+str(EstMPPTimeSaving)))\n", " display(HTML('--------------------------------- '))\n", " return results\n", "\n", "#Function to Display results returned from our MPP process of a given location. This uses Pandas Dataframes to display\n", "#Data in a table and folium to map the location\n", "def displayMPPResults(results,locallonLat,locationName):\n", " lonLat = locallonLat\n", " series = pd.DataFrame(results)\n", " series = series[::-1]\n", " series = series[['date','cloud','ndvi']]\n", " redColors = sb.light_palette(\"blue\",as_cmap=True)\n", " greenColors = sb.light_palette(\"green\",as_cmap=True)\n", " display(HTML('Location: '+locationName))\n", " display(HTML('Average NDVI: '+str(series['ndvi'].mean())))\n", " display(HTML('Average Cloud %: '+str(series['cloud'].mean())))\n", " with pd.option_context('display.max_rows',None,'display.max_columns',None):\n", " display(series.style.background_gradient(cmap=redColors,subset=['cloud']).highlight_max(subset=['cloud']).set_caption(locationName).hide_index().background_gradient(cmap=greenColors,subset=['ndvi']).highlight_max(subset=['ndvi']))\n", " SceneMap = None\n", " SceneMap = folium.Map(location=[float(lonLat.split(',')[1]),float(lonLat.split(',')[0])],zoom_start=4,width=400,height=300)\n", " folium.Marker(location=[float(lonLat.split(',')[1]),float(lonLat.split(',')[0])]).add_to(SceneMap)\n", " display(SceneMap)\n", "\n", "# Function to Plot all NDVI results returned from a given location dataset \n", "def plotMPPResults(results,locallonLat,locationName):\n", " timeSeries = filter(None,results)\n", " #Use DataFrame to remove duplicate dates. Two datapoints with same date will cause fitTime to contain invalid datapoints for plotting\n", " dfTS = pd.DataFrame.from_dict(timeSeries)\n", " dfTSU = dfTS.drop_duplicates(subset='date',keep='last',inplace=False)\n", " timeSeries = dfTSU.to_dict('records')\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 normalize 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", " # 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(locationName, fontsize=20)\n", " plt.grid(True)\n", " plt.ylim([-.1,1.0])\n", " plt.legend(fontsize=14)\n", " plt.show()\n", " lonLat = locallonLat\n", " SceneMap = None\n", " SceneMap = folium.Map(location=[float(lonLat.split(',')[1]),float(lonLat.split(',')[0])],zoom_start=4,width=400,height=300)\n", " folium.Marker(location=[float(lonLat.split(',')[1]),float(lonLat.split(',')[0])]).add_to(SceneMap)\n", " display(SceneMap)\n", " display(HTML('--------------------------------------------------------- '))\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Get Scene Lists\n", "The First thing we are going to do is call **getSceneList** for each of our locations of interest. We will use the data in the returned lists to get specific images for processing.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "scenesColumbiaCoffee = getSceneList(ColumbiaCoffee)\n", "scenesCaNapaValley = getSceneList(CaNapaValley)\n", "scenesDeathValley = getSceneList(DeathValley)\n", "scenesSulaVineyard = getSceneList(SulaVineyard)\n", "\n", "display(HTML('Total Columbia Coffee Scenes: ' + str(len(scenesColumbiaCoffee)) + ''))\n", "display(HTML('Total Napa Valley Scenes: ' + str(len(scenesCaNapaValley)) + ''))\n", "display(HTML('Total Death Valley Scenes: ' + str(len(scenesDeathValley)) + ''))\n", "display(HTML('Total Sula Vineyard Scenes: ' + str(len(scenesSulaVineyard)) + ''))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Get familiar with the data\n", "\n", "Now that we have our scene lists, let's take a look at a single date in the first location. We will be calling **diaplaySingleScene** with a random index on our first location's scene list. Running this cell a few times will give you an idea of the differences in the image clarity and level of cloud cover over time for the same location and also a feel for how long it takes to process a single image. Later in this session we will be using PyWren and Lambda to process ALL images in a given scene at once.\n", "\n", "## Let's have a look at a Columbian Coffee Field" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lonLat=ColumbiaCoffee\n", "ColumbiaCoffeeEstimatedNonMPP = displaySingleScene(random.choice(scenesColumbiaCoffee),ColumbiaCoffee,ColumbiaCoffeeLabel,len(scenesColumbiaCoffee))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Review Results\n", "Have a look at the output above and get a feel for the type of data we will be working with. Take note of the time it took to process a single image, we have also used this time to calculate an estimate of how long it may take to process the entire image set for this scene one at a time. Notice as the cloud percentage increases, the NDVI decreases. This is because we have less light and less reflectivity/absorption from possible vegetation in the scene. As we stated above, we are not masking for clouds in our calculations so in theory there may be a scene with a very good NDVI score being skewed slightly by cloud cover." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Let's take a look at Napa Valley" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lonLat=CaNapaValley\n", "CaNapaValleyEstimatedNonMPP = displaySingleScene(random.choice(scenesCaNapaValley),CaNapaValley, CaNapaValleyLabel,len(scenesCaNapaValley))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Now Death Valley" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lonLat=DeathValley\n", "DeathValleyEstimatedNonMPP = displaySingleScene(random.choice(scenesDeathValley),DeathValley, DeathValleyLabel,len(scenesDeathValley))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice in Death Valley, we get low NDVI even when the cloud cover percentage is low. This is because there isn't much if any vegetation to report back to the LandSat Sensor Array either reflected or absorbed light levels.\n", "\n", "## And Finally Sula Vineyard, India" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "lonLat=SulaVineyard\n", "SulaVineyardEstimatedNonMPP = displaySingleScene(random.choice(scenesSulaVineyard),SulaVineyard, SulaVineyardLabel,len(scenesSulaVineyard))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Review the results from all above single images runs\n", "As you can see each of these single processes took on average 2 seconds to complete. If we had 1000 images, this would take over 30 min to process. Now that we have a good baseline for timing, let's move on to the MPP part of this session." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Massively Parallel Processing 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 copied it to your S3 bucket during the readme setup for this session.\n", "\n", "Running the cell below will update the Lambda Function that the Pywren installer created. Before running this, make sure the **LambdaName** and **AWSRegion** variables are defined in the initialization cell at the top of this notebook and match those used when configuring Pywren-Setup." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lambdaclient = boto3.client('lambda', AWSRegion)\n", "\n", "response = lambdaclient.update_function_code(\n", " FunctionName=LambdaName,\n", " Publish=True,\n", " S3Bucket=LambdaS3Bucket,\n", " S3Key=LambdaS3Key\n", ")\n", "\n", "response = lambdaclient.update_function_configuration(\n", " FunctionName=LambdaName,\n", " Environment={\n", " 'Variables': {\n", " 'GDAL_DATA': '/var/task/lib/gdal'\n", " }\n", " }\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Now let's create a CloudWatch Dashboard to watch our lambda function metrics\n", "The code below will create a Dashboard in CloudWatch in the configured region for the configured function and load the metrics for **Concurrent Executions**, **Throttles**, **Invocations**, **Duration**, and **Errors** \n", "Before continuing, please open your console http://console.aws.amazon.com and locate the newly created dashboard. \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cloudwatch = boto3.client(\"cloudwatch\")\n", "\n", "widgets = []\n", "\n", "widget = {\n", " \"type\": \"metric\",\n", " \"x\": 0,\n", " \"y\": 0,\n", " \"width\": 24,\n", " \"height\": 3,\n", " \"properties\": {\n", " \"view\": \"singleValue\",\n", " \"metrics\": [\n", " [ \"AWS/Lambda\", \"ConcurrentExecutions\" ],\n", " [ \".\", \"Throttles\", \"FunctionName\", LambdaName ],\n", " [ \".\", \"Invocations\", \"FunctionName\", LambdaName ],\n", " [ \".\", \"Duration\", \"FunctionName\", LambdaName ],\n", " [ \".\", \"Errors\", \"FunctionName\", LambdaName ]\n", " ],\n", " \"view\": \"singleValue\",\n", " \"region\": AWSRegion,\n", " \"stat\": \"Sum\",\n", " \"period\": 300,\n", " \"stacked\": False\n", " }\n", " }\n", "\n", "widgets.append(widget)\n", "\n", "body = {'widgets' : widgets}\n", "body_j = json.dumps(body)\n", "\n", "cloudwatch.put_dashboard(DashboardName = \"ARC329-Reinvent2018-\"+LambdaName,\n", " DashboardBody=body_j)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ok, we have now updated our AWS Lambda Function and created a CloudWatch Dashboard to monitor our Lambda executions let's move on\n", "As you may have noticed, each of our locations has over 100 scenes with two locations having more than 300 scenes. This is a good use-case for Pywren as it will allow us to have AWS Lambda perform the calculation of NDVI for us in parallel. Running the calculation in AWS will also be more efficiant when pulling images from the LandSat-8 data set stored in S3. \n", "\n", "Ok let's run the full scene list for the Columbian Coffee Triangle to generate our dataset for future processing." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "lonLat = ColumbiaCoffee\n", "resultsColumbiaCoffee = processScenesMPP(scenesColumbiaCoffee,ColumbiaCoffee,ColumbiaCoffeeLabel,ColumbiaCoffeeEstimatedNonMPP)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Now let's get our dataset for Napa Valley" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lonLat = CaNapaValley\n", "resultsCaNapaValley = processScenesMPP(scenesCaNapaValley,CaNapaValley,CaNapaValleyLabel,CaNapaValleyEstimatedNonMPP)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Now let's get our dataset for Death Valley\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lonLat = DeathValley\n", "resultsDeathValley = processScenesMPP(scenesDeathValley,DeathValley,DeathValleyLabel,DeathValleyEstimatedNonMPP)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Now let's get our dataset for Sula Vineyard" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "lonLat = SulaVineyard\n", "resultsSulaVineyard = processScenesMPP(scenesSulaVineyard,SulaVineyard,SulaVineyardLabel,SulaVineyardEstimatedNonMPP)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Now Let's Display Results\n", "Let's take a look at the dataset for Columbia Coffee Triangle. We are using a Pandas DataFrame with some style formatting applied to create a heatmap for each data point of interest. We use a blue pallet to indicate the level of cloud cover where a darker blue indicates a higher amount of cloud cover. We use a green pallet to indicate the NDVI level from -1 to 1 where the darker the gradient, the better the NDVI reading. Notice that where the clouds are darker blue, the corresponding NDVI is lighter green. The opposite is true as well, lighter clouds corresponds to a darker NDVI. We are also highlighting Yellow the MAX value in each column.\n", "The details of this function **displayMPPResults** are defined in the initilization cell." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "displayMPPResults(resultsColumbiaCoffee,ColumbiaCoffee,ColumbiaCoffeeLabel)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Let's take a look at Napa Valley CA" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "displayMPPResults(resultsCaNapaValley,CaNapaValley,CaNapaValleyLabel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Review Results for Death Valley USA" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "displayMPPResults(resultsDeathValley,DeathValley,DeathValleyLabel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Review Results for Sula Vineyard India" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "displayMPPResults(resultsSulaVineyard,SulaVineyard,SulaVineyardLabel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While the DataFrame table visualization is not bad, we can use the [Matplotlib](https://matplotlib.org/) library instead to plot our data over time. We call the function **plotNDVIData** which is defined in the initialization cell at the beginning of this notebook.\n", "\n", "## Let's Plot the Columbian Coffee Triangle" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotMPPResults(resultsColumbiaCoffee,ColumbiaCoffee,ColumbiaCoffeeLabel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Now Let's Plot of Napa Valley CA" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "plotMPPResults(resultsCaNapaValley,CaNapaValley,CaNapaValleyLabel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Now we will Plot of Death Valley USA" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotMPPResults(resultsDeathValley,DeathValley,DeathValleyLabel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finally let's take a look at the Plot of Sula Vineyard " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotMPPResults(resultsSulaVineyard,SulaVineyard,SulaVineyardLabel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Challenge : Define a location of your choice and plot the NDVI results" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.15" } }, "nbformat": 4, "nbformat_minor": 2 }