flood-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "<img src=\"images/Tko46mp.png\" width=\"1000\" style=\"display:block; margin:auto\" />\n",
    "\n",
    "<br>\n",
    "\n",
    "In mid-June 2022, Pakistan began experiencing unusually heavy rainfall. Within weeks, the relentless downpours had led to widespread flooding across significant parts of the country. Pakistan's climate is characterized by a distinct monsoon season, with the majority of its annual rainfall, approximately 70-75%, occurring between June and September. However, the intensity of these monsoon rains varies annually, with some years witnessing far higher rainfall than usual.\n",
    "\n",
    "In 2022, the monsoon season was exceptionally severe. Described by U.N. Secretary-General António Guterres as a \"monsoon on steroids,” Pakistan received about ten times more rain than its average. Some regions were particularly hard-hit, receiving 600-800% percent of their typical August rainfall. At the peak of this crisis, Sindh, the country's most populous province, saw an overwhelming 15 inches of rain in just one day. As the summer progressed, these waters accumulated in the low-lying areas surrounding the Indus River, transforming vast stretches of farmland into virtual lakes and isolating numerous villages.\n",
    "\n",
    "The catastrophic floods resulted in a tragic loss of life and widespread devastation. A total of 1,739 individuals lost their lives, including 647 children, and there were 12,867 reported injuries. The disaster left over 2.1 million people without homes, marking it as the most severe flood in Pakistan since 2010, which had a death toll of nearly 2,000. These floods were also the most catastrophic globally since the South Asian floods in 2020. \n",
    "\n",
    "Synthetic Aperture Radar (SAR) offers several significant advantages for flood mapping, making it a highly effective tool in disaster management and environmental monitoring. Unlike optical sensors, SAR can penetrate through clouds and is not affected by weather conditions. This feature is particularly crucial for flood mapping, as heavy cloud cover often accompanies flooding events. Additionally, SAR systems can operate effectively both day and night, ensuring continuous monitoring of flood situations. Lastly, SAR is particularly sensitive to smooth surfaces like water, which appear dark on SAR images. This contrast makes it easier to differentiate between flooded and non-flooded areas, enhancing the accuracy of flood extent mapping.\n",
    "\n",
    "The following tutorial is based upon a [UN-SPIDER Google Earth Engine (GEE) JavaScript Script](https://www.un-spider.org/advisory-support/recommended-practices/recommended-practice-google-earth-engine-flood-mapping/step-by-step) for flood mapping. It has been translated to Python, heavily adapted, and expanded. \n",
    "\n",
    "Let's get started... \n",
    "\n",
    "<hr>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Import required packages\n",
    "\n",
    "Let's start by importing the necessary libraries. \n",
    "\n",
    "1. **Earth Engine Python API (`ee`)** enables interaction with the Google Earth Engine servers, allowing for complex geospatial computations and access to a vast catalogue of satellite imagery and geospatial datasets. It allows Google Earth Engine to be used on your local machine without the need to utilise the web-based JavaScript IDE. \n",
    "\n",
    "2. **Geemap** facilitates interactive mapping with Google Earth Engine. It's particularly user-friendly for those familiar with Jupyter notebooks and is built on ipyleaflet and folium libraries, providing robust tools for Earth Engine data visualisation and analysis. More info is available here: https://geemap.org/ \n",
    "\n",
    "3. **Geopandas** extends the datatypes used by pandas to allow spatial operations on geometric types, making it a powerful library for geospatial analysis.\n",
    "\n",
    "**Note:** Before running these import statements, you should ensure that these libraries are installed in your Python environment. You can install them using pip or conda if they are not already installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pip install earthengine-api geemap geopandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import required libraries\n",
    "import ee\n",
    "import geemap\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Authenticate/initiate Earth Engine Python API\n",
    "\n",
    "For this tutorial, you'll need an Earth Engine account. If you don't already have one, you can register at [Google Earth Engine Registration](https://code.earthengine.google.com/register)\n",
    "\n",
    "The initial execution of the code will require you to authenticate access to the Earth Engine API. The steps are: \n",
    "\n",
    "1. **Run the Authentication Code:** When you run the authentication code segment for the first time, a new browser window will pop up.\n",
    "\n",
    "2. **Follow the Authentication Steps:** In the browser, you'll be guided through a series of steps to authenticate your Earth Engine account. This involves signing in with your Google account credentials and granting necessary permissions.\n",
    "\n",
    "3. **Copy the Authentication Code:** On successful authentication, you'll receive an authentication code. Copy this code.\n",
    "\n",
    "4. **Paste the Code in the notebook:** Return to this notebook. A prompt or input box will appear requesting the authentication code. Paste the code you copied into this box and submit.\n",
    "\n",
    "Once these steps are complete, the Earth Engine API will be authenticated, and you'll be able to run the code segments in this tutorial. This is a one-time process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Authenticate/initiate Earth Engine Python API\n",
    "try:\n",
    "    # Attempt to initialize the Earth Engine API.\n",
    "    ee.Initialize()\n",
    "    print(\"Earth Engine API successfully initialized.\")\n",
    "except:\n",
    "    # If initialization fails, authenticate and then initialize.\n",
    "    print(\"You are not authenticated. Running authentication now...\")\n",
    "    ee.Authenticate()\n",
    "    ee.Initialize()\n",
    "    print(\"Authentication successful. Earth Engine API successfully initialized.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define AOI and load ADM2 (division) boundaries\n",
    "\n",
    "Here we define our Area of Interest (AOI) using a bounding box. This AOI was chosen based on reports of which regions were most affected by the 2022 foods. Our AOI is focused on Sindh, the country's most populous province, which saw overwhelming rainfall. \n",
    "\n",
    "The AOI:\n",
    "\n",
    "<img src=\"images/Np3zgQl.png\" width=\"700\" />\n",
    "\n",
    "Administrative boundary data for ADM2 level (divisions) are imported from a Google Earth Engine (GEE) asset and clipped to our AOI. The purpose of this data is to facilitate the extraction of division-level flood statistics in subsequent analysis steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define AOI bounding box\n",
    "aoi = ee.Geometry.Polygon([\n",
    "    [66.60903973281872,23.261545926289028],\n",
    "    [70.48721356094372,23.261545926289028],\n",
    "    [70.48721356094372,29.286410056079323],\n",
    "    [66.60903973281872,29.286410056079323]])\n",
    "\n",
    "# Import Pakistan adm2 boundaries from GEE assets (to be used for flood statistics)\n",
    "adm2 = (ee.FeatureCollection('projects/ee-swcoughlan/assets/pak_adm2')\n",
    "          .filterBounds(aoi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Define date ranges for before and after the flood\n",
    "\n",
    "The objective here is to capture comparative snapshots before and after the flood to evaluate the extent of inundation. For this, I've selected month-long intervals as reference periods for the `before` and `after` scenarios. Typically, March through May are dry months in Pakistan, preceding the onset of the monsoon season, which extends from early June to September. The 2022 floods, notably starting around June 14th, persisted until September's end.\n",
    "\n",
    "Ensure the ranges are long enough to acquire a sufficient number of tiles for your AOI. This will vary depending on AOI size and location, and you might have to play with the dates. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Before the flood\n",
    "before_start= '2022-03-15'\n",
    "before_end='2022-04-15'\n",
    "\n",
    "# After the flood\n",
    "after_start='2022-08-15'\n",
    "after_end='2022-09-15'\n",
    "\n",
    "# Print date range\n",
    "print(\"Flooding between:\" + \" \" + after_start + \" and \" + after_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Define water threshold and Sentinel-1 parameters\n",
    "\n",
    "Set up some key parameters for the analysis:\n",
    "\n",
    "`water_threshold` is used as a threshold for identifying possible surface water in the SAR imagery. Later, we'll calculate a difference image between 'before' and 'after' and any pixels in that image with a value above this threshold are likely to be classified as water. The value choosen is a general value which works relatively well in dry & wet season. \n",
    "\n",
    "`polarization` refers to the polarization of the radar signal used in the Sentinel-1 imagery. Radar signals can be polarized vertically (V) or horizontally (H). 'VH' means that the signal is transmitted vertically and received horizontally. Choice of polarization can affect the ability to detect certain features. 'VH' polarization is often used in flood detection because it's sensitive to changes in surface roughness and can help distinguish water from land.\n",
    "\n",
    "`orbital_pass` refers to the direction of Sentinel-1's orbit at the time the image was captured. Satellites can capture images on the ascending (from south to north) or descending (from north to south) part of their orbit. Choice of orbital pass can affect the viewing geometry, so for this analysis it's important to only use one orbital direction. You can change this parameter if your image collection doesn't have many images. Some regions have more ascending images than descending and vice versa. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the threshold for identifying possible surface water.\n",
    "# This is a general value which works relatively well in dry & wet season.\n",
    "water_threshold = 1.19\n",
    "\n",
    "# Specify polarization(VH or VV) \n",
    "polarization = 'VH'\n",
    "\n",
    "# Specify orbital pass (ASCENDING or DESCENDING)\n",
    "orbital_pass = 'ASCENDING'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Required functions\n",
    "\n",
    "Let's create some necessary functions to be used in subsequent steps. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6.1. Cloud mask\n",
    "\n",
    "This function is designed to eliminate or 'mask' cloud cover from Sentinel-2 images. We'll generate a composite image from Sentinel-2 data, collected over a span of one month, to facilitate a visual comparison between our SAR flood data and the optical imagery obtained from Sentinel-2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to mask clouds using the Sentinel-2 QA band\n",
    "def mask_s2_clouds(image):\n",
    "    qa = image.select('QA60')\n",
    "\n",
    "    # Bits 10 and 11 are clouds and cirrus, respectively.\n",
    "    cloud_bit_mask = 1 << 10\n",
    "    cirrus_bit_mask = 1 << 11\n",
    "\n",
    "    # Both flags should be set to zero, indicating clear conditions.\n",
    "    mask = qa.bitwiseAnd(cloud_bit_mask).eq(0) \\\n",
    "        .And(qa.bitwiseAnd(cirrus_bit_mask).eq(0))\n",
    "\n",
    "    return image.updateMask(mask).divide(1e4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6.2. Normalized Difference Water Index (NDWI) \n",
    "\n",
    "This function adds a Normalized Difference Vegetation Index (NDWI) band to the Sentinel-2 imagery. NDWI is primarily used to monitor vegetation and water content, and is particularly effective at quantifing the presence and extent of water bodies. NDWI is calculated by comparing the reflectance values of near-infrared (NIR) and shortwave infrared (SWIR) bands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to add Normalized Difference Vegetation Index (NDWI) to Sentinel-2 imagery\n",
    "def add_ndwi(image):\n",
    "\n",
    "    # Calculating NDWI using the normalizedDifference function\n",
    "    ndwi = image.normalizedDifference(['B3', 'B8']).rename('NDWI')\n",
    "\n",
    "    # Add band to the image\n",
    "    return image.addBands([ndwi])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6.3. Water mask\n",
    "\n",
    "This function creates an estimated surface water mask\n",
    "\n",
    "A water mask is created from the before/after difference image, then gradually refined using ancilliary data sets. A binary mask is initially created where pixels with a value greater than the threshold are considered water. Next, the flood mask is refined using the Global Surface Water dataset from the Joint Research Centre (JRC). Perennial water bodies (present more than 10 months of the year) are removed from the flood analysis. The flood mask is further refined by removing pixels connected to 8 or fewer neighbors. This helps to remove isolated pixels that may be noise. Areas with more than 5 percent slope and more than 100m elevation are excluded using a Digital Elevation Model (DEM) from the World Wildlife Fund (WWF) HydroSHEDS dataset. This is based on the assumption that areas with high slope or elevation are less likely to be flooded. Lastly, the flood extent area is calculated by multiplying the flooded mask with the pixel area. The areas of flooded pixels within the AOI are summed and converted to square kilometers and hectares. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to create estimated surface water mask\n",
    "def create_water_mask(img, threshold, scale_value):\n",
    "    # Create water mask image using VV/VH difference (estimated surface water) threshold\n",
    "    water_mask = img.select('VH').gt(threshold)\n",
    "\n",
    "    # Refine flood result using additional datasets\n",
    "    swater = ee.Image('JRC/GSW1_0/GlobalSurfaceWater').select('seasonality')\n",
    "    swater_mask = swater.gte(10).updateMask(swater.gte(10))\n",
    "    \n",
    "    # Flooded layer where perennial water bodies (water > 10 mo/yr) are assigned a 0 value\n",
    "    flooded_mask = water_mask.where(swater_mask, 0)\n",
    "\n",
    "    # Final flooded area without pixels in perennial waterbodies\n",
    "    flooded = flooded_mask.updateMask(flooded_mask)\n",
    "\n",
    "    # Compute connectivity of pixels to eliminate those connected to 8 or fewer neighbors\n",
    "    connections = flooded.connectedPixelCount()    \n",
    "    flooded = flooded.updateMask(connections.gte(8))\n",
    "\n",
    "    # Mask out areas with > 5% slope and > 100m elevation using a Digital Elevation Model \n",
    "    dem = ee.Image('WWF/HydroSHEDS/03VFDEM')\n",
    "    elevation = dem.select('b1')\n",
    "    terrain = ee.Algorithms.Terrain(dem)\n",
    "    slope = terrain.select('slope')\n",
    "    flooded = flooded.updateMask(slope.lt(5))\n",
    "    flooded = flooded.updateMask(elevation.lt(100))\n",
    "\n",
    "    # Calculate flood extent area\n",
    "    flood_pixelarea = flooded.select('VH').multiply(ee.Image.pixelArea())\n",
    "\n",
    "    # Sum the areas of flooded pixels\n",
    "    flood_stats = flood_pixelarea.reduceRegion(\n",
    "        reducer=ee.Reducer.sum(),\n",
    "        geometry=aoi,\n",
    "        scale=10,  # native Sentinel-1 resolution \n",
    "        bestEffort=True\n",
    "    )\n",
    "\n",
    "    # Convert the flood extent to square km and hectares\n",
    "    flood_area_km2 = flood_stats.getNumber('VH').divide(1e6).round()\n",
    "    flood_area_ha = flood_stats.getNumber('VH').divide(1e4).round()\n",
    "    \n",
    "    return flooded, flood_area_km2, flood_area_ha\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6.4. Split AOI\n",
    "\n",
    "This function segments our large AOI into smaller, more manageable tiles. This is intended to avoid processing the whole AOI in one go, which would be computationally intensive and lead to errors during processing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to split AOI into smaller tiles.\n",
    "def split_aoi_into_tiles(aoi, tile_size):\n",
    "    \n",
    "    tiles = []\n",
    "\n",
    "    # Get the bounds of the AOI and extract the coordinates\n",
    "    bounds = aoi.bounds()\n",
    "    coords = bounds.coordinates().getInfo()[0]\n",
    "\n",
    "    # Extract the min and max coordinates\n",
    "    min_lon = min(coord[0] for coord in coords)\n",
    "    max_lon = max(coord[0] for coord in coords)\n",
    "    min_lat = min(coord[1] for coord in coords)\n",
    "    max_lat = max(coord[1] for coord in coords)\n",
    "\n",
    "    current_lat = min_lat\n",
    "    while current_lat < max_lat:\n",
    "        current_lon = min_lon\n",
    "        while current_lon < max_lon:\n",
    "            upper_left = [current_lon, current_lat]\n",
    "            lower_right = [min(current_lon + tile_size, max_lon),\n",
    "                           min(current_lat + tile_size, max_lat)]\n",
    "\n",
    "            # Create a rectangle tile\n",
    "            tile_coords = [\n",
    "                upper_left,\n",
    "                [lower_right[0], upper_left[1]],\n",
    "                lower_right,\n",
    "                [upper_left[0], lower_right[1]],\n",
    "                upper_left  # Closing the loop\n",
    "            ]\n",
    "\n",
    "            tile = ee.Geometry.Polygon([tile_coords])\n",
    "            tiles.append(tile)\n",
    "            current_lon += tile_size\n",
    "\n",
    "        current_lat += tile_size\n",
    "\n",
    "    return tiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Load Sentinel-2 \n",
    "\n",
    "Load the Sentinel-2 imagery for our AOI over a one-month period. Images are filtered to ensure that cloud cover is below 50%. Cloud masking and NDWI functions are applied. \n",
    "\n",
    "A composite image is created by averaging the collection and clipping it to the exact dimensions of the AOI. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Sentinel-2 image collection\n",
    "s2 = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')\n",
    "        .filterDate('2022-08-15', '2022-09-15')\n",
    "        .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 50))\n",
    "        .map(mask_s2_clouds)\n",
    "        .map(add_ndwi)\n",
    "        .filterBounds(aoi)\n",
    "        )\n",
    "\n",
    "# Create a composite and clip to AOI\n",
    "s2_aoi = s2.mean().clip(aoi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Load Sentinel-1 collection, create before & after subsets, apply speckle filter, and calculate the before/after difference. \n",
    "\n",
    "Load Sentinel-1 images for our AOI. Images are filtered based on instrument mode, polarization, orbital pass, resolution, and geographic bounds. \n",
    "\n",
    "We create two subsets of images, one for the period before and one for after the flood event. Both subsets are processed into mosaics, clipped to the AOI, and smoothed using a speckle filter. \n",
    "\n",
    "Finally, we calculate the change between the 'before' and 'after' periods by dividing the smoothed images to create our difference image (to which our flood mask function will be applied)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Sentinel-1 image collection\n",
    "s1 = (ee.ImageCollection('COPERNICUS/S1_GRD')\n",
    "        .filter(ee.Filter.eq('instrumentMode', 'IW'))\n",
    "        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', polarization))\n",
    "        .filter(ee.Filter.eq('orbitProperties_pass', orbital_pass))\n",
    "        .filter(ee.Filter.eq('resolution_meters', 10))\n",
    "        .filterBounds(aoi)\n",
    "        .map(terrain_correction) # Apply terrain correction\n",
    "        )\n",
    "\n",
    "# Filter by date\n",
    "before_collection = s1.select('VH').filterDate(before_start, before_end)\n",
    "after_collection = s1.select('VH').filterDate(after_start, after_end)\n",
    "\n",
    "# Function to extract date from metadata\n",
    "def dates(imgcol):\n",
    "    range = imgcol.reduceColumns(ee.Reducer.minMax(), [\"system:time_start\"])\n",
    "    printed = (ee.String('from ')\n",
    "               .cat(ee.Date(range.get('min')).format('YYYY-MM-dd'))\n",
    "               .cat(' to ')\n",
    "               .cat(ee.Date(range.get('max')).format('YYYY-MM-dd')))\n",
    "    return printed.getInfo()\n",
    "\n",
    "# Print dates of before images\n",
    "before_count = before_collection.size().getInfo()\n",
    "print(f\"Tiles selected: Before Flood ({before_count})\", dates(before_collection))\n",
    "\n",
    "# Print dates of after images\n",
    "after_count = after_collection.size().getInfo()\n",
    "print(f\"Tiles selected: After Flood ({after_count})\", dates(after_collection))\n",
    "\n",
    "# Create a mosaic of selected tiles\n",
    "before = before_collection.mosaic().clip(aoi)\n",
    "after = after_collection.mosaic().clip(aoi)\n",
    "\n",
    "# Apply speckle filter\n",
    "smoothing_radius = 50\n",
    "before_smooth = before.focal_mean(smoothing_radius, 'circle', 'meters')\n",
    "after_smooth = after.focal_mean(smoothing_radius, 'circle', 'meters')\n",
    "\n",
    "# Calculate difference between before and after\n",
    "difference = after_smooth.divide(before_smooth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9. Create water mask and calculate estimated flooded area\n",
    "\n",
    "Here we call the `create_water_mask` function to identify potentially flooded areas, using the difference image created in the last step. The function returns the flood mask and the estimated flood area in both square kilometers and hectares. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call the function to create the water mask\n",
    "flooded_area = create_water_mask(difference, water_threshold, 10)\n",
    "\n",
    "# Get the flooded area in km2 and hectares\n",
    "flooded_km2 = flooded_area[1].getInfo()\n",
    "flooded_ha = flooded_area[2].getInfo()\n",
    "\n",
    "flooded = flooded_area[0]\n",
    "\n",
    "# Print the results\n",
    "print('Total potential flooded area (km2):', flooded_km2)\n",
    "print('Total potential flooded area (ha):', flooded_ha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10. Split the AOI into smaller tiles for processing\n",
    "\n",
    "Defines a tile size of 10x10km and applies our `split_aoi_into_tiles` function to the AOI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define tile size in meters\n",
    "tile_size = 10000 # 10km x 10km\n",
    "\n",
    "# Implements function to split AOI into smaller tiles\n",
    "tiles = split_aoi_into_tiles(aoi, tile_size) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 11. Calculate the estimated number of affected people\n",
    "\n",
    "It would be useful to have an estimate of the number of people affected by flooding within our AOI. We can use the Global Human Settlement Layer (GHSL) Population Density layer for 2020 for population data. We must first reproject our flood layer to the same projection as the GHSL layer, in order to maintain spatial accuracy. \n",
    "\n",
    "Using this reprojected flood layer, we create a new layer indicating the population exposed to flooding by masking the population layer with the flood layer. We then call the `pop_within_tile` function to calculate the sum of the exposed population within each subdivided tile of the AOI, aggregating the values to determine the total exposed population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load JRC GHS-POP R2023A Population Density layer 100m (2020)\n",
    "# Number of people per cell is returned and aggregated to the AOI\n",
    "population_count = (ee.Image('JRC/GHSL/P2023A/GHS_POP/2020')\n",
    "                      .clip(aoi)\n",
    "                      )\n",
    "\n",
    "# Get GHSL projection\n",
    "ghs_projection = population_count.projection()\n",
    "\n",
    "# Reproject flood layer to GHSL scale\n",
    "flooded_rpj = flooded.reproject(crs = ghs_projection)\n",
    "\n",
    "# Create a raster showing exposed population only using the resampled flood layer\n",
    "population_exposed = population_count.updateMask(flooded_rpj).updateMask(population_count)\n",
    "\n",
    "# Function to calculate the sum of exposed population within each tile\n",
    "def pop_within_tile(tile):\n",
    "    stats = population_exposed.reduceRegion(\n",
    "        reducer = ee.Reducer.sum(),\n",
    "        geometry = tile,\n",
    "        scale = 100,\n",
    "        maxPixels = 1e9,\n",
    "        bestEffort = True\n",
    "    )\n",
    "    return stats.getNumber('population_count')\n",
    "\n",
    "# Apply sum function to each tile\n",
    "exposed_sums = [pop_within_tile(tile) for tile in tiles] \n",
    "\n",
    "# Convert each ee.Number to Python number\n",
    "exposed_sums_info = [num.getInfo() for num in exposed_sums] \n",
    "\n",
    "# Combine results from all tiles\n",
    "total_exposed_population = round(sum(exposed_sums_info))\n",
    "\n",
    "# Print result\n",
    "print('Estimated number of people exposed:', total_exposed_population)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 12. Calculate the estimated affected cropland\n",
    "\n",
    "To assess the impact of flooding on agricultural areas within our AOI, we can utilise the Copernicus Global Land Cover Layer for 2019. Cropland areas are isolated and a mask that highlights agricultural zones is created. The flood layer is again reprojected to align with the land cover projection for spatial accuracy. Using this reprojected flood layer, a new raster is generated to show the cropland areas potentially affected by flooding. The total area of this exposed cropland is calculated in both square kilometers and hectares by summing the pixels within the flood-affected zones. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Copernicus Global Land Cover Layers: CGLS-LC100 Collection 3 100m (2015-2019)\n",
    "# Select 'discrete classification' band\n",
    "landcover = (ee.Image('COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019')\n",
    "               .select(\"discrete_classification\")\n",
    "               .clip(aoi)\n",
    "               )\n",
    "\n",
    "# Extract cropland pixels using the class 'cultivated and managed vegetation / agriculture' (value = 40)\n",
    "crop_mask = landcover.eq(40)\n",
    "cropland = landcover.updateMask(crop_mask)\n",
    "\n",
    "# Get landcover projection\n",
    "lc_projection = landcover.projection()\n",
    "\n",
    "# Reproject flood layer to GHSL scale\n",
    "flooded_rpj = flooded.reproject(crs = lc_projection)\n",
    "\n",
    "# Create raster of exposed cropland using the resampled flood layer as a mask\n",
    "cropland_exposed = flooded_rpj.updateMask(cropland)\n",
    "\n",
    "# Get pixel area of exposed cropland and urban layers\n",
    "cropland_pixelarea = cropland_exposed.multiply(ee.Image.pixelArea())\n",
    "\n",
    "# Sum pixels of affected cropland layer\n",
    "crop_stats = cropland_pixelarea.reduceRegion(\n",
    "    reducer=ee.Reducer.sum(),\n",
    "    geometry=aoi,\n",
    "    scale=100,\n",
    "    maxPixels=1e9\n",
    ")\n",
    "\n",
    "# Convert area to km2 and hectares\n",
    "crop_area_km2 = crop_stats.getNumber('VH').divide(1000000).round()\n",
    "crop_area_ha = crop_stats.getNumber('VH').divide(10000).round()\n",
    "\n",
    "# Assign to variables\n",
    "crop_km2 = crop_area_km2.getInfo()\n",
    "crop_ha = crop_area_ha.getInfo()\n",
    "\n",
    "# Print result\n",
    "print('Estimated cropland exposed (km2):', crop_km2)\n",
    "print('Estimated cropland exposed (ha):', crop_ha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13. Calculate the estimated affected urban area\n",
    "\n",
    "To assess the impact of flooding on urban areas, we can make use of the Copernicus Global Land Cover Layer for 2019. We isolate urban areas by applying a mask for the 'urban/built-up' class (value 50). The flood layer is reprojected to align with the urban landcover layer projection. A new raster that highlights urban areas affected by the flood is created. The sum of the exposed urban areas within each tile of the AOI is calculated, combining these results to estimate the total urban area affected. This total is then converted to square kilometers and hectares. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Copernicus Global Land Cover Layers: CGLS-LC100 Collection 3 100m (2015-2019)\n",
    "# Select 'discrete classification' band\n",
    "landcover_urban = (ee.Image('COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019')\n",
    "               .select(\"discrete_classification\")\n",
    "               .clip(aoi)\n",
    "               )\n",
    "\n",
    "\n",
    "# Extract urban pixels using the class 'urban / built up' (value = 50)\n",
    "urban_mask = landcover_urban.eq(50)\n",
    "urban = landcover_urban.updateMask(urban_mask)\n",
    "\n",
    "# Get landcover projection\n",
    "urban_projection = landcover_urban.projection()\n",
    "\n",
    "# Reproject flood layer to GHSL scale\n",
    "flooded_rpj = flooded.reproject(crs = urban_projection)\n",
    "\n",
    "# Create raster of exposed urban areas using the resampled flood layer as a mask\n",
    "urban_exposed = flooded_rpj.updateMask(urban)\n",
    "\n",
    "# Get pixel area of exposed urban layer\n",
    "urban_pixelarea = urban_exposed.multiply(ee.Image.pixelArea())\n",
    "\n",
    "# Function to calculate the sum of exposed urban area within each tile\n",
    "def urban_within_tile(tile):\n",
    "    stats = urban_pixelarea.reduceRegion(\n",
    "        reducer = ee.Reducer.sum(),\n",
    "        geometry = tile,\n",
    "        scale = 100,\n",
    "        maxPixels = 1e9,\n",
    "        bestEffort = True\n",
    "    )\n",
    "    return stats.getNumber('VH')\n",
    "\n",
    "# Apply sum function to each tile\n",
    "exposed_sums = [urban_within_tile(tile) for tile in tiles] \n",
    "\n",
    "# Convert each ee.Number to Python number\n",
    "exposed_sums_info = [num.getInfo() for num in exposed_sums] \n",
    "\n",
    "# Combine results from all tiles\n",
    "total_exposed_urban = round(sum(exposed_sums_info))\n",
    "\n",
    "# Convert area to km2 and hectares\n",
    "total_exposed_urban_km2 = round(total_exposed_urban / 1e6)\n",
    "total_exposed_urban_ha = round(total_exposed_urban / 1e4)\n",
    "\n",
    "# Print result\n",
    "print('Estimated urban area exposed (km2):', total_exposed_urban_km2)\n",
    "print('Estimated urban area exposed (ha):', total_exposed_urban_ha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 14. Visualise the results using geemap\n",
    "\n",
    "Now we can sets up and display a multi-layered map using geemap, providing a comprehensive visual analysis of the flood. We define distinct visualization parameters for the various layers. The following layers are added to the map: Sentinel-2 RGB and NDWI, Sentinel-1 (before and after the flood), the difference between these images, a composite of before and after, and layers representing flooded areas, exposed population, cropland, and urban areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define visualisation parameters\n",
    "s1_vis = {\n",
    "    'min': -25, \n",
    "    'max': 0, \n",
    "    'gamma': 2.85\n",
    "    }\n",
    "\n",
    "s2_rgb = {\n",
    "    'min': 0.0,\n",
    "    'max': 0.3,\n",
    "    'bands': ['B4', 'B3', 'B2'],\n",
    "    }\n",
    "\n",
    "s2_ndwi = {\n",
    "    'min': -1.0, \n",
    "    'max': 1.0, \n",
    "    'palette': ['#008000','#FFFFFF','#0000CC'], \n",
    "    'bands': ['NDWI']\n",
    "    }\n",
    "\n",
    "# Define opacity for flood layer\n",
    "opacity = 0.6 \n",
    "\n",
    "# Create a geemap map object\n",
    "Map = geemap.Map(basemap = 'CartoDB.Positron', height='1200px')\n",
    "\n",
    "# Display layers (layers can be loaded activated or deactivated using '0' or '1')\n",
    "Map.centerObject(aoi,8)\n",
    "Map.addLayer(s2_aoi, s2_rgb, 'S2 RGB', 0)\n",
    "Map.addLayer(s2_aoi, s2_ndwi, 'S2 NDWI', 1)\n",
    "Map.addLayer(before_smooth, s1_vis, 'Before Flood', 0)\n",
    "Map.addLayer(after_smooth, s1_vis, 'After Flood', 1)\n",
    "Map.addLayer(difference, {'min': 0, 'max': 2}, 'B/A Difference', 0)\n",
    "Map.addLayer(before_smooth.addBands(after_smooth).addBands(before_smooth), {'min': -25, 'max': 0, 'gamma': 2.0}, 'B/A Composite', 0)\n",
    "Map.addLayer(flooded, {'palette':\"0000FF\"}, 'Flooded Areas', 1, opacity)\n",
    "Map.addLayer(population_exposed, {'palette':\"FF5F15\"}, 'Exposed Population', 0)\n",
    "Map.addLayer(cropland_exposed, {'palette':\"f096ff\"}, 'Exposed Cropland', 0)\n",
    "Map.addLayer(urban_exposed, {'palette':\"fa0000\"}, 'Exposed Urban Area', 0)\n",
    "\n",
    "# Show map\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 15. Print out flood statistics\n",
    "\n",
    "Does exactly what it says on the tin... prints out all the flood statistics from the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print flood statistics\n",
    "print('Flood Statistics:')\n",
    "print('Estimated total flooded area (km2):', flooded_km2)\n",
    "print('Estimated total flooded area (ha):', flooded_ha)\n",
    "print('Estimated number of people exposed:', total_exposed_population)\n",
    "print('Estimated cropland exposed (km2):', crop_km2)\n",
    "print('Estimated cropland exposed (ha):', crop_ha)\n",
    "print('Estimated urban area exposed (km2):', total_exposed_urban_km2)\n",
    "print('Estimated urban area exposed (ha):', total_exposed_urban_ha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 16. Calculate statistics per administrative division\n",
    "\n",
    "We've successfully calculated flood stats for our AOI, but it would be nice to have statistics for each ADM2 division within the AOI. To do this we make use of the GEE asset we loaded in step 3. The function `calculate_flood_extent` computes the flood extent within each administrative division."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 16.1. Calculate estimated flooding extent per administrative division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the ADM2 administrative boundaries\n",
    "divisions = adm2\n",
    "\n",
    "# Function to calculate flooding extent within each division\n",
    "def calculate_flood_extent(feature):\n",
    "\n",
    "    # Calculate flood extent area\n",
    "    flood_pixelarea = flooded.select('VH').multiply(ee.Image.pixelArea())\n",
    "\n",
    "    # Clip the flooded layer to the current administrative boundary\n",
    "    flooded_clipped = flood_pixelarea.clip(feature.geometry())\n",
    "\n",
    "    # Calculate the sum of flooded pixels within the boundary\n",
    "    flood_stats = flooded_clipped.reduceRegion(\n",
    "        reducer=ee.Reducer.sum(),\n",
    "        geometry=feature.geometry(),\n",
    "        scale=10,\n",
    "        maxPixels=1e9,\n",
    "        bestEffort=True\n",
    "    )\n",
    "    # Return the feature with the added property of flood extent\n",
    "    return feature.set(flood_stats)\n",
    "\n",
    "# Apply the function to each feature in the administrative boundaries collection\n",
    "flood_extent_by_division = divisions.map(calculate_flood_extent)\n",
    "\n",
    "# Fetch results\n",
    "flood_result = flood_extent_by_division.getInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 16.2. Calculate estimated exposed population per administrative division\n",
    "\n",
    "The same is done for exposed population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to calculate the sum of exposed population within each division\n",
    "def calculate_pop_exposed(feature):\n",
    "\n",
    "    # Clip the population layer to the current administrative boundary\n",
    "    pop_clipped = population_exposed.clip(feature.geometry())\n",
    "\n",
    "    pop_stats = pop_clipped.reduceRegion(\n",
    "        reducer = ee.Reducer.sum(),\n",
    "        geometry = feature.geometry(),\n",
    "        scale = 100,\n",
    "        maxPixels = 1e9,\n",
    "        bestEffort = True\n",
    "    )\n",
    "    return feature.set(pop_stats)\n",
    "\n",
    "# Apply the function to each feature in the administrative boundaries collection\n",
    "pop_exposed_by_division = divisions.map(calculate_pop_exposed)\n",
    "\n",
    "# Fetch results\n",
    "pop_result = pop_exposed_by_division.getInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 16.3. Calculate estimated exposed cropland per administrative division\n",
    "\n",
    "The same is done for exposed cropland."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to calculate the sum of exposed cropland within each division\n",
    "def calculate_cropland_exposed(feature):\n",
    "    crop_stats = cropland_pixelarea.reduceRegion(\n",
    "    reducer=ee.Reducer.sum(),\n",
    "    geometry=feature.geometry(),\n",
    "    scale=100,\n",
    "    maxPixels=1e9\n",
    "    )\n",
    "    return feature.set(crop_stats)\n",
    "\n",
    "# Apply the function to each feature in the administrative boundaries collection\n",
    "cropland_exposed_by_division = divisions.map(calculate_cropland_exposed)\n",
    "\n",
    "# Fetch results\n",
    "cropland_result = cropland_exposed_by_division.getInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 16.4. Calculate estimated exposed urban area per administrative division\n",
    "\n",
    "And lastly, the same is done for exposed urban areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to calculate the sum of exposed urban area within each divsion\n",
    "def calculate_urban_exposed(feature):\n",
    "    urban_stats = urban_pixelarea.reduceRegion(\n",
    "        reducer = ee.Reducer.sum(),\n",
    "        geometry = feature.geometry(),\n",
    "        scale = 100,\n",
    "        maxPixels = 1e9,\n",
    "        bestEffort = True\n",
    "    )\n",
    "    return feature.set(urban_stats)\n",
    "\n",
    "# Apply the function to each feature in the administrative boundaries collection\n",
    "urban_exposed_by_division = divisions.map(calculate_urban_exposed)\n",
    "\n",
    "# Fetch results\n",
    "urban_result = urban_exposed_by_division.getInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 16.5. Create pandas dataframe from results\n",
    "\n",
    "We can compile the results from the previous steps into a pandas dataframe to allow for easy analysis and visualisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract data from results\n",
    "data = {}\n",
    "for feature in flood_result['features']:\n",
    "    division = feature['properties']['ADM2_EN']\n",
    "    flood_extent = feature['properties']['VH']\n",
    "    adm1_en = feature['properties']['ADM1_EN']\n",
    "    adm1_pcode = feature['properties']['ADM1_PCODE']\n",
    "    adm2_pcode = feature['properties']['ADM2_PCODE']\n",
    "    data[division] = {\n",
    "        'ADM2_PCODE': adm2_pcode,\n",
    "        'ADM1_EN': adm1_en,\n",
    "        'ADM1_PCODE': adm1_pcode,\n",
    "        'flood_extent_km2': round(flood_extent / 1e6),\n",
    "        'flood_extent_ha': round(flood_extent / 1e4)\n",
    "    }\n",
    "\n",
    "for feature in pop_result['features']:\n",
    "    division = feature['properties']['ADM2_EN']\n",
    "    pop_exposed = feature['properties']['population_count']\n",
    "    data[division]['exposed_population'] = round(pop_exposed)\n",
    "\n",
    "for feature in cropland_result['features']:\n",
    "    division = feature['properties']['ADM2_EN']\n",
    "    crop_exposed = feature['properties']['VH']\n",
    "    data[division].update({\n",
    "        'exposed_cropland_km2': round(crop_exposed / 1e6),\n",
    "        'exposed_cropland_ha': round(crop_exposed / 1e4)  \n",
    "    })\n",
    "\n",
    "for feature in urban_result['features']:\n",
    "    division = feature['properties']['ADM2_EN']\n",
    "    urban_exposed = feature['properties']['VH']\n",
    "    data[division].update({\n",
    "        'exposed_urban_km2': round(urban_exposed / 1e6),\n",
    "        'exposed_urban_ha': round(urban_exposed / 1e4)  \n",
    "    })\n",
    "\n",
    "# Create dataframe\n",
    "df = pd.DataFrame.from_dict(data, orient='index')\n",
    "df.reset_index(inplace=True)\n",
    "df.rename(columns={'index': 'ADM2_EN'}, inplace=True)\n",
    "\n",
    "# Show the dataframe\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 17. Create charts\n",
    "\n",
    "Now that we've successfully compiled our results into a dataframe, we can delve deeper into our analysis by creating some insightful charts. These visualisations will focus on understanding the impact of the flooding across administrative divisions, providing a clearer picture at a localized level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 17.1. Chart estimated flooding extent by division\n",
    "\n",
    "Let's create a horizontal bar chart depicting estimated flood extent across ADM2 administrative divisions. The dataframe we created containing flooding statistics is filtered to focus on divisions with significant flooding extent (>=20km²). This filtered data is then sorted in ascending order of flood extent for clarity in the chart. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter out divisions with flood extent less than 20km²\n",
    "flood_filtered_df = df[df['flood_extent_km2'] >= 20]\n",
    "\n",
    "# Sort the filtered dataframe in descending order of flood extent\n",
    "flood_sorted_df = flood_filtered_df.sort_values(by='flood_extent_km2', ascending=True)\n",
    "\n",
    "# Extract the sorted data\n",
    "divisions = flood_sorted_df['ADM2_EN']\n",
    "flood_extents = flood_sorted_df['flood_extent_km2']\n",
    "\n",
    "# Create a horizontal bar chart\n",
    "plt.figure(figsize=(16, 10))\n",
    "bars = plt.barh(divisions, flood_extents, color='#0076be')\n",
    "\n",
    "# Make axis labels and title\n",
    "plt.ylabel('Administrative Division', fontweight='bold')\n",
    "plt.xlabel('Estimated Flooded Area (km²)', fontweight='bold')\n",
    "plt.title('Estimated Flooding Extent by Administrative Division (km²)', fontweight='bold')\n",
    "\n",
    "# Adjust x and y axis ticks and add grid\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_major_locator(ticker.MaxNLocator(10))\n",
    "plt.grid(True, which='both', axis='x', linestyle='--', linewidth=0.2, color='gray')\n",
    "\n",
    "# Add labels to the bars\n",
    "for bar in bars:\n",
    "    plt.text(bar.get_width(), bar.get_y() + bar.get_height()/2, round(bar.get_width(), 2), \n",
    "             va='center', ha='left')\n",
    "    \n",
    "# Add a note to the graph\n",
    "note = \"Note: Data represents divisions with >= 20 km² flooding extent\"\n",
    "plt.text(0.56, -5, note, fontsize=10, ha='center', va='bottom', color='gray')\n",
    "\n",
    "# Show the chart\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 17.2. Chart estimated exposed population by division\n",
    "\n",
    "Let's do the same for exposed population. The dataframe containing exposed population statistics is filtered to focus on divisions with more than 100 people exposed. This filtered data is then sorted in ascending order of population exposure for clarity in the chart. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter out divisions with exposed population of less than 100\n",
    "pop_filtered_df = df[df['exposed_population'] >= 100]\n",
    "\n",
    "# Sort the filtered dataframe in descending order of exposed cropland\n",
    "pop_sorted_df = pop_filtered_df.sort_values(by='exposed_population', ascending=True)\n",
    "\n",
    "# Extract the sorted data\n",
    "divisions = pop_sorted_df['ADM2_EN']\n",
    "exposed_pop = pop_sorted_df['exposed_population']\n",
    "\n",
    "# Create a horizontal bar chart\n",
    "plt.figure(figsize=(16, 10))\n",
    "bars = plt.barh(divisions, exposed_pop, color='#f9a73e')\n",
    "\n",
    "# Make axis labels and title\n",
    "plt.ylabel('Administrative Division', fontweight='bold')\n",
    "plt.xlabel('Estimated Exposed Population', fontweight='bold')\n",
    "plt.title('Estimated Exposed Population by Administrative Division', fontweight='bold')\n",
    "\n",
    "# Adjust x and y axis ticks and add grid\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_major_locator(ticker.MaxNLocator(10))\n",
    "plt.grid(True, which='both', axis='x', linestyle='--', linewidth=0.2, color='gray')\n",
    "\n",
    "# Add labels to the bars\n",
    "for bar in bars:\n",
    "    plt.text(bar.get_width(), bar.get_y() + bar.get_height()/2, round(bar.get_width(), 2), \n",
    "             va='center', ha='left')\n",
    "    \n",
    "# Add a note to the graph\n",
    "note = \"Note: Data represents divisions with >= 100 people exposed\"\n",
    "plt.text(0.56, -5, note, fontsize=10, ha='center', va='bottom', color='gray')\n",
    "\n",
    "# Show the chart\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 17.3. Chart estimated exposed cropland by division\n",
    "\n",
    "Again, let's do the same for exposed cropland."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter out divisions with less than 20km² exposed cropland\n",
    "crop_filtered_df = df[df['exposed_cropland_km2'] >= 20]\n",
    "\n",
    "# Sort the filtered dataframe in descending order of exposed cropland\n",
    "crop_sorted_df = crop_filtered_df.sort_values(by='exposed_cropland_km2', ascending=True)\n",
    "\n",
    "# Extract the sorted data\n",
    "divisions = crop_sorted_df['ADM2_EN']\n",
    "flood_extents = crop_sorted_df['exposed_cropland_km2']\n",
    "\n",
    "# Create a horizontal bar chart\n",
    "plt.figure(figsize=(16, 10))\n",
    "bars = plt.barh(divisions, flood_extents, color='#48bf91')\n",
    "\n",
    "# Make axis labels and title\n",
    "plt.ylabel('Administrative Division', fontweight='bold')\n",
    "plt.xlabel('Estimated Exposed Cropland (km²)', fontweight='bold')\n",
    "plt.title('Estimated Exposed Cropland by Administrative Division (km²)', fontweight='bold')\n",
    "\n",
    "# Adjust x and y axis ticks and add grid\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_major_locator(ticker.MaxNLocator(10))\n",
    "plt.grid(True, which='both', axis='x', linestyle='--', linewidth=0.2, color='gray')\n",
    "\n",
    "# Add labels to the bars\n",
    "for bar in bars:\n",
    "    plt.text(bar.get_width(), bar.get_y() + bar.get_height()/2, round(bar.get_width(), 2), \n",
    "             va='center', ha='left')\n",
    "    \n",
    "# Add a note to the graph\n",
    "note = \"Note: Data represents divisions with >= 20 km² of exposed cropland\"\n",
    "plt.text(0.56, -5, note, fontsize=10, ha='center', va='bottom', color='gray')\n",
    "\n",
    "# Show the chart\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 17.4. Chart the top 10 most affected divisions in terms of flood extent and exposed cropland\n",
    "\n",
    "Here we create a horizontal bar chart that compares flood extents and exposed cropland areas in the top 10 most affected administrative divisions. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort the dataframe by flood extent and select the top 10\n",
    "top10_sorted_df = df.sort_values(by='flood_extent_km2', ascending=False).head(10)\n",
    "\n",
    "# Extract the sorted data\n",
    "divisions = top10_sorted_df['ADM2_EN']\n",
    "flood_extents = top10_sorted_df['flood_extent_km2']\n",
    "cropland_extents = top10_sorted_df['exposed_cropland_km2']\n",
    "\n",
    "# Number of divisions\n",
    "n_divisions = len(divisions)\n",
    "\n",
    "# Create a horizontal bar chart\n",
    "fig, ax = plt.subplots(figsize=(16, 10))\n",
    "\n",
    "# Calculate bar positions\n",
    "bar_width = 0.4\n",
    "indices = np.arange(n_divisions)\n",
    "\n",
    "# Plot the bars for each metric\n",
    "ax.barh(indices - bar_width/2, flood_extents, height=bar_width, color='#0076be', label='Flood Extent')\n",
    "ax.barh(indices + bar_width/2, cropland_extents, height=bar_width, color='#48bf91', label='Exposed Cropland')\n",
    "\n",
    "# Set y-axis labels, ticks and grid\n",
    "ax.set(yticks=indices, yticklabels=divisions)\n",
    "ax.invert_yaxis() \n",
    "ax.xaxis.set_major_locator(ticker.MaxNLocator(10))\n",
    "plt.grid(True, which='both', axis='x', linestyle='--', linewidth=0.2, color='gray')\n",
    "\n",
    "# Make axis labels and title\n",
    "plt.xlabel('Area (km²)', fontweight='bold')\n",
    "plt.title('10 Most Affected Divisions: Comparison of Flood Extent and Exposed Cropland (km²)', fontweight='bold')\n",
    "plt.legend()\n",
    "\n",
    "# Show the chart\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 18. Join results dataframe to ADM2 geodataframe\n",
    "\n",
    "Let's create a geodataframe with our flooding statistics dataframe. To do this we import a geojson file with the adm2 boundaries and join our flooding statistics dataframe based on the `ADM2_PCODE` columns in both data sets. We can also do a bit of data wrangling to clean up the resulting geodataframe. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to GeoJSON file\n",
    "gjsn_path = \"inputs/pak_adm2.geojson\"\n",
    "\n",
    "# Read GeoJSON file using Geopandas\n",
    "gdf = gpd.read_file(gjsn_path)\n",
    "\n",
    "# Join the DataFrame and GeoDataFrame\n",
    "joined_df = pd.merge(df, gdf, on='ADM2_PCODE', how='left')\n",
    "\n",
    "# Drop unnecessary columns\n",
    "drop_list = ['ADM2_REF', 'ADM2ALT1EN', 'ADM2ALT2EN', 'date', 'validOn', 'validTo']\n",
    "joined_df.drop(columns=drop_list, inplace=True, errors='ignore')\n",
    "joined_df = joined_df[joined_df.columns.drop(list(joined_df.filter(regex='_y')))]\n",
    "\n",
    "# Rename columns\n",
    "joined_df.rename(columns={'ADM2_EN_x': 'ADM2_EN', 'ADM1_EN_x': 'ADM1_EN', 'ADM1_PCODE_x': 'ADM1_PCODE'}, inplace=True)\n",
    "\n",
    "# Convert joined DataFrame to GeoDataFrame\n",
    "final_gdf = gpd.GeoDataFrame(joined_df, geometry='geometry')\n",
    "\n",
    "# Display GeoDataFrame\n",
    "final_gdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 19. Export to GeoJSON and CSV\n",
    "\n",
    "We can export our results to a geojson file and a CSV for further analysis or map-making in GIS software such as QGIS or ArcGIS. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export GeoDataFrame to GeoJSON\n",
    "final_gdf.to_file(\"outputs/pak_adm2_flooding.geojson\", driver='GeoJSON')\n",
    "\n",
    "# Drop unnecessary columns for export to CSV\n",
    "drop_list = ['Shape_Leng', 'geometry', 'Shape_Area']\n",
    "joined_df.drop(columns=drop_list, inplace=True, errors='ignore')\n",
    "final_csv = joined_df\n",
    "\n",
    "# Export to CSV\n",
    "final_csv.to_csv(\"outputs/pak_adm2_flooding.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20. Export images to GeoTIFF\n",
    "\n",
    "Finally, let's set up a series of exports to Google Drive, converting our flood extent, cropland, population, and urban images into GeoTIFF format. These exports will take quite a long time. Check your Drive account in a few hours. You can view the export progress in the [GEE tasks panel](https://code.earthengine.google.com/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export flood extent to Google Drive as GeoTIFF\n",
    "export_task1 = ee.batch.Export.image.toDrive(image=flooded,\n",
    "                                            description='Export Flood Extent',\n",
    "                                            folder='Pakistan_Flood_Analysis_GeoTIFFs',\n",
    "                                            fileNamePrefix='Flood_Extent',\n",
    "                                            scale=10,\n",
    "                                            region=aoi, \n",
    "                                            maxPixels=1e13,\n",
    "                                            crs='EPSG:4326',\n",
    "                                            fileFormat='GeoTIFF',\n",
    "                                            formatOptions={\n",
    "                                                    'cloudOptimized': True\n",
    "                                                })\n",
    "\n",
    "# Export population exposed to Google Drive as GeoTIFF\n",
    "export_task2 = ee.batch.Export.image.toDrive(image=population_exposed,\n",
    "                                            description='Export Population Exposed',\n",
    "                                            folder='Pakistan_Flood_Analysis_GeoTIFFs',\n",
    "                                            fileNamePrefix='Population_Exposed',\n",
    "                                            scale=100,\n",
    "                                            region=aoi, \n",
    "                                            maxPixels=1e13,\n",
    "                                            crs='EPSG:4326',\n",
    "                                            fileFormat='GeoTIFF',\n",
    "                                            formatOptions={\n",
    "                                                    'cloudOptimized': True\n",
    "                                                })\n",
    "\n",
    "# Export cropland exposed to Google Drive as GeoTIFF\n",
    "export_task3 = ee.batch.Export.image.toDrive(image=cropland_exposed,\n",
    "                                            description='Export Cropland Exposed',\n",
    "                                            folder='Pakistan_Flood_Analysis_GeoTIFFs',\n",
    "                                            fileNamePrefix='Cropland_Exposed',\n",
    "                                            scale=100,\n",
    "                                            region=aoi, \n",
    "                                            maxPixels=1e13,\n",
    "                                            crs='EPSG:4326',\n",
    "                                            fileFormat='GeoTIFF',\n",
    "                                            formatOptions={\n",
    "                                                    'cloudOptimized': True\n",
    "                                                })\n",
    "\n",
    "# Export urban exposed to Google Drive as GeoTIFF\n",
    "export_task4 = ee.batch.Export.image.toDrive(image=urban_exposed,\n",
    "                                            description='Export Urban Exposed',\n",
    "                                            folder='Pakistan_Flood_Analysis_GeoTIFFs',\n",
    "                                            fileNamePrefix='Urban_Exposed',\n",
    "                                            scale=100,\n",
    "                                            region=aoi, \n",
    "                                            maxPixels=1e13,\n",
    "                                            crs='EPSG:4326',\n",
    "                                            fileFormat='GeoTIFF',\n",
    "                                            formatOptions={\n",
    "                                                    'cloudOptimized': True\n",
    "                                                })\n",
    "\n",
    "# Start the export tasks\n",
    "export_task1.start()\n",
    "export_task2.start()\n",
    "export_task3.start()\n",
    "export_task4.start()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "eo",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}