05b_Google_Earth_Engine.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "52786534",
   "metadata": {},
   "source": [
    "# Getting started with python API for Google Earth Engine (GEE)\n",
    "\n",
    "## What is GEE and why use it?\n",
    "(quote from the official website)\n",
    "____________________________________________________________________________________________\n",
    "Google Earth Engine is a geospatial processing service. With Earth Engine, you can perform geospatial processing at scale, powered by Google Cloud Platform. The purpose of Earth Engine is to:\n",
    "\n",
    "    * Provide an interactive platform for geospatial algorithm development at scale\n",
    "    * Enable high-impact, data-driven science\n",
    "    * Make substantive progress on global challenges that involve large geospatial datasets\n",
    "____________________________________________________________________________________________\n",
    "\n",
    "More information on Google Earth Engine can be found here: https://developers.google.com/earth-engine\n",
    "\n",
    "### A couple of important points are summarized below: \n",
    "\n",
    "__-->__ GEE hosts a ***multi-petabyte catalog of satellite imagery and other geospatial datasets*** on their servers ***(1 petabyte = 1024 terabytes, or a million gigabytes)*`. More information can be found here: https://developers.google.com/earth-engine/datasets. \n",
    "\n",
    "__-->__ GEE uses Java script API but also has a python API. Both APIs aaccess the same server-side functionality. However, there are syntax differences between the two programming langauges. More information and examples can be found here: https://developers.google.com/earth-engine/guides/python_install#syntax\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Today's exercise: Case study - Hurricane Maria (2017) impact on Puerto Rico \n",
    "\n",
    "\n",
    "_**IMPORTANT:**_ This script ported from Google Earth Engine (GEE) Java script and adaptaed to python jupyter notebook.\n",
    "\n",
    "_**--> Source & citation provided below.**_\n",
    "\n",
    "_**Source of original GEE script:**_ NASA-ARSET-training - Satellite Observations for Analyzing Natural Hazards on Small Island Nations\n",
    "https://appliedsciences.nasa.gov/join-mission/training/english/arset-satellite-observations-analyzing-natural-hazards-small-island\n",
    "\n",
    "_**Full citation:**_\n",
    "Podest, E.; McCartney, S.; Mehta, A.; Englander, J.; Hamlington, B.; Stanley, T. (2021). Satellite Observations for Analyzing Natural Hazards on Small Island Nations. NASA Applied Remote Sensing Training Program (ARSET). https://appliedsciences.nasa.gov/join-mission/training/english/arset-satellite-observations-analyzing-natural-hazards-small-island\n",
    "\n",
    "\n",
    "\n",
    "_**--> Revisions and edits by Atef AMRICHE, Jeff Liu 2022.**_\n",
    "\n",
    "**-----------------------------------------------------------------------------------**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a46d9536",
   "metadata": {},
   "source": [
    "# VERY IMPORTANT:\n",
    "____________________________________________________________________________________________\n",
    "____________________________________________________________________________________________\n",
    "General information is formatted with an arrow at the beginning. \n",
    "\n",
    "**--> general information** about a certain tool or library or dataset ... and link provided at the end where available.\n",
    "\n",
    "\n",
    "python script analysis steps are presented in bullet points. Important information is highlighted in **bold** and code is notated in `preformatted text`. example below:\n",
    "- step 1 python analysis: import study area **\"name of study area\"** shapefile.\n",
    "to do this we use the command `example python command`. Other important information about this command or type of command can be found here (*example link*)\n",
    "- step 2 python analysis: use the command `example python command` to do the following operation.\n",
    "____________________________________________________________________________________________\n",
    "____________________________________________________________________________________________"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cdda9ab",
   "metadata": {},
   "source": [
    "# ... Now let's begin !!\n",
    "\n",
    "## The following libraries are imported to use in this demo.\n",
    "\n",
    "- `ee` is the earth-engine python library that is installed & authenticated uzing a personal/professional Google Earth Engine account before use. More information can be found here: https://developers.google.com/earth-engine/guides/python_install\n",
    "\n",
    "\n",
    "- `geemap` is a python package that allows for interactive use of Google Earth Engine mapping and data analysis functions. More information can be found here: https://geemap.org/\n",
    "\n",
    "    The specific function `geemap` allows for creating interactive GEE maps that have a __base layer__ (google maps / terrain / other) plus a __specific result(s)__ overlayed on top (ex: satellite image, boundary of the study area, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f0a1698",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import necessary libraries\n",
    "import ee\n",
    "import geemap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50c08826-20a3-47a3-a322-144501742169",
   "metadata": {},
   "outputs": [],
   "source": [
    "# you only need to run the first time running the notebook, or if you haven't run it in more than a week\n",
    "# comment it out once you've run it the first time\n",
    "ee.Authenticate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6075d91-8042-4290-9f4e-2124866f7c24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize earth engine\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd3805cc",
   "metadata": {},
   "source": [
    "- __User input - event start and end dates -__ Here, we specify a date interval to search for satellite images before and after the hurricane. This must be an interval because satellite images are __NOT__ necessarily acquired on the same dates each year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7864bfae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set dates for pre- and post-hurricane event\n",
    "pre_event_beg = '2017-07-15'\n",
    "pre_event_final = '2017-07-31'\n",
    "\n",
    "post_event_beg = '2017-09-19'\n",
    "post_event_final = '2017-09-25'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "627bb3a2",
   "metadata": {},
   "source": [
    "## Load the polygon for the study area \"Puerto Rico\"\n",
    "\n",
    "__-->__ ***The FAO GAUL*** (Food and Agriculture Organization - Global Administrative Unit Layers) simplified 500m 2015 is a ***geospatial dataset*** that is part of the GEE collection. More information can be found here: https://developers.google.com/earth-engine/datasets/catalog/FAO_GAUL_2015_level0?hl=en\n",
    "\n",
    "This dataset contains country boundaries as well as other information. It can be loaded using the command `ee.FeatureCollection`. Then, we utilize the filtering function to examine the column `ADM0_NAME`, which contains the country name. We apply the GEE filter `Equal` using the command `ee.Filter.eq` to extract the polygon of our study area **\"Puerto Rico\"**. \n",
    "\n",
    "More information on GEE filters and how to use them can be found here: https://developers.google.com/s/results/earth-engine?q=ee.filter\n",
    "\n",
    "The specific information about the GEE filter \"equal\" and how to use it can be found here: https://developers.google.com/earth-engine/apidocs/ee-filter-eq?hl=en\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d9099c6-2d16-4d6d-b985-6d4ac89d9e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the shapefile for your area of interest\n",
    "roi = ee.FeatureCollection(\"FAO/GAUL_SIMPLIFIED_500m/2015/level0\").filter(ee.Filter.eq('ADM0_NAME', 'Puerto Rico'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18478cfa-efa2-45dc-89f6-21e2ae0df1b5",
   "metadata": {},
   "source": [
    "## Load a map and add our results\n",
    "Here, we specify the center of the map coordinates and the zoom level. This would vary depending on the study area location and size, and can be determined through trial and error. We use the command `geemap.Map` with the default Google maps as a base layer.\n",
    "    - For other locations, first determine the center coordinates, then\n",
    "    - Zoom less (smaller zoom number) for large areas, or zoom more (larger zoom number) for small areas\n",
    "\n",
    "Following that, we use the command `Map.addLayer` to display the `\"roi\"` layer , which corresponds to the polygon of Puerto Rico created in the previous step. This command allows for the selection of several display options. In our case, we only used the 'color' to assign the gray color to the polygon of our study area. More information about this command can be found here: https://developers.google.com/earth-engine/apidocs/map-addlayer?hl=en\n",
    "\n",
    "At the end, we simply use the command `Map` to actually display our output map in a separate window below the python code cell. Some general information about interactive maps and map formatting is provided below as it is used repeatedly throughout this exercise.\n",
    "\n",
    "\n",
    "__-->__ Each time we want to generate a map output, we first use the command `Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)`. This would create a new map with the default Google maps base layer showing our study area. Then we use the command `Map.addLayer` to specify which layers (results) to display and what display parameters are used (color, minimum and maximum values, etc). Finally we use the command `Map` to actually display our output map below the python code cell. \n",
    "\n",
    "__-->__ In each map window, we can zoom in and out as it is an interactive map. \n",
    "\n",
    "__-->__ Multiple layers or results can be displayed in a single map. In the **top right corner**, we see an icon showing multiple squares on top of each other (representing layers). We can hover over that icon to view and toggle which layers to display.\n",
    "\n",
    "__-->__ There are other buttons available on the left side of the map. Most of these we will not use in this exercise. They are:\n",
    "\n",
    "    The first two buttons are +/- to zoom in and out. The third button is to go **full screen**. \n",
    "    We can press `esc` to exit full screen.\n",
    "\n",
    "    The next five buttons **(not used in this exercise)** are to draw various shapes \n",
    "    (line, polygon, square/rectangle, circle, point, circle marker. \n",
    "\n",
    "    The next two buttons allow for **editing** and **deleting** existing shapes.\n",
    "\n",
    "    The last button is a map search button (similar to google maps search). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ccba6fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Use Google maps as the default basemap & display the study area polygon on it.\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(roi, {'color': 'gray'}, 'Study Area')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07581d51",
   "metadata": {},
   "source": [
    "### OPTIONAL - Convert Feature-Collection to GeoPandas-Dataframe\n",
    "This can be done using the command ` geemap.ee_to_geopandas`  (example results shown below for the `roi`  layer)\n",
    "\n",
    "\n",
    "__-->__ ***Using the geemap library*** - This offers may conversion commands to and other functionality to convert data between GEE variable types to common python geospatial and other variable types. Useful information can be found here https://geemap.org/common/?h=ee_to_#geemap.common.ee_num_round"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "974d9bcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpd_dataframe = geemap.ee_to_geopandas(roi, selectors=None, verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8afe0bc",
   "metadata": {},
   "source": [
    "## Load SAR images from the European satellite Sentinel-1.\n",
    "We do this separately before the event and then after the event for our study area. \n",
    "\n",
    "__-->__ ***The COPERNICUS/S1_GRD*** are ***satellite images from Sentinel-1*** that are part of the GEE collection. More information can be found here: https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD?hl=en\n",
    "\n",
    "This dataset contains historical satellite images since 2014 (ongoing satellite mission). It can be loaded using the command `ee.ImageCollection`. Then, we apply multiple filters using the command `.filter`  to select a specific SAR product, orbit, resolution, boundary, dates, and bands. More specifically:\n",
    "\n",
    "    The command \"ee.Filter.eq\" allows us to extract the \"instrument mode\" equal to \"IW\" (Interferometric Wide-swath)\n",
    "    The command \"ee.Filter.eq\" allows us to extract the \"orbit\" equal to \"Ascending\"\n",
    "    \n",
    "    The command \".filterMetadata\" allows us to extract the \"spatial resolution\" equal to \"10m\"\n",
    "    \n",
    "    The command \".filterBounds\" allows us to clip the image to our \"roi\" study area\n",
    "\n",
    "    The command \".filterDate\" allows us to extract available satellite images between a start anad end data.\n",
    "    Here, we used the dates we defined above to extract available images before and after the hurricane.\n",
    "    \n",
    "    The command \".select\" allows us to extract the \"specific bands (VV, VH)\" for this SAR product\n",
    "    \n",
    "There is a large number of filter/selection commands in GEE. More information on these commands and how to use them can be found here: https://developers.google.com/s/results/earth-engine?q=ee.filter\n",
    "\n",
    "A given study area is not necessarily covered by a single SAR image. Therefore, after extracting the available SAR images before and after the hurricane, we need to generate a single image mosaic covering our entire study area. For this, we use the command `.mosaic` and apply it to our Sentinel-1 collections before and after the hurricane. We use the command `.Clip` agaain to make have the final mosaic image clipped to our study area.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a76308e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Sentinel-1 C-band SAR Ground Range collection (log scale): pre-event\n",
    "S1collection_1 = ee.ImageCollection('COPERNICUS/S1_GRD') \\\n",
    "                   .filter(ee.Filter.eq('instrumentMode', 'IW')) \\\n",
    "                   .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING')) \\\n",
    "                   .filterMetadata('resolution_meters', 'equals' , 10) \\\n",
    "                   .filterBounds(roi) \\\n",
    "                   .filterDate(pre_event_beg, pre_event_final) \\\n",
    "                   .select('VV', 'VH')\n",
    "\n",
    "# Load Sentinel-1 C-band SAR Ground Range collection (log scale): post-event\n",
    "S1collection_2 = ee.ImageCollection('COPERNICUS/S1_GRD') \\\n",
    "                   .filter(ee.Filter.eq('instrumentMode', 'IW')) \\\n",
    "                   .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING')) \\\n",
    "                   .filterMetadata('resolution_meters', 'equals' , 10) \\\n",
    "                   .filterBounds(roi) \\\n",
    "                   .filterDate(post_event_beg, post_event_final) \\\n",
    "                   .select('VV', 'VH')\n",
    "\n",
    "# Create a mosaic for each collection\n",
    "S1_1 = S1collection_1.mosaic().clip(roi)\n",
    "S1_2 = S1collection_2.mosaic().clip(roi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5d524d7",
   "metadata": {},
   "source": [
    "## Preprocess the SAR images\n",
    "This step is common for SAR images. Basically, we apply a smoothing filter because SAR signal is grainy. In this case, we use the command `.focal_mean`  which calculates the average values of pixels within a certain radius **(in our case 30m)** and assign that value to the center pixel.\n",
    "\n",
    "### Load a map and add our results \n",
    "\n",
    "First, we specify the center of the map coordinates and the zoom level (same step as before).\n",
    "We use the command `geemap.Map` with the default Google maps as a base layer. \n",
    "\n",
    "Then, we add our SAR images (4 separate layers) using the command \".addLayer\" four times. \n",
    "- 1st argument represents the variable we want to add (ex: S1_1 for Sentinel-1 collection 1)\n",
    "- 2nd argument represents the visualization parameters (ex: minimum value = -25, maximum value = -5)\n",
    "The min/max values change from one image to another, so we choose fixed values for consistent colors.\n",
    "- 3rd argument is the name of the layer to be displayed on the map window.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "656965f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply speckle filter\n",
    "SMOOTHING_RADIUS = 30\n",
    "pre_event_filtered = S1_1.focal_mean(SMOOTHING_RADIUS, 'circle', 'meters')\n",
    "post_event_filtered = S1_2.focal_mean(SMOOTHING_RADIUS, 'circle', 'meters')\n",
    "\n",
    "# Add images to \"Layers\" in order to visualize them\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(S1_1, {'min':-25, 'max':-5}, 'S1-1')\n",
    "Map.addLayer(S1_2, {'min':-25, 'max':-5}, 'S1-2')\n",
    "Map.addLayer(pre_event_filtered, {'min':-25, 'max':-5}, 'S1-1-Filt')\n",
    "Map.addLayer(post_event_filtered, {'min':-25, 'max':-5}, 'S1-2-Filt')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e298b5a4",
   "metadata": {},
   "source": [
    "## Combine the SAR image bands to generate a 3-band image for visualization\n",
    "Here, we use the command `.cat` (short for concatenate) to combine the same SAR band (filtered VV or VH) and generate a 3-band image and visualize it in RGB colors. We use the command `.select` (same as before) to choose a specific band.\n",
    "\n",
    "\n",
    "### Load a map and add our results\n",
    "\n",
    "First, we specify the center of the map coordinates and the zoom level (same step as before).\n",
    "We use the command `geemap.Map` with the default Google maps as a base layer. \n",
    "\n",
    "Then, we add our SAR images (2 separate layers) using the command \".addLayer\" twice. \n",
    "- 1st argument represents the variable we want to add (new RGB images in this case).\n",
    "- 2nd argument represents the visualization parameters (different values for VV and VH RGB images).\n",
    "- 3rd argument is the name of the layer to be displayed on the map window.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a5fbed4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create RGB images\n",
    "S1_VV_RGB = ee.Image.cat(pre_event_filtered.select('VV'), post_event_filtered.select('VV'), pre_event_filtered.select('VV'))\n",
    "S1_VH_RGB = ee.Image.cat(pre_event_filtered.select('VH'), post_event_filtered.select('VH'), pre_event_filtered.select('VH'))\n",
    "\n",
    "# Add images to \"Layers\" in order to visualize them\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(S1_VV_RGB, {'min':-18, 'max':0}, 'S1-VV-RGB')\n",
    "Map.addLayer(S1_VH_RGB, {'min':-25, 'max':-5}, 'S1-VH-RGB')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e3b625c",
   "metadata": {},
   "source": [
    "## Compare before and after hurricane images\n",
    "We compute the ratio between the before and after images\n",
    "\n",
    "\n",
    "### Load a map and add our results \n",
    "\n",
    "First, we specify the center of the map coordinates and the zoom level (same step as before).\n",
    "We use the command `geemap.Map` with the default Google maps as a base layer. \n",
    "\n",
    "Then, we add our SAR difference image layer using the command \".addLayer\". \n",
    "- 1st argument represents the variable we want to add (new difference image in this case).\n",
    "- 2nd argument represents the visualization parameters (min & max values).\n",
    "- 3rd argument is the name of the layer to be displayed on the map window.\n",
    "\n",
    "\n",
    "## Apply a threshold to the SAR difference image \n",
    "This is also a common practice in image analysis. This way, we separate flooded areas and produce a binary image (1 and 0 values). These values represent SAR difference within our threshold interval or outside our threshold interval. We use the command `.gt`  which stands for \"Greater-Than\" and the command `.lt`  which stands for \"Less-Than\".\n",
    "\n",
    "\n",
    "### Load a map and add our results\n",
    "\n",
    "Because we are loading this result in the same map, we do NOT specify center coordinates and zoom again.\n",
    "\n",
    "Then, we add our SAR difference image layer using the command \".addLayer\". \n",
    "- 1st argument represents the variable we want to add (new difference image in this case).\n",
    "Here, we used the binary inundation image as a mask \".updateMask\" to display only the values = 1\n",
    "- 2nd argument represents the visualization parameters (color palette \"Green-Blue\" or GnBu in this case).\n",
    "- 3rd argument is the name of the layer to be displayed on the map window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5f07e67",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute ratio between before and after images\n",
    "differenceVH = post_event_filtered.select('VH').divide(pre_event_filtered.select('VH'))\n",
    "\n",
    "# Add images to \"Layers\" in order to visualize them\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(differenceVH, {'min': 0, 'max':2}, 'difference VH filtered')\n",
    "\n",
    "# Apply a threshold - based on difference image values\n",
    "UPPER_THRESHOLD = 1.15\n",
    "LOWER_THRESHOLD = 0.7\n",
    "inundation1 = differenceVH.gt(UPPER_THRESHOLD).Or(differenceVH.lt(LOWER_THRESHOLD))\n",
    "\n",
    "# Add images to \"Layers\" in order to visualize them\n",
    "Map.addLayer(inundation1.updateMask(inundation1),\n",
    "             {'palette':\"GnBu\"},'Flooded Areas - RAW')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36fb2937",
   "metadata": {},
   "source": [
    "## Calculate pixel connectivity and remove isolated pixels\n",
    "Isolated pixels are most likely outliers that should be removed. We do this in two steps:\n",
    "\n",
    "- Calculate connected pixels count using the command \".connectedPixelCount\".\n",
    "This is applied to the original inundation image.\n",
    "- Generate a new binary innundation image that excludes small patches from the previous inundation image.\n",
    "This is done by excluding pixels that are connected to less than 8 other pixels.  \n",
    "For this, we use the command \".updateMask\" and combine it with the filter \".gte\" (Greater-Than-or-Equal) \n",
    "\n",
    "\n",
    "### Load a new blank map and add our results\n",
    "Here, we display the study area polygon to have better contrast when visualizing the flooded area results.\n",
    "\n",
    "- First, we set the map center coordinates and zoom level (same step as before).\n",
    "- Add \"roi\" layer (color gray)\n",
    "- Add the new \"inundation\" layer with no isolated pixels (color palette \"GnBu\"). \n",
    "The new binary inundation image is used as a mask with the command `.updateMask` \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d7a78d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate pixel connectivity and remove those connected by less than 8 pixels.\n",
    "connections = inundation1.connectedPixelCount()\n",
    "inundation2 = inundation1.updateMask(connections.gte(8))\n",
    "\n",
    "# Add images to \"Layers\" in order to visualize them\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(roi, {'color': 'gray'}, 'Study Area')\n",
    "Map.addLayer(inundation2.updateMask(inundation2), {'palette':\"GnBu\"},'Flooded Areas - non-sparse')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "381f3193",
   "metadata": {},
   "source": [
    "## Remove misclassified pixel where the slope is higher than 5%\n",
    "Water collection in slopes higher than 5% is unlikely since it flows downhill, and means that these inundated pixels are most likely misclassified and should be removed. To do this, we need to load a new dataset: **Digital Elevation Model (DEM)** and compute the slope. Then we can use the slope to mask out areas with values higher than 5%.\n",
    "\n",
    "__-->__ ***USGS/SRTMGL1_003*** global elevation raster image is part of the GEE collection. More information can be found here: https://developers.google.com/earth-engine/datasets/catalog/USGS_SRTMGL1_003?hl=en\n",
    "\n",
    "- Load the global elevation data using the command \"ee.Image\"\n",
    "- Use the GEE algorithms \"terrain to compute slope, aspect, and hillshade using the DEM\n",
    "\n",
    "__-->__ ` ee.Algorithms`  is a collection of image analysis and data processing functions. These can be general mathematical functions or functions that are specific product (satellite image or geospatial dataset). More information about the specific **terrain** algorithm can be found here: https://developers.google.com/earth-engine/apidocs/ee-algorithms-terrain?hl=en\n",
    "\n",
    "- Select the \"slope\" layer from the terrain results using the command `.select`\n",
    "- Generate a new binary innundation image that excludes slopes > 5%.\n",
    "This is done using the command \".updateMask\" combined with the command `.lt` which stands for Less-Than. \n",
    "\n",
    "\n",
    "### Load a new blank map and add our results\n",
    "Here, we display the study area polygon to have better contrast when visualizing the flooded area results.\n",
    "\n",
    "- First, we set the map center coordinates and zoom level (same step as before).\n",
    "- Add the \"DEM\" layer (specific min/max values)\n",
    "- Add \"roi\" layer (color gray)\n",
    "- Add the new \"inundation-slope-corrected\" layer (color palette `GnBu` (Green-Blue))\n",
    "The new binary inundation image is used as a mask with the command `.updateMask` \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "302241ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove misclassified pixels in areas with slopes greater than 5%\n",
    "srtm = ee.Image('USGS/SRTMGL1_003')\n",
    "terrain = ee.Algorithms.Terrain(srtm)\n",
    "slope = terrain.select('slope')\n",
    "inundation3 = inundation2.updateMask(slope.lt(5))\n",
    "\n",
    "# Add images to \"Layers\" in order to visualize them\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(srtm, {'min':0, 'max':1000}, 'SRTM')\n",
    "Map.addLayer(roi, {'color': 'gray'}, 'Study Area')\n",
    "Map.addLayer(inundation3.updateMask(inundation3),{'palette':\"GnBu\"},'Flooded Areas - slope-adjusted')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e666ceb",
   "metadata": {},
   "source": [
    "- **Calculate inundation extent** - First, we multiply the latest binary inundation image by the pixel area value. For this, we use the command `.multiply`  and combine it with the function `ee.Image.pixelArea`  (described below).\n",
    "\n",
    "__-->__ ` ee.Image.pixelArea`  is function that Generate an image in which the value of each pixel is the area of that pixel in square meters. The returned image has a single band called **area** (https://developers.google.com/earth-engine/apidocs/ee-image-pixelarea?hl=en).\n",
    "\n",
    "\n",
    "- __Calculate the total area of inundated pixels -__ Here, we use the function `.reduceRegion`  to compute the sum (area) of all the inundated pixels.\n",
    "\n",
    "__-->__ ` ee.Image.reduceRegion`  is function that Apply a reducer to all the pixels in a specific region. A reducer is a statistical metric such as **mean, median, standard deviation, etc.** This function allows for several other options which can be seen using this link: https://developers.google.com/earth-engine/apidocs/ee-image-reduceregion?hl=en\n",
    "\n",
    "The options used in this example are:\n",
    "\n",
    "- `reducer = ee.Reducer.sum()` --> to select the Sum statistics \n",
    "- `geometry = roi` --> to limit the extent to our study area\n",
    "- `scale = 10` --> which corresponds to our SAR image native spatial resolution\n",
    "- `maxPixels = 1e9` --> this option is commented out. \n",
    "It is used to request that GEE processess a larger number of pixels.\n",
    "The request is not always possible to do because GEE has some limits on processing power and request size.\n",
    "- `bestEffort = True` --> is set to reduce the processing time.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75f182e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate inundation extent. Create a raster that contains information on pixel area\n",
    "inundation_area_pixel = inundation3.multiply(ee.Image.pixelArea())\n",
    "\n",
    "# Sum the area covered by inundated pixels\n",
    "# 'bestEffort: True' to reduce processing time. Note - for more accurate results set\n",
    "# bestEffort to False and increase 'maxPixels'.\n",
    "inundation_stats = inundation_area_pixel.reduceRegion(\n",
    "                                         reducer= ee.Reducer.sum(),\n",
    "                                         geometry= roi,\n",
    "                                         scale= 10, # native resolution\n",
    "                                        #  maxPixels= 1e9,\n",
    "                                         bestEffort= True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c266db1c",
   "metadata": {},
   "source": [
    "## Convert inundated area to hectares\n",
    "We use the `inundated_stats`  variable (type dictionary) from the previous step. \n",
    "\n",
    "__-->__ The function `.getNumber`  may be applied to different types of datasets/variables as can be seen here https://developers.google.com/s/results/earth-engine?hl=en&q=getNumber\n",
    "\n",
    "The function `Dictionary.getNumber`  is used in this case requires one argument, which is the `key`  (value \"VH\").\n",
    "\n",
    "Then, we use the function `.divide`  and the value 10000, which like the name indicates divides the input value (inundation_stats), which is in square meters, by 10000 to get the result in hectares.\n",
    "\n",
    "Finally, we use the function `.round`  to round up the result in hectares.\n",
    "\n",
    "\n",
    "### Print the results to display\n",
    "Here, we simply use the python `print`  command. To obtain the actual value (hectares), we must use the GEE command `.getInfo()`  as seen below. Additionally, the result is an integer that must be converted to a python string for printing using the command ` str()` \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea8237b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert inundated extent to hectares\n",
    "inundation_area_ha = inundation_stats \\\n",
    "                          .getNumber(\"VH\") \\\n",
    "                          .divide(10000) \\\n",
    "                          .round()\n",
    "\n",
    "# print the results to display\n",
    "print(f'Calculation result --> Estimated flood extent from Hurricane Maria:\\n {inundation_area_ha.getInfo()} Hectares')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fafa959",
   "metadata": {},
   "source": [
    "# Analysis of Multispectral imaging from Sentinel-2\n",
    "## Create a function for cloud masking specific to Sentinel-2 images\n",
    "This function is simply converted from Javascript syntax to python. It already exists on GEE documentation here https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2?hl=en#bands\n",
    "\n",
    "In summary, we can use the **pixel quality (QA) band** of a Sentinel-2 image to identify **cloudy pixels** and mask them out. (The QA band varies from one satellite to another and a different function would be needed for each one). For Sentinel-2, we use the band ` QA10`  and the data-bits 10 and 11 to identify clouds and cirrus pixels. We create a mask by using the command `.bitwiseAnd`  and the filter `.eq`  to identify pixels where the data bits 10 and 11 are equal to Zero. Finally, the cloud mask is applied to the input image of the function (which would be a Sentinel-2 image) using the command `.updateMask`  and the result is divided by 10000 to rescale the reflectance values.\n",
    "\n",
    "### Set visualization parameters for S2 imagery\n",
    "This is done by creating a dictionary with keys representing different parameters and the values of each key is assigned by the user. In this example we define the following map visualization parameters (specifically for Sentinel-2 images:\n",
    "     \n",
    "- The selected bands are B4, B3, and B2 (represent Red, Green, and Blue)\n",
    "- The minimum and maximum reflectance values of 0, and 0.4 (two different parameters)\n",
    "- Gamma values of 0.98, 1.1 and 1 for the respective bands\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80dd81ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to mask clouds using the Sentinel-2 QA band\n",
    "def maskS2clouds(image):\n",
    "    qa = image.select('QA60')\n",
    "\n",
    "    # Bits 10 and 11 are clouds and cirrus, respectively.\n",
    "    cloudBitMask = 1 << 10\n",
    "    cirrusBitMask = 1 << 11\n",
    "\n",
    "    # Both flags should be set to zero, indicating clear conditions.\n",
    "    mask = qa.bitwiseAnd(cloudBitMask).eq(0) \\\n",
    "            .And(qa.bitwiseAnd(cirrusBitMask).eq(0))\n",
    "\n",
    "    return image.updateMask(mask).divide(10000)\n",
    "\n",
    "# Set visualization parameters for S2 imagery\n",
    "vizParams = {\n",
    "  'bands': ['B4', 'B3', 'B2'],\n",
    "  'min': 0,\n",
    "  'max': 0.4,\n",
    "  'gamma': [0.98, 1.1, 1]\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "224bddcf",
   "metadata": {},
   "source": [
    "## Load post-hurricane Sentinel-2 Top-Of-Atmosphere (TOA) imagery\n",
    "We do this separately before the event and then after the event for our study area. \n",
    "\n",
    "__-->__ ***The COPERNICUS/S1_GRD*** are ***satellite images from Sentinel-2*** that are part of the GEE collection. More information can be found here: https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2?hl=en#description\n",
    "\n",
    "This dataset contains historical satellite images since 2015 (ongoing satellite mission). It can be loaded using the command `ee.ImageCollection`. Then, we apply multiple filters using the command `.filter`  to select a specific cloud percentage threshold, dates, and boundary. More specifically:\n",
    "\n",
    "The command `ee.Filter.lt` allows us to select images with cloud cover less-than 10%\n",
    "\n",
    "The command `.filterDate` allows us to extract available satellite images between a start anad end data.\n",
    "Here, we used `2017-10-01 to 2018-10-31` time series to minimize data gap caused by cloudy pixels.\n",
    "\n",
    "The command `.filterBounds` allows us to clip the image to our \"roi\" study area\n",
    "\n",
    "The command `.map` allows us to apply the cloud masking function (defined above) to remove cloudy pixels\n",
    "\n",
    "\n",
    "### Calculate the median pixel value from the Sentinel-2 time series \n",
    "This is done using the command `.median()`  applied to the time series variable. Next, we clip the median image using the `.clip`  command and our study area polygon variable ` (roi)`  that was created in a previous step. \n",
    "\n",
    "\n",
    "### Load a new blank map and add our results\n",
    "\n",
    "- First, we set the map center coordinates and zoom level (same step as before).\n",
    "- Add the \"median image\" layer using the visualization parameters defined in the previous step.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18354a20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a variable for post-hurricane Sentinel-2 TOA imagery\n",
    "# and filter by cloud cover (<10%), date, roi, and apply cloud mask\n",
    "post_Maria = ee.ImageCollection('COPERNICUS/S2') \\\n",
    "               .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10)) \\\n",
    "               .filterDate('2017-10-01', '2018-10-31') \\\n",
    "               .filterBounds(roi) \\\n",
    "               .map(maskS2clouds)\n",
    "\n",
    "# Take median pixel value from time series and add as a Layer\n",
    "medianpixels1 = post_Maria.median()\n",
    "medianpixelsclipped1 = medianpixels1.clip(roi)\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(medianpixelsclipped1, vizParams, 'Post-Maria S2')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0436a072",
   "metadata": {},
   "source": [
    "- **Repeat the previous step to load pre-hurricane Sentinel-2 Top-Of-Atmosphere (TOA) imagery** - This time, we use a different dates interval for the time series: '2015-08-31', '2017-08-31'. The only other change is the variable names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15c75893",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a variable for pre-hurricane Sentinel-2 TOA imagery\n",
    "# and filter by cloud cover (<10%), date, roi, and apply cloud mask\n",
    "pre_Maria = ee.ImageCollection('COPERNICUS/S2') \\\n",
    "              .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10)) \\\n",
    "              .filterDate('2015-08-31', '2017-08-31') \\\n",
    "              .filterBounds(roi) \\\n",
    "              .map(maskS2clouds)\n",
    "\n",
    "# Take median pixel value from time series and add as a Layer\n",
    "medianpixels2 = pre_Maria.median()\n",
    "medianpixelsclipped2 = medianpixels2.clip(roi)\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(medianpixelsclipped2, vizParams, 'Pre-Maria S2')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1eb58fb4",
   "metadata": {},
   "source": [
    "## Calculate post-hurricane NDVI image\n",
    "The formula for NDVI (Normalized Difference Vegetation Index) is: \n",
    "\n",
    "$$\\frac{(NIR - Red)}{(NIR + Red)}$$\n",
    "\n",
    "This is done using the command `.normalizedDifference`  and the bands ` B8`  (Near InfraRed or NIR) and ` B4`  (Red) as input. More information on this function can be found here https://developers.google.com/earth-engine/apidocs/ee-image-normalizeddifference?hl=en\n",
    "\n",
    "__--> What is an index (in Remote Sensing)?__ An index is a formula computed with the data from multi-spectral imagery bands. The values of these indices can be used to identify different types of material. NDVI is a common one used for identifying vegetation. For a list of common indices used, see the [index database](https://www.indexdatabase.de/). You can see all of the indexes that can be calculated from Sentinel-2A data in [this page](https://www.indexdatabase.de/db/is.php?sensor_id=96)\n",
    "\n",
    "- **Set visualization parameters for NDVI images** - This is done by creating a dictionary with keys representing different parameters and the values of each key is assigned by the user. In this example we define the following map visualization parameters (specifically for Sentinel-2 images:\n",
    "\n",
    "      - The minimum and maximum reflectance values of -0.2, and 0.8 (two different parameters)\n",
    "      - Select a custom color palette with ['blue', 'white', 'yellow', 'green']\n",
    "\n",
    "- **clip the NDVI image** using the `.clip`  command and our study area polygon variable ` (roi)`  that was created in a previous step. \n",
    "\n",
    "\n",
    "- __Load a new blank map and add our results__\n",
    "\n",
    "      - First, we set the map center coordinates and zoom level (same step as before).\n",
    "      - Add the \"NDVI\" layer using the new visualization parameters defined in this step.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b16e585",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create post-hurricane NDVI\n",
    "post_ndvi = medianpixels1.normalizedDifference(['B8', 'B4'])\n",
    "\n",
    "# Set visualization parameters for NDVI\n",
    "visParams_ndvi = {'min': -0.2, 'max': 0.8, 'palette': ['blue', 'white', 'yellow', 'green']}\n",
    "\n",
    "# Clip post-hurricane NDVI to roi\n",
    "post_ndvi_clip = post_ndvi.clip(roi)\n",
    "\n",
    "# Add post-hurricane NDVI as a Layer\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(post_ndvi_clip, visParams_ndvi, 'Post-Maria NDVI')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29fd6e41",
   "metadata": {},
   "source": [
    "### Repeat the previous step to create pre-hurricane NDVI image\n",
    "This time, we use the pre-hurricane time series median image as input. The only other change is the variable names. The same visualization parameters defined in the previous step are used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "550b9b4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create pre-hurricane NDVI\n",
    "pre_ndvi = medianpixels2.normalizedDifference(['B8', 'B4'])\n",
    "pre_ndvi_clip = pre_ndvi.clip(roi)\n",
    "\n",
    "# Add pre-hurricane NDVI as a Layer\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(pre_ndvi_clip, visParams_ndvi, 'Pre-Maria NDVI')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dee5c1c-bd55-4c53-99c0-a8d93c0d1d16",
   "metadata": {},
   "source": [
    "## Manually computing the index\n",
    "__-->__ The **normalized difference** form $\\frac{x-y}{x+y}$ is common enough that `ee` provides a function for it. However, you can also manually compute an index using the `.add, .subtract, .multiply` and `.divide` functions with the respective bands. This is useful when calculating more complex index formulae. Here we show how to manually compute the **NDVI**, and we can see that we get the same result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2739eac7-6d7e-4fdb-8d6b-e5c98e618fef",
   "metadata": {},
   "outputs": [],
   "source": [
    "manual_post_ndvi = medianpixels1.select('B8')\\\n",
    "                                .subtract(medianpixels1.select('B4'))\\\n",
    "                                .divide(medianpixels1.select('B8')\\\n",
    "                                        .add(medianpixels1.select('B4'))).clip(roi)\n",
    "\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(manual_post_ndvi, visParams_ndvi, 'Post-Maria NDVI - manual computation')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acb8f356",
   "metadata": {},
   "source": [
    "## Calculate NDVI difference\n",
    "We use the command `.subtract`  applied to the pre-hurricane NDVI image and subtract from it the post-hurricane NDVI image. \n",
    "\n",
    "\n",
    "### Load a new blank map and add our results\n",
    "\n",
    "- First, we set the map center coordinates and zoom level (same step as before).\n",
    "- Add the \"NDVI difference\" layer using visualization parameters specifically defined in this step.\n",
    "  * min, max values of -0.2 and 0.2; everything below -0.2 is mapped to the min value, and everything above 0.2 is mapped to the max\n",
    "  * a custom color palette with ['red', 'yellow', 'green']: min value is visualized as red, max is as gree, and linearly interpolated in between\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eaa24a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create NDVI Difference Map\n",
    "NDVI_diff = post_ndvi_clip.subtract(pre_ndvi_clip)\n",
    "\n",
    "# Add NDVI difference as a Layer\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(NDVI_diff, {'min': -0.2, 'max': 0.2, 'palette': ['red', 'yellow', 'green']}, 'NDVI Difference')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a030ecb",
   "metadata": {},
   "source": [
    "## Estimating number of people affected by flooding\n",
    "### Load the JRC Global Human Settlement Popluation Density layer (250m)\n",
    "This layer provides an estimate of the number of people per cell - This is done using the command ` ee.Image`. Then, we clip the image using our study area polygon ` roi`.\n",
    "\n",
    "\n",
    "__-->__ **The JRC Global Human Settlement Popluation Density layer** is part of the GEE collection. More information can be found here: https://developers.google.com/earth-engine/datasets/catalog/JRC_GHSL_P2016_POP_GPW_GLOBE_V1?hl=en\n",
    "\n",
    "\n",
    "### Exclude population outside non-flooded areas\n",
    "Here, we use the command `.updateMask` with our inundation mask to exclude non-flooded pixels.\n",
    "\n",
    "\n",
    "### Calculate the total population affected by inundation\n",
    "Here, we use the function `.reduceRegion`  to compute the sum (area) of all the inundated pixels. The options used in this example are:\n",
    "\n",
    "- `reducer = ee.Reducer.sum()` --> to select the Sum statistics \n",
    "- `geometry = roi` --> to limit the extent to our study area\n",
    "- `scale = 250` --> which corresponds to the population layer spatial resolution in meters\n",
    "- `maxPixels = 1e9` --> is the maximum number of pixels to consider in this GEE request\n",
    "\n",
    "### print the results to display\n",
    "Here, we simply use the python print command. To obtain the actual value (population), we must use the GEE command `.getInfo()` as seen below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebf60c28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import JRC Global Human Settlement Popluation Density layer (250m) - number of people per cell\n",
    "population_count = ee.Image('JRC/GHSL/P2016/POP_GPW_GLOBE_V1/2015').clip(roi)\n",
    "\n",
    "# Create a raster showing exposed population only using the resampled flood layer\n",
    "population_exposed = population_count.updateMask(inundation3).updateMask(population_count)\n",
    "\n",
    "# Sum pixel values of exposed population raster\n",
    "stats = population_exposed.reduceRegion(reducer= ee.Reducer.sum(),\n",
    "                                        geometry= roi,\n",
    "                                        scale= 250,\n",
    "                                        maxPixels=1e9)\n",
    "\n",
    "# Get number of exposed people as integer\n",
    "number_pop_exposed = stats.getNumber('population_count').round()\n",
    "print(f'Calculation result --> Estimated Exposed population from Hurricane Maria:\\n{number_pop_exposed.getInfo()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6316d574",
   "metadata": {},
   "source": [
    "## Estimate effect of flooding on land usage type\n",
    "### Load the Copernicus Global Land Service (CGLS) land cover (100m)\n",
    "This is done using the command `ee.Image`. Then, we select the band `discrete classification` (which provides integer values representing the different land cover categories). Finally, we clip the image using our study area polygon `roi`.\n",
    "\n",
    "\n",
    "__-->__ **The Copernicus Global Land Service (CGLS) land cover (100m) layer** is part of the GEE collection. More information can be found here: https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_Landcover_100m_Proba-V-C3_Global?hl=en\n",
    "\n",
    "\n",
    "### Load a new blank map and add our results\n",
    "- First, we set the map center coordinates and zoom level (same step as before).\n",
    "- Add the \"NDVI difference\" layer using the default visualization parameters \"empty dictionary {}\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fb326a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a variable for Copernicus Global Land Service (CGLS) land cover (100m)\n",
    "# Clip land cover to roi and add as a Layer\n",
    "lc = ee.Image(\"COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019\") \\\n",
    "       .select('discrete_classification').clip(roi)\n",
    "\n",
    "# Add pre-hurricane NDVI as a Layer\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(lc, {}, \"Land Cover\")\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccd25aa9-3978-4f43-a7d4-ebecc8d5ed80",
   "metadata": {},
   "source": [
    "To get a sense of what each classification represents, you can get their `properties` using the `.get` function. Here we construct a dictionary of the class id's and their descriptions from the dataset's [various properties](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_Landcover_100m_Proba-V-C3_Global?hl=en#image-properties)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d863aee-1312-4756-8751-069f06df3b1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "lc_class_dict = {class_id: label for class_id, label in zip(lc.get('discrete_classification_class_values').getInfo(),\n",
    "                                                  lc.get('discrete_classification_class_names').getInfo())}\n",
    "lc_class_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da3bcec5",
   "metadata": {},
   "source": [
    "## Extract cropland from the land cover map\n",
    "First, we generate a binary mask with the land cover value of 40, which corresponds to cropland (Cultivated and managed vegetation / agriculture). This is done using the command `.eq`  to apply the mask \"equal-to\". Then we apply this mask to the land cover layer itself using the command `.updateMask`  similar to previous steps.\n",
    "\n",
    "\n",
    "### Extract the affected cropland \n",
    "Here, we use the command `.updateMask`  applied to the latest inundation variable (slope-corrected, removed isolated pixels) to get only cropland pixels that were affected by inundation.\n",
    "\n",
    "\n",
    "### Calculate the affected cropland extent\n",
    "First, we multiply the affected cropland image by the pixel area value. For this, we use the command `.multiply`  and combine it with the function ` ee.Image.pixelArea`.\n",
    "\n",
    "\n",
    "### Calculate the total area of affected cropland pixels\n",
    "Here, we use the function `.reduceRegion`  to compute the sum (area) of all the inundated cropland pixels. The options used in this example are:\n",
    "\n",
    "- `reducer = ee.Reducer.sum()` --> to select the Sum statistics \n",
    "- `geometry = roi` --> to limit the extent to our study area\n",
    "- `scale = 100` --> which corresponds to the native spatial resolution of the land cover image\n",
    "- `maxPixels = 1e9` --> is the maximum number of pixels in this GEE request\n",
    "\n",
    "\n",
    "### Convert inundated area to hectares\n",
    "We use the `crop_stats`  variable (dictionary) and the function `.getNumber`  to extract the `key`  value \"VH\". Then, we use the function `.divide`  and the value 10000 to convert the result from square meters to hectares. Finally, we use the function `.round`  to round up the result in hectares.\n",
    "\n",
    "\n",
    "- __print the results to display__ - Here, we simply use the python print command. To obtain the actual value (hectares), we must use the GEE command .getInfo() as seen below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3975e392",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract only cropland pixels from CGLS using class value equal to 40\n",
    "# (e.g., cultivated and managed vegetation/agriculture)\n",
    "cropmask = lc.eq(40)\n",
    "cropland = lc.updateMask(cropmask)\n",
    "\n",
    "# Extract the affected cropland using the flood layer\n",
    "cropland_affected = inundation3.updateMask(cropland)\n",
    "\n",
    "# Get pixel area of affected cropland layer\n",
    "crop_pixelarea = cropland_affected.multiply(ee.Image.pixelArea());\n",
    "\n",
    "# Sum pixels of affected cropland layer\n",
    "crop_stats = crop_pixelarea.reduceRegion(\n",
    "                              reducer = ee.Reducer.sum(), #sum all pixels with area information\n",
    "                              geometry = roi,\n",
    "                              scale= 100,\n",
    "                              maxPixels= 1e9)\n",
    "\n",
    "# Convert area to hectares\n",
    "crop_area_ha = crop_stats \\\n",
    "                  .getNumber(\"VH\") \\\n",
    "                  .divide(10000) \\\n",
    "                  .round()\n",
    "\n",
    "# Print results to Console\n",
    "print(f'Calculation result --> Estimated Inundated Cropland area from Hurricane Maria:\\n{crop_area_ha.getInfo()} Hectares')\n",
    "#print (crop_area_ha, 'Hectares of Inundated Cropland')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f5b504",
   "metadata": {},
   "source": [
    "## Visualize before/after\n",
    "### Set map visualization parameters\n",
    "This is done for the cropland layers before and after the hurricane. The selected Min and Max are the same (0 - 14). The only difference is the color (Yellow: before hurricane, and Red: after hurricane).\n",
    "\n",
    "\n",
    "### Load a new blank map and add our results \n",
    "\n",
    "- First, we set the map center coordinates and zoom level (same step as before).\n",
    "- Then, we load the study area polygon in Gray color to improve contrast and better see the results.\n",
    "- Add the two \"cropland\" layers (before and after hurricane) using their respective visualization parameters.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32e76f6b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Set cropland visualization parameters\n",
    "croplandVis_before = {'min': 0, 'max': 14.0, 'palette': ['Yellow'],}\n",
    "croplandVis_after = {'min': 0, 'max': 14.0, 'palette': ['Red'],}\n",
    "\n",
    "# Add cropland to map\n",
    "Map = geemap.Map(center=[18.2208, -66.5901], zoom = 9)\n",
    "Map.addLayer(roi, {'color': 'gray'}, 'Study Area')\n",
    "Map.addLayer(cropland, croplandVis_before, 'Cropland')\n",
    "\n",
    "# Add affected cropland to map\n",
    "Map.addLayer(cropland_affected, croplandVis_after, 'Affected Cropland')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28be9e9b-a037-4ca1-8872-61969086aae9",
   "metadata": {},
   "source": [
    "## Save a file from GEE\n",
    "If you want to have a copy of the file, you can save it to a `.tiff` to open in another program, such as QGIS or in python with `rasterio`.\n",
    "\n",
    "We can use the `geemap.ee_export_image` [function](https://geemap.org/notebooks/11_export_image/#download-an-eeimage) to do so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5823b140-3490-49f4-9a7d-544eee986f6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "geemap.ee_export_image(cropland_affected, filename='affected_cropland.tif', scale=100, region=roi.geometry())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6f838ae",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Challenge\n",
    "\n",
    "Using this script as a template, find another major event and examine its impacts. \n",
    "\n",
    "- Remember that the DEM, the population layer, and the land cover layer all have **global coverage**. \n",
    "\n",
    "- Remember that **Sentinel-1** SAR data is only available **since 2014** and that **Sentinel-2** data us only available **since 2015**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "239f6fb5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}