tutorial/06_Validate_Selected_LGAs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "13e0fadc-9090-44fe-b1e8-4ccc8e66dd07",
   "metadata": {},
   "source": [
    "# Validate Selected LGAs\n",
    "\n",
    "The workflow in this notebook demonstrates how to load and a section of a generated flood map and plot it over the relevant satellite imagery.\n",
    "\n",
    "For this tutorial, we will load the data for the EMSR586 Copernicus activation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af1c5398-4318-46e0-bb7b-2f7ff38cf3c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Necessary imports\n",
    "import os\n",
    "os.environ['USE_PYGEOS'] = '0'\n",
    "import gc\n",
    "import numpy as np\n",
    "from itertools import product\n",
    "from datetime import datetime\n",
    "from zoneinfo import ZoneInfo\n",
    "from dotenv import load_dotenv\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors\n",
    "import ee\n",
    "from georeader.readers import ee_query\n",
    "import folium\n",
    "import geemap.foliumap as geemap\n",
    "import shapely\n",
    "from shapely.geometry import box, Point, MultiPoint\n",
    "from shapely.ops import nearest_points\n",
    "\n",
    "\n",
    "from ml4floods.data import utils\n",
    "from ml4floods.visualization import plot_utils\n",
    "\n",
    "from db_utils import DB\n",
    "\n",
    "# Uncomment this to suppress deprecation warnings\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning) \n",
    "from shapely.errors import ShapelyDeprecationWarning\n",
    "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning) \n",
    "\n",
    "# Set bucket will not be requester pays\n",
    "utils.REQUESTER_PAYS_DEFAULT = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40bfce80",
   "metadata": {},
   "source": [
    "## Load environment and project details\n",
    "\n",
    "As with the other notebooks, we load credentials and project details from a hidden ```.env``` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36c177f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load environment variables (including path to credentials) from '.env' file\n",
    "env_file_path = \"../.env\"\n",
    "\n",
    "# Uncomment for alternative version for Windows (r\"\" indicates raw string)\n",
    "#env_file_path = r\"C:/Users/User/floodmapper/.env\"\n",
    "\n",
    "assert load_dotenv(dotenv_path=env_file_path) == True, \"[ERR] failed to load environment!\"\n",
    "assert \"GOOGLE_APPLICATION_CREDENTIALS\" in os.environ, \"[ERR] missing $GOOGLE_APPLICATION_CREDENTIAL!\"\n",
    "assert \"GS_USER_PROJECT\" in os.environ, \"[ERR] missing $GS_USER_PROJECT!\"\n",
    "key_file_path = os.environ[\"GOOGLE_APPLICATION_CREDENTIALS\"]\n",
    "assert os.path.exists(key_file_path), f\"[ERR] Google credential key file does not exist: \\n{key_file_path} \"\n",
    "assert \"ML4FLOODS_BASE_DIR\" in os.environ, \"[ERR] missing $ML4FLOODS_BASE_DIR!\"\n",
    "base_path = os.environ[\"ML4FLOODS_BASE_DIR\"]\n",
    "assert os.path.exists(base_path), f\"[ERR] base path does not exist: \\n{base_path} \"\n",
    "bucket_name = os.environ[\"BUCKET_URI\"]\n",
    "assert bucket_name is not None and bucket_name != \"\", f\"Bucket name not defined {bucket_name}\"\n",
    "print(\"[INFO] Successfully loaded FloodMapper environment.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f74f609f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Connect to the database (point to the .env file for credentials)\n",
    "db_conn = DB(env_file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d91dc02",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the Google Earth Engine connection.\n",
    "# Follow instructions on login prompt, if required.\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39922e40",
   "metadata": {},
   "source": [
    "## Retrieve the session parameters from the database\n",
    "\n",
    "Edit the following cell to set the session name and retrieve the flood map parameters, including mapping grid and affected LGAs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0605733b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# EDIT THIS CELL: All work is conducted under a unique session name\n",
    "session_name = \"EMSR586\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c232acb-cd0b-4f41-b61e-de5e788b99d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "tz = ZoneInfo(\"UTC\")\n",
    "midnight = datetime.min.time()\n",
    "\n",
    "# Query the floodmapping parameters from the DB\n",
    "query = (f\"SELECT flood_date_start, flood_date_end, \"\n",
    "         f\"ref_date_start, ref_date_end, bucket_uri \"\n",
    "         f\"FROM session_info \"\n",
    "         f\"WHERE session = %s;\")\n",
    "data = (session_name,)\n",
    "session_df = db_conn.run_query(query, data, fetch=True)\n",
    "flood_start_date = session_df.iloc[0][\"flood_date_start\"]\n",
    "flood_start_date = datetime.combine(flood_start_date, midnight).replace(tzinfo=tz)\n",
    "flood_end_date = session_df.iloc[0][\"flood_date_end\"]\n",
    "flood_end_date = datetime.combine(flood_end_date, midnight).replace(tzinfo=tz)\n",
    "#ref_start_date = session_df.iloc[0][\"ref_date_start\"]\n",
    "#ref_start_date = datetime.combine(ref_start_date, midnight).replace(tzinfo=tz)\n",
    "#ref_end_date = session_df.iloc[0][\"ref_date_end\"]\n",
    "#ref_end_date = datetime.combine(ref_end_date, midnight).replace(tzinfo=tz)\n",
    "bucket_uri = session_df.iloc[0][\"bucket_uri\"]\n",
    "\n",
    "# Query the selected grid positions and LGAs\n",
    "query = (f\"SELECT sp.patch_name, ST_AsText(gr.geometry), gr.lga_name22 \"\n",
    "         f\"FROM session_patches sp \"\n",
    "         f\"INNER JOIN grid_loc gr \"\n",
    "         f\"ON sp.patch_name = gr.patch_name \"\n",
    "         f\"WHERE sp.session = %s ;\")\n",
    "data = (session_name,)\n",
    "grid_sel_df = db_conn.run_query(query, data, fetch=True)\n",
    "\n",
    "# Format the results into a correct GeoDataFrame\n",
    "grid_sel_df['geometry'] = gpd.GeoSeries.from_wkt(grid_sel_df['st_astext'])\n",
    "grid_sel_df.drop(['st_astext'], axis=1, inplace = True)\n",
    "grid_sel_gdf = gpd.GeoDataFrame(grid_sel_df, geometry='geometry', crs=\"EPSG:4326\")\n",
    "grid_sel_gdf.drop_duplicates(subset=[\"patch_name\"], inplace=True)\n",
    "print(f\"[INFO] {len(grid_sel_gdf)} grid patches selected.\")\n",
    "\n",
    "# Query the affected LGA shapes\n",
    "lgas_sel_lst = grid_sel_df.lga_name22.unique().tolist()\n",
    "query = (f\"SELECT DISTINCT lga_name22, ST_AsText(geometry_col) \"\n",
    "         f\"FROM lgas_info \"\n",
    "         f\"WHERE lga_name22 IN %s ;\")\n",
    "data = (tuple(lgas_sel_lst),)\n",
    "lgas_sel_df = db_conn.run_query(query, data, fetch=True)\n",
    "\n",
    "# Format the results into a correct GeoDataFrame\n",
    "lgas_sel_df['geometry'] = gpd.GeoSeries.from_wkt(lgas_sel_df['st_astext'])\n",
    "lgas_sel_df.drop(['st_astext'], axis=1, inplace = True)\n",
    "lgas_sel_gdf = gpd.GeoDataFrame(lgas_sel_df, geometry='geometry', crs=\"EPSG:4326\")\n",
    "print(f\"[INFO] {len(lgas_sel_gdf)} LGAs affected.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6ca080f",
   "metadata": {},
   "source": [
    "## Plot the mapping grid and overlapping LGAs \n",
    "\n",
    "Use this map as an aid to choose which LGA to validate. Display LGA names by hovering the mouse cursor over the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef7eabac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Overplot the LGAs on the selected mapping grid\n",
    "m = grid_sel_gdf.explore(style_kwds={\"fillOpacity\": 0.00}, color=\"red\", \n",
    "                         name=\"Grid Patches\", highlight=False)\n",
    "lgas_sel_gdf.explore(m=m, style_kwds={\"fillOpacity\": 0.01}, name=\"LGAs\")\n",
    "folium.LayerControl(collapsed=False).add_to(m)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bb59b24",
   "metadata": {},
   "source": [
    "## Choose the LGAs to validate\n",
    "\n",
    "Edit the next cell to choose a subset of LGAs over which to load the flood extent map. Set an empty list to load the full flood map.\n",
    "\n",
    "**WARNING: loading too large an area will likely fail due to memory limits.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9f366ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# EDIT THIS CELL: Set the list of LGAs to load or leave blank [] to display whole map\n",
    "selected_lgas = []\n",
    "#selected_lgas = [\"Newcastle\", \"Maitland\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc2fb14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the outer bounds of the chosen area\n",
    "if len(selected_lgas) == 0:\n",
    "    # If no LGAs specified, select the bounds of all patches\n",
    "    print(\"[INFO] Selecting bounds of all patches.\")\n",
    "    boundary_box = box(*grid_sel_gdf.total_bounds)\n",
    "else:\n",
    "    # Or select the bounds of chosen LGAs\n",
    "    print(\"[INFO] Selecting bounds of chosen LGAs.\")\n",
    "    lgas_filtered_gdf = lgas_sel_gdf[lgas_sel_gdf.lga_name22.isin(selected_lgas)]\n",
    "    boundary_box = box(*lgas_filtered_gdf.total_bounds)\n",
    "boundary_gdf = gpd.GeoDataFrame(geometry=[boundary_box], crs=grid_sel_gdf.crs)\n",
    "print(\"[INFO] crs is\", grid_sel_gdf.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "185b3d05",
   "metadata": {},
   "source": [
    "**Select the patch nearest the centre to plot a time-series.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b054631",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use patch centroids to select patch_name\n",
    "grid_sel_gdf[\"cents\"] = grid_sel_gdf.centroid\n",
    "patch_cents = MultiPoint(grid_sel_gdf[\"cents\"].to_list())\n",
    "nearest_point = nearest_points(boundary_gdf.centroid, patch_cents)[1]\n",
    "nearest_idx = grid_sel_gdf[\"cents\"] == nearest_point.values[0]\n",
    "nearest_patch_name = grid_sel_gdf[nearest_idx].patch_name.values[0]\n",
    "print(f\"[INFO] Closest patch to centre is {nearest_patch_name}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a98b37f1",
   "metadata": {},
   "source": [
    "## Visualise data availability and cloud-cover"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18044ab6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the date and cloud cover of available imagery\n",
    "query = (f\"SELECT image_id, satellite, cloud_probability, solarday, status \"\n",
    "         f\"FROM image_downloads \"\n",
    "         f\"WHERE patch_name = %s;\")\n",
    "data = (nearest_patch_name,)\n",
    "downloads_df = db_conn.run_query(query, data, fetch=True)\n",
    "print(f\"[INFO] Found {len(downloads_df)} images for {nearest_patch_name}.\")\n",
    "\n",
    "# Plot the time-series\n",
    "fig, ax = plt.subplots(1,1, figsize=(15,5))\n",
    "for sat, downloaded in product([\"Landsat\",\"S2\"], [1, 0]):\n",
    "    label=f\"{sat} ({'' if downloaded else 'NOT '}Downloaded)\"\n",
    "    color = \"C0\" if sat ==\"S2\" else \"C1\"\n",
    "    marker = \"o\" if downloaded else \"x\"\n",
    "    selection_mask = (downloads_df.satellite == sat) & (downloads_df.status == downloaded)\n",
    "    ax.scatter(x=downloads_df[selection_mask].solarday, \n",
    "               y=downloads_df[selection_mask].cloud_probability,\n",
    "               label=label, c=color, marker=marker, s=100)\n",
    "    \n",
    "# Format plot to look nice\n",
    "ax.legend()\n",
    "downloads_df.plot(x=\"solarday\", y=\"cloud_probability\", ax=ax,legend=None)\n",
    "ax.axvspan(flood_start_date, flood_end_date, alpha=0.05)\n",
    "datespan = abs(flood_end_date - flood_start_date)\n",
    "ax.set_xlim(flood_start_date - datespan *.1, flood_end_date + datespan *.1)\n",
    "plt.ylabel(\"Mean Cloud Coverage (%)\")\n",
    "plt.xlabel(\"Date\")\n",
    "plt.title(f\"Available Data During Flood for Patch {nearest_patch_name}.\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "994ea01d",
   "metadata": {},
   "source": [
    "## Load the floodmap of the selected area\n",
    "\n",
    "First load the floodmap and clip to the selected area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d96d53b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Query the floodmap name on GCP\n",
    "query = (f\"SELECT data_path FROM postproc_spatial \"\n",
    "         f\"WHERE session = %s \"\n",
    "         f\"AND mode = %s;\")\n",
    "data = (session_name, \"flood\")\n",
    "tmp_df = db_conn.run_query(query, data, fetch=True)\n",
    "floodmap_path = tmp_df.iloc[0][\"data_path\"]\n",
    "\n",
    "# Load the flood map from GCP\n",
    "# Convert to projected CRS EPSG3857\n",
    "print(f\"[INFO] Loading the floodmap:\\n\\t{floodmap_path}\")\n",
    "floodmap = utils.read_geojson_from_gcp(floodmap_path).to_crs(epsg=3857)\n",
    "\n",
    "# Clip the floodmap at the boundary and delete the original\n",
    "floodmap_clip = floodmap.clip(boundary_gdf.to_crs(epsg=3857))\n",
    "del floodmap\n",
    "gc.collect()\n",
    "\n",
    "# Plot the floodmap\n",
    "plot_utils.plot_floodmap(floodmap_clip) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4529b5df",
   "metadata": {},
   "source": [
    "## Query the available satellite imagery during the flood period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea8bf893",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Format the boundary into a multipolygon\n",
    "region = boundary_gdf.unary_union\n",
    "\n",
    "# Run a GEE query for Landsat and Sentinel-2 data.\n",
    "#  producttype can be 'both', 'S2', \"Landsat\", \"L8\" or \"L9\".\n",
    "#  add_s2cloudless adds a column that indicates if the s2cloudless image is available.\n",
    "flood_images_gee, flood_collection = ee_query.query(\n",
    "    area=region, \n",
    "    date_start=flood_start_date, \n",
    "    date_end=flood_end_date,                                                   \n",
    "    producttype=\"both\", \n",
    "    return_collection=True, \n",
    "    add_s2cloudless=True)\n",
    "\n",
    "# Print data about the available images\n",
    "num_images = flood_images_gee.shape[0]\n",
    "print(f\"[INFO] Found {num_images} flooding images on archive.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a08ee7c",
   "metadata": {},
   "source": [
    "## Plot the floodmap over the satellite data\n",
    "\n",
    "**Building the map in the next cell may take some time.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0531aef8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter the map for 'flood_trace' and 'water' polygons only\n",
    "floodmap = floodmap_clip.loc[floodmap_clip['class'].apply(lambda x: x in ['flood_trace', 'water'])]\n",
    "categories = floodmap[\"class\"].unique()\n",
    "print(\"[INFO] Categories in file: \", categories)\n",
    "\n",
    "# Set the colourmap for the maps\n",
    "COLORS = {\n",
    "    'cloud': \"gray\",\n",
    "    'flood_trace': \"turquoise\",\n",
    "    'water': \"blue\"\n",
    "}\n",
    "cmap = matplotlib.colors.ListedColormap([COLORS[b] for b in categories])\n",
    "\n",
    "# Intialise the OpenStreetMap base layer\n",
    "m = geemap.Map(location=region.centroid.coords[0][-1::-1], zoom_start=10)\n",
    "\n",
    "# Add the satellite data for each day\n",
    "for (day, satellite), images_day in flood_images_gee.groupby([\"solarday\", \"satellite\"]):    \n",
    "    image_col_day_sat = flood_collection.filter(ee.Filter.inList(\"title\", images_day.index.tolist()))    \n",
    "    bands = [\"B11\",\"B8\",\"B4\"] if satellite.startswith(\"S2\") else [\"B6\",\"B5\",\"B4\"]\n",
    "    m.addLayer(image_col_day_sat, \n",
    "               {\"min\":0, \"max\":3000 if satellite.startswith(\"S2\") else 0.3, \"bands\": bands},\n",
    "               f\"{satellite}: {day}\",\n",
    "               False)\n",
    "    \n",
    "# Overlay the floodmap and the boundary\n",
    "floodmap.explore(m=m, column=\"class\", cmap=cmap, categories=categories, \n",
    "                 name=\"Flood Extent Map\", style_kwds={\"fillOpacity\": 0.3, \"weight\": 1, \"opacity\": 0.5})\n",
    "boundary_gdf.explore(m=m, style_kwds={\"fillOpacity\": 0.0}, color=\"black\", \n",
    "                     name=\"Boundary\", highlight=False)\n",
    "\n",
    "# Add the layer control\n",
    "folium.LayerControl(collapsed=False).add_to(m)\n",
    "print(\"[INFO] Map ready: execute next cell to display.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3752296f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the map\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad377991",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (floodmapper)",
   "language": "python",
   "name": "floodmapper"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}