Weeks/Week10/W10-02-Remote Sensing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "smart-bracket",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction to Remote Sensing with Python\n",
    "\n",
    ">\"...the time is now right and urgent to apply space technology towards the solution of many pressing natural resources problems being compounded by population and industrial growth.\"\n",
    "\n",
    "- [Stewart Udall, Secretary of the Interior, September 21, 1966](https://prd-wret.s3-us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/1966.09.21-DOI-Udall-Earth%20Resources%20Studied%20From%20Space.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "continuous-screen",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A remote sensing crash course\n",
    "\n",
    "<img src=\"images/remote.png\" width=600>\n",
    "\n",
    "https://youtu.be/DGE-N8_LQBo\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "modified-delay",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Landsat\n",
    "\n",
    "Landsat is the first \"civilian Earth observation satellite\" launched in 1972, a collaboration between the Department of the Interior, NASA, and the Department of Agriculture.\n",
    "\n",
    "At over 40 years, the Landsat series of satellites provides the longest temporal record of moderate resolution multispectral data of the Earth’s surface on a global basis. The Landsat record has remained remarkably unbroken, proving a unique resource to assist a broad range of specialists in managing the world’s food, water, forests, and other natural resources for a growing world population. It is a record unmatched in quality, detail, coverage, and value. Source: USGS\n",
    "\n",
    "- [NASA's Landsat page](https://www.nasa.gov/mission_pages/landsat/overview/index.html)\n",
    "- [USGS's Landsat page](https://www.usgs.gov/core-science-systems/nli/landsat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "musical-edgar",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://prd-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/styles/full_width/public/thumbnails/image/5262-Landsat-timeline-FNL-2.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3601aeb6",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Remote sensing applications in urban planning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c1100a-c44f-41c1-90bc-66e3cf630331",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- spatial and temporal analysis\n",
    "- urban development\n",
    "- urban sprawl\n",
    "- agriculture, forestry\n",
    "- human impact on climate change\n",
    "- emergency response, disaster relief"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "particular-particular",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How to access landsat data\n",
    "\n",
    "- [USGS Earth Explorer](https://earthexplorer.usgs.gov/)\n",
    "- [Amazon AWS Open Data](https://registry.opendata.aws/landsat-8/)\n",
    "- [Google Earth Engine](https://developers.google.com/earth-engine/datasets/catalog/landsat-8)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ethical-submission",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Finding satellite image of the Bobcat Fire\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/c/c1/Bobcat_Fire%2C_Los_Angeles%2C_San_Gabriel_Mountains.jpg/800px-Bobcat_Fire%2C_Los_Angeles%2C_San_Gabriel_Mountains.jpg\">\n",
    "\n",
    "The Bobcat Fire was a fire that started on September 6, 2020 as part of the 2020 California wildfire season. By December 18, it was fully contained and had burned 115,796 acres (46,861 ha) in the central San Gabriel Mountains, in and around the Angeles National Forest. It was one of the largest fires on record in Los Angeles County to date.\n",
    "\n",
    "- https://en.wikipedia.org/wiki/Bobcat_Fire"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "advised-victoria",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In this lab, we will do the following:\n",
    "\n",
    "- Use Google Earth Engine's Python API\n",
    "- define an AOI (area of interest) in the Central San Gabriel Mountains\n",
    "- import multiple Landsat images from the dates before and after the Bobcat fire\n",
    "- determine which images are best for analysis\n",
    "- create various NDVI map outputs to assess the amount of fire damage"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "retired-flash",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "racial-likelihood",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# the regulars\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# work with dates\n",
    "from datetime import datetime as dt\n",
    "\n",
    "# allow images to display\n",
    "from IPython.display import Image\n",
    "\n",
    "# earth engine\n",
    "import ee"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "otherwise-figure",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Authenticate Earth Engine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "existing-generic",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Trigger the authentication flow.\n",
    "ee.Authenticate()\n",
    "\n",
    "# Initialize the library.\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "returning-bracket",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For this lab, we will use Google Earth Engine's \"USGS Landsat 8 Surface Reflectance Tier 1\"\n",
    "\n",
    "- [EE Landsat 8](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_SR)\n",
    "\n",
    "You can also:\n",
    "\n",
    "- `.filterBounds()` method allows you to filter by location\n",
    "- `.filterDate()` allows you to filter by date"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "seventh-theme",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Define filters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "legitimate-phrase",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# coordinates of the bobcat fire\n",
    "lat = 34.289\n",
    "lon = -117.977\n",
    "\n",
    "# point of interest as an ee.Geometry\n",
    "poi = ee.Geometry.Point(lon,lat)\n",
    "\n",
    "# start date of range to filter for\n",
    "start_date = '2020-08-01'\n",
    "\n",
    "# end date\n",
    "end_date = '2020-12-31'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "reverse-immune",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Get Landsat 8 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "portuguese-asset",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# get the satellite data\n",
    "# ee.ImageCollection(\"LANDSAT/LC08/C02/T1_L2\")\n",
    "# landsat = ee.ImageCollection(\"LANDSAT/LC08/C02/T1_L2\")\\\n",
    "landsat = ee.ImageCollection(\"LANDSAT/LC08/C01/T1_SR\")\\\n",
    "            .filterBounds(poi)\\\n",
    "            .filterDate(start_date,end_date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "invisible-amplifier",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# how many images did we get?\n",
    "print('Total number:', landsat.size().getInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "jewish-dakota",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# information about the first image in our collection\n",
    "landsat.first().getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "velvet-transport",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# what about cloud cover of our first image?\n",
    "landsat.first().get('CLOUD_COVER').getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "seventh-tucson",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# when was this image taken?\n",
    "date = ee.Date(landsat.first().get('system:time_start')).getInfo()['value']\n",
    "time = date/1000\n",
    "date = dt.utcfromtimestamp(time).strftime('%Y-%m-%d %H:%M:%S')\n",
    "date"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ranking-architecture",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bands\n",
    "\n",
    "<img src=\"https://prd-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/thumbnails/image/Landsat%208%20band%20designations.jpg\" width=800>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "surprising-bleeding",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# what bands did we get?\n",
    "landsat.first().bandNames().getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "inclusive-blame",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# put the images in a list\n",
    "landsat_list = landsat.toList(landsat.size());"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "compact-childhood",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Display satellite image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "recreational-feeding",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set some parameters for image\n",
    "parameters = {\n",
    "                'min': 0,\n",
    "                'max': 2500,\n",
    "                'dimensions': 600,\n",
    "                'bands': ['B4', 'B3', 'B2']\n",
    "             }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "subject-mongolia",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# loop through each image and display it\n",
    "for i in range(landsat.size().getInfo()):\n",
    "\n",
    "    # when was this image taken?\n",
    "    date = ee.Date(ee.Image(landsat_list.get(i)).get('system:time_start'))\n",
    "    time = date.getInfo()['value']/1000\n",
    "    date = dt.utcfromtimestamp(time).strftime('%Y-%m-%d %H:%M:%S')\n",
    "    \n",
    "    # cloud cover\n",
    "    cloud = ee.Image(landsat_list.get(i)).get('CLOUD_COVER').getInfo()\n",
    "    \n",
    "    print('Image #',i,date,'Cloud cover:',cloud)\n",
    "    \n",
    "    display(Image(url = ee.Image(landsat_list.get(i)).getThumbUrl(parameters)))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "parliamentary-correspondence",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Selecting images, zooming in\n",
    "Now that we have inspected our collection of images, we can pick and choose which ones are relevant for our study. Ideally, we want to have images for before and after the fire to be able to assess the level of damage.\n",
    "\n",
    "We also want to create an ROI (region of interest) and zoom in to the area of interest. We do so by appying a 20km buffer around our POI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "blank-professional",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# create a list of images we want (before, during, after)\n",
    "landsat_sequence = [1,3,4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acceptable-tribe",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Define a region of interest with a buffer zone of 20 km\n",
    "roi = poi.buffer(20000) # meters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "played-definition",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "parameters = {'min': 0,\n",
    "              'max': 2500,\n",
    "              'dimensions': 512,\n",
    "              'bands': ['B4', 'B3', 'B2'],\n",
    "              'region':roi\n",
    "              }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "graduate-consent",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "for i in landsat_sequence:\n",
    "    # when was this image taken?\n",
    "    date = ee.Date(ee.Image(landsat_list.get(i)).get('system:time_start'))\n",
    "    time = date.getInfo()['value']/1000\n",
    "    date = dt.utcfromtimestamp(time).strftime('%Y-%m-%d %H:%M:%S')\n",
    "\n",
    "    # cloud cover\n",
    "    cloud = ee.Image(landsat_list.get(i)).get('CLOUD_COVER').getInfo()\n",
    "    \n",
    "    print('Image #',i,date,'Cloud cover:',cloud)\n",
    "    \n",
    "    display(Image(url = ee.Image(landsat_list.get(i)).getThumbUrl(parameters)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abstract-sound",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "Take a moment to inspect the three images above. What do they tell you? What do they NOT tell you?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "increasing-rebel",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Normalized Difference Vegetation Index (NDVI)\n",
    "\n",
    "The normalized difference vegetation index (NDVI) is a simple graphical indicator that can be used to analyze remote sensing measurements, often from a space platform, assessing whether or not the target being observed contains live green vegetation.\n",
    "\n",
    "<img src=\"https://eos.com/wp-content/uploads/2020/11/how-NDVI-works-EN.jpg\" width=800>\n",
    "\n",
    "- [EOS](https://eos.com/ndvi/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "surprised-elder",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# ndvi palette: red is low, green is high vegetation\n",
    "palette = ['red', 'yellow', 'green']\n",
    "\n",
    "ndvi_parameters = {'min': 0,\n",
    "                   'max': 1,\n",
    "                   'dimensions': 512,\n",
    "                   'palette': palette,\n",
    "                   'region': roi}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "reasonable-stroke",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "for i in landsat_sequence:\n",
    "\n",
    "    # when was this image taken?\n",
    "    date = ee.Date(ee.Image(landsat_list.get(i)).get('system:time_start'))\n",
    "    time = date.getInfo()['value']/1000\n",
    "    date = dt.utcfromtimestamp(time).strftime('%Y-%m-%d %H:%M:%S')\n",
    "    \n",
    "    # print some information\n",
    "    print('Image #',i,date)\n",
    "    \n",
    "    # display the image\n",
    "    display(Image(url=ee.Image(landsat_list.get(i)).normalizedDifference(['B5', 'B4']).getThumbUrl(ndvi_parameters)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exceptional-elimination",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Folium\n",
    "\n",
    "Of course, we can't end the lab without an interactive map."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9e946d8",
   "metadata": {},
   "source": [
    "Google earth engine works with folium:\n",
    "- https://developers.google.com/earth-engine/tutorials/community/intro-to-python-api-guiattard?hl=en#interactive_mapping_using_folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "expensive-replica",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hearing-shakespeare",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Google function that allows ee layers on folium\n",
    "def add_ee_layer(self, ee_image_object, vis_params, name):\n",
    "    \"\"\"Adds a method for displaying Earth Engine image tiles to folium map.\"\"\"\n",
    "    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)\n",
    "    folium.raster_layers.TileLayer(\n",
    "        tiles=map_id_dict['tile_fetcher'].url_format,\n",
    "        attr='Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
    "        name=name,\n",
    "        overlay=True,\n",
    "        control=True\n",
    "    ).add_to(self)\n",
    "\n",
    "# Add Earth Engine drawing method to folium\n",
    "folium.Map.add_ee_layer = add_ee_layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "breeding-composer",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Create a map\n",
    "my_map = folium.Map(location=[lat, lon], zoom_start=10)\n",
    "\n",
    "# Add a layer for each satellite image of interest (before, during and after)\n",
    "for i in landsat_sequence:\n",
    "\n",
    "    # when was this image taken?\n",
    "    date = ee.Date(ee.Image(landsat_list.get(i)).get('system:time_start'))\n",
    "    time = date.getInfo()['value']/1000\n",
    "    date = dt.utcfromtimestamp(time).strftime('%Y-%m-%d %H:%M:%S')\n",
    "\n",
    "    my_map.add_ee_layer(ee.Image(landsat_list.get(i)).normalizedDifference(['B5', 'B4']), \n",
    "                        ndvi_parameters, \n",
    "                        name=date)\n",
    "    \n",
    "# Add a layer control panel to the map\n",
    "folium.LayerControl(collapsed = False).add_to(my_map)\n",
    "\n",
    "# Display the map.\n",
    "display(my_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "surgical-schema",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Save the folium map as an html file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "similar-effects",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_map.save('bobcat.html')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "specialized-dallas",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Resources\n",
    "\n",
    "- [Google's Python Tutorial](https://developers.google.com/earth-engine/tutorials/community/intro-to-python-api-guiattard)\n",
    "- [Earth Lab](https://www.earthdatascience.org/courses/use-data-open-source-python/multispectral-remote-sensing/landsat-in-Python/)\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "271.188px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}