3. gen-flood-raster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Config Cell\n",
    "## [Jump to outputs](#Outputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add shapefile\n",
    "shp = \"input/NER_ADM0.geojson\"\n",
    "buffer = 0.2  # boundary buffer\n",
    "cell_size = 0.09  # grid cell size\n",
    "\n",
    "# Run 1. aoi-threshold.ipynb to get the value of th_aoi and store it here.\n",
    "th_aoi = -21\n",
    "# Set test_cells value if you need to test with initial 'n' cells instead of the entire grid\n",
    "test_cells = []  # [start, end] Run from start(excl) to end(incl).\n",
    "# Set run_checks to True if you want to run checks\n",
    "run_checks = False\n",
    "\n",
    "# Create folder ids by copying the ID from the g-drive folder url. If not using a particular ID, set as None\n",
    "F_MN_FID = \"1q7S4tYQ9MHRDJGpyIwRqEF_ViiAS8Qb6\"\n",
    "F_MD_FID = None\n",
    "PRF_MN_FID = None\n",
    "PRF_MD_FID = None\n",
    "POF_MN_FID = None\n",
    "POF_MD_FID = None\n",
    "FE_MN_FID = None\n",
    "FE_MD_FID = None\n",
    "MF_MX_FID = None\n",
    "\n",
    "ERR_FOLDER_ID = \"168qr6gMAY0q53Oq7t5SuaArEDSy2wxq4\"\n",
    "\n",
    "process_name = \"NER_10_1123-1124\"\n",
    "\n",
    "# Define main time period of analysis\n",
    "timerange = (\"2023-11\", \"2024-11\")\n",
    "# Define sub-periods of analysis - should be within main time period\n",
    "pre_flood = []\n",
    "flood = [\"2023-11\", \"2024-11\"]\n",
    "post_flood = []\n",
    "max_flood = []\n",
    "# Analysis periods\n",
    "a_periods = [\"flood\"]\n",
    "# Required Analysis measurements\n",
    "a_measures = [\"mean\"]  # a_measures = [\"mean\", \"median\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Init Requirements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# !pip install python-dotenv\n",
    "# load_dotenv()\n",
    "# from dotenv import load_dotenv\n",
    "\n",
    "import os, glob, warnings, datacube, rasterio, folium, json\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import rioxarray as rio\n",
    "import pandas as pd\n",
    "from rasterio.merge import merge\n",
    "from rasterio.plot import show\n",
    "from shapely.geometry import Point\n",
    "from shapely.geometry import Polygon\n",
    "from datetime import date\n",
    "from dateutil.relativedelta import relativedelta\n",
    "\n",
    "\n",
    "from scipy.ndimage import uniform_filter\n",
    "from scipy.ndimage import variance\n",
    "from skimage.filters import threshold_minimum\n",
    "from datacube.utils.geometry import Geometry\n",
    "\n",
    "from deafrica_tools.spatial import xr_rasterize\n",
    "from deafrica_tools.datahandling import load_ard\n",
    "from deafrica_tools.plotting import display_map, rgb\n",
    "from deafrica_tools.areaofinterest import define_area\n",
    "\n",
    "from typing import Literal\n",
    "\n",
    "from IPython.display import clear_output\n",
    "from IPython.display import display\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Init Google Drive (GCS) module\n",
    "from tools.gdrive import GDrive\n",
    "\n",
    "gd = GDrive()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gd.get_storage()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connect to the datacube\n",
    "\n",
    "Connect to the datacube so we can access DEA data.\n",
    "The `app` parameter is a unique name for the analysis which is based on the notebook file name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dc = datacube.Datacube(app=\"Radar_water_detection\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter and Classifier Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to apply lee filtering on S1 image. Speckle Filter\n",
    "def lee_filter(da, size):\n",
    "    \"\"\"\n",
    "    Apply lee filter of specified window size.\n",
    "    Adapted from https://stackoverflow.com/questions/39785970/speckle-lee-filter-in-python\n",
    "\n",
    "    \"\"\"\n",
    "    img = da.values\n",
    "    img_mean = uniform_filter(img, size)\n",
    "    img_sqr_mean = uniform_filter(img**2, size)\n",
    "    img_variance = img_sqr_mean - img_mean**2\n",
    "\n",
    "    overall_variance = variance(img)\n",
    "\n",
    "    img_weights = img_variance / (img_variance + overall_variance)\n",
    "    img_output = img_mean + img_weights * (img - img_mean)\n",
    "\n",
    "    return img_output\n",
    "\n",
    "\n",
    "# Classifier Function\n",
    "def S1_water_classifier(da, threshold):\n",
    "    water_data_array = da < threshold\n",
    "    return water_data_array.to_dataset(name=\"s1_water\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operational Funtions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_months(timerange: tuple) -> list[date]:\n",
    "    \"\"\"\n",
    "    Converts the timerange into a list of months.\n",
    "\n",
    "    Parameters:\n",
    "    timerange (\"YYYY-MM\", \"YYYY-MM\"): tuple, required\n",
    "        Timerange of the analysis. Start month to end month.\n",
    "\n",
    "    Returns:\n",
    "    months: list[date]\n",
    "    \"\"\"\n",
    "    syear, smonth = timerange[0].split(\"-\")\n",
    "    syear = int(syear)\n",
    "    smonth = int(smonth)\n",
    "    start = date(syear, smonth, 30)\n",
    "\n",
    "    eyear, emonth = timerange[1].split(\"-\")\n",
    "    eyear = int(eyear)\n",
    "    emonth = int(emonth)\n",
    "    end = date(eyear, emonth, 30)\n",
    "    months = []\n",
    "    while start <= end:\n",
    "        months.append(start.strftime(\"%Y-%m\"))\n",
    "        start = start + relativedelta(months=+1)\n",
    "\n",
    "    return months"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate raster outputs\n",
    "def gen_output(\n",
    "    DS: xr.DataArray,\n",
    "    i: list,\n",
    "    poly: Geometry,\n",
    "    cell: int,\n",
    "    aoi_m: list,\n",
    "    period: Literal[\"preflood\", \"flood\", \"postflood\", \"flood_extents\", \"maxflood\"],\n",
    "    measure: Literal[\"mean\", \"median\", \"max\"],\n",
    ") -> list:\n",
    "    \"\"\"\n",
    "    Generates file names and uploads tp Google Drive\n",
    "\n",
    "    Parameters:\n",
    "    DS: xr.DataArray, required\n",
    "        Dataset to be converted into raster output.\n",
    "    i: list, required\n",
    "        Index of centroid list.\n",
    "    poly: Geometry, required\n",
    "        Polygon feature used to create the DS.\n",
    "    cell: int, required\n",
    "        Cell number being processed.\n",
    "    aoi_m: list, required\n",
    "        List of geojson feature collection (cells) that make up the entire grid.\n",
    "    period: Literal[\"preflood\", \"flood\", \"postflood\", \"flood_extents\"], required\n",
    "        Time period of the process - PRE_FLOOD, FLOOD.\n",
    "    measure: Literal[\"mean\", \"median\"], required\n",
    "        Central tendency measurement of the DS - mean, median.\n",
    "\n",
    "    Returns:\n",
    "    err: list\n",
    "        Error log list of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload.\n",
    "    \"\"\"\n",
    "\n",
    "    err = []\n",
    "\n",
    "    pm_dict = {\n",
    "        \"preflood_mean\": [pre_flood, PRF_MN_FID],\n",
    "        \"preflood_median\": [pre_flood, PRF_MD_FID],\n",
    "        \"flood_mean\": [flood, F_MN_FID],\n",
    "        \"flood_median\": [flood, F_MD_FID],\n",
    "        \"postflood_mean\": [post_flood, POF_MN_FID],\n",
    "        \"postflood_median\": [post_flood, POF_MD_FID],\n",
    "        \"flood_extents_mean\": [timerange, FE_MN_FID],\n",
    "        \"flood_extents_median\": [timerange, FE_MD_FID],\n",
    "        \"maxflood_max\": [timerange, MF_MX_FID],\n",
    "    }\n",
    "\n",
    "    poly_gdf = gpd.GeoDataFrame(geometry=[poly], crs=poly.crs)\n",
    "    lat_range = (poly_gdf.total_bounds[1], poly_gdf.total_bounds[3])\n",
    "    lon_range = (poly_gdf.total_bounds[0], poly_gdf.total_bounds[2])\n",
    "    g = poly_gdf.centroid\n",
    "\n",
    "    data_val = \"CELL_\" + str(i[2]) + \"_{}_{}\".format(period.upper(), measure.upper())\n",
    "    data_name = data_val + \".tif\"\n",
    "    data_out = \"output/{}/\".format(period) + data_name\n",
    "    DS.rio.to_raster(data_out)\n",
    "\n",
    "    # preflood meta\n",
    "    data_dict = {\n",
    "        \"GRID_CELL_ID\": i[2],\n",
    "        \"start_time\": pm_dict[\"{}_{}\".format(period, measure)][0][0],  # pre_flood[0]\n",
    "        \"end_time\": pm_dict[\"{}_{}\".format(period, measure)][0][-1],  # pre_flood[-1]\n",
    "        \"lat\": lat_range,\n",
    "        \"lon\": lon_range,\n",
    "        \"centroid\": \"{}, {}\".format(g.y[0], g.x[0]),\n",
    "        \"crs\": str(poly.crs),\n",
    "    }\n",
    "\n",
    "    text_flie_name = data_val + \"_META.json\"\n",
    "    data_meta_path = \"output/{}/\".format(period) + text_flie_name\n",
    "    with open(data_meta_path, \"w\") as f:\n",
    "        json.dump(data_dict, f)\n",
    "\n",
    "    try:\n",
    "        gd.upload_files(\n",
    "            [data_out, data_meta_path], pm_dict[\"{}_{}\".format(period, measure)][1]\n",
    "        )  # pm_dict['preflood_mean'][1] = FOLDER_ID\n",
    "    except Exception as e:\n",
    "        err.append([g.x[0], g.y[0], i[2], \"U-PRF\"])\n",
    "        print(\n",
    "            \"\\n\\n\"\n",
    "            + \"\\033[31m\"\n",
    "            + \"ERROR UPLOADING GRID CELL ID {} NO.  {}/{} CENTROID ({}, {}). LOGGED CENTROID INFO in e_log\".format(\n",
    "                i[2], cell, len(aoi_m), round(g.y[0], 5), round(g.x[0], 5)\n",
    "            )\n",
    "            + \"\\033[0m\"\n",
    "        )\n",
    "        print(\"UPLOAD ERROR: {}\".format(e))\n",
    "\n",
    "    return err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_elog(e_log: list) -> list:\n",
    "    \"\"\"\n",
    "    Writes a the error log json file and uploads it to the google drive folder ID, if specified.\n",
    "\n",
    "    Parameters:\n",
    "    e_log: list, required\n",
    "        Error log list of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload.\n",
    "\n",
    "    Returns:\n",
    "    e_log: list\n",
    "        Error log list having the same values as the input parameter\n",
    "    \"\"\"\n",
    "    e_log = np.array(e_log)\n",
    "    with open(\"error_centroids.json\", \"w\") as filehandle:\n",
    "        json.dump(e_log.tolist(), filehandle)\n",
    "\n",
    "    # read error log from disk\n",
    "    with open(\"error_centroids.json\") as f:\n",
    "        e_log = json.load(f)\n",
    "    for e in e_log:\n",
    "        e[0] = float(e[0])\n",
    "        e[1] = float(e[1])\n",
    "        e[2] = int(e[2])\n",
    "\n",
    "    try:\n",
    "        gd.upload_files([\"error_centroids.json\"], ERR_FOLDER_ID, False)\n",
    "    except Exception as e:\n",
    "        print(\"FAILED TO UPLOAD ERROR LOG FILE REASON:{}\".format(e))\n",
    "\n",
    "    return e_log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_table():\n",
    "    # Init pd.DataFrame\n",
    "    headers = [\"cell_id\", \"total\", \"lat\", \"lon\", \"datasets\"]\n",
    "    months = get_months(timerange)\n",
    "    headers.extend(months)\n",
    "    count_table = pd.DataFrame(columns=headers)\n",
    "    return count_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "\n",
    "# Iterate through the input grid\n",
    "def iterate_grid(aoi_m: list, c: list, count_table) -> list:\n",
    "    \"\"\"\n",
    "    Iterates through every feature (cell) in the AOI grid.\n",
    "\n",
    "    Parameters:\n",
    "    aoi_m: list, required\n",
    "        List of geojson feature collection (cells) that make up the entire grid.\n",
    "    c: list, required\n",
    "        List of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload, if executio fails\n",
    "    count_table: pd.DataFrame, optional\n",
    "        Water pixel count table.\n",
    "\n",
    "    Returns:\n",
    "    e_log: list\n",
    "        Error log list of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload.\n",
    "    \"\"\"\n",
    "\n",
    "    e_log = []\n",
    "    cell = 1\n",
    "    months = months = get_months(timerange)\n",
    "\n",
    "    # Run the main iterator\n",
    "    for aoi, i, t in zip(\n",
    "        aoi_m, c, tqdm(range(len(c)), initial=1, desc=\"Grid Processed\")\n",
    "    ):\n",
    "        geopolygon = Geometry(aoi[\"features\"][0][\"geometry\"], crs=\"epsg:4326\")\n",
    "        geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)\n",
    "        g = geopolygon_gdf.centroid\n",
    "        print(\n",
    "            \"\\n\\n\"\n",
    "            + \"\\033[32m\"\n",
    "            + \"PROCESSING GRID CELL ID {} NO. {}/{} CENTROID ({}, {})\".format(\n",
    "                i[2], cell, len(aoi_m), round(g.y[0], 5), round(g.x[0], 5)\n",
    "            )\n",
    "            + \"\\033[0m\"\n",
    "        )\n",
    "\n",
    "        # Get the latitude and longitude range of the geopolygon\n",
    "        lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])\n",
    "        lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])\n",
    "\n",
    "        # Load Sentinel1 data\n",
    "        try:\n",
    "            S1 = load_ard(\n",
    "                dc=dc,\n",
    "                products=[\"s1_rtc\"],\n",
    "                # measurements=[\"vv\", \"vh\"],\n",
    "                measurements=[\"vh\"],\n",
    "                y=lat_range,\n",
    "                x=lon_range,\n",
    "                time=timerange,\n",
    "                output_crs=\"EPSG:6933\",\n",
    "                resolution=(-20, 20),\n",
    "                group_by=\"solar_day\",\n",
    "                dtype=\"native\",\n",
    "            )\n",
    "        except Exception as e:\n",
    "            # Log error aoi centroids and keep looping\n",
    "            e_log.append([g.x[0], g.y[0], i[2], \"P\"])\n",
    "            print(\n",
    "                \"\\n\\n\"\n",
    "                + \"\\033[31m\"\n",
    "                + \"ERROR PROCESSING GRID CELL {}/{} CENTROID ({}, {}). LOGGED CENTROID INFO in e_log\".format(\n",
    "                    i[2], len(aoi_m), round(g.y[0], 5), round(g.x[0], 5)\n",
    "                )\n",
    "                + \"\\033[0m\"\n",
    "            )\n",
    "            print(\"PROCESS ERROR: {}\".format(e))\n",
    "            cell += 1\n",
    "            continue\n",
    "\n",
    "        datasets = S1.time\n",
    "        dn = len(datasets.time)\n",
    "\n",
    "        # Initialize row for the count_table\n",
    "        count_row = [i[2], len(aoi_m), round(g.y[0], 5), round(g.x[0], 5), dn]\n",
    "\n",
    "        # timesteps = [2, 4, 6, 9, 11]\n",
    "\n",
    "        # The lee filter above doesn't handle null values\n",
    "        # We therefore set null values to 0 before applying the filter\n",
    "        valid = np.isfinite(S1)\n",
    "        S1 = S1.where(valid, 0)\n",
    "\n",
    "        # Create a new entry in dataset corresponding to filtered VV and VH data\n",
    "        S1[\"filtered_vh\"] = S1.vh.groupby(\"time\").apply(lee_filter, size=7)\n",
    "\n",
    "        # Null pixels should remain null\n",
    "        S1[\"filtered_vh\"] = S1.filtered_vh.where(valid.vh)\n",
    "\n",
    "        # Convert the digital numbers to dB\n",
    "        S1[\"filtered_vh\"] = 10 * np.log10(S1.filtered_vh)\n",
    "\n",
    "        threshold_vh = th_aoi\n",
    "\n",
    "        S1[\"water\"] = S1_water_classifier(S1.filtered_vh, threshold_vh).s1_water\n",
    "        S1Water = S1.water\n",
    "        S1_BIN = S1Water.where(S1Water > 0)\n",
    "        FS1 = S1_BIN\n",
    "        PRFS1 = S1_BIN\n",
    "\n",
    "        # Generate CSV\n",
    "        print(\"Generating row data...\")\n",
    "        for month in months:\n",
    "            S1_month = S1_BIN.sel(time=[month], method=\"nearest\").mean(dim=\"time\")\n",
    "            # Get count of all pixels that are water and not np.nan\n",
    "            count = np.count_nonzero(~np.isnan(S1_month.values))\n",
    "            count_row.append(count)\n",
    "        print(count_row)\n",
    "        print(\"Adding to table\")\n",
    "        count_table = pd.concat(\n",
    "            [pd.DataFrame([count_row], columns=count_table.columns), count_table],\n",
    "            ignore_index=True,\n",
    "        )\n",
    "        count_table.to_csv(\n",
    "            \"output/csv/{}.csv\".format(process_name), index=False, mode=\"w\"\n",
    "        )\n",
    "\n",
    "        # Creating outputs\n",
    "        # Export to raster - upload to g-drive - delete from sandbox\n",
    "        print(\"Analysis periods:\")\n",
    "        for p in a_periods:\n",
    "            print(p)\n",
    "\n",
    "        print(\"\\nUploading...\")\n",
    "        # -------------- maxflood ----------------\n",
    "        if \"max_flood\" in a_periods:\n",
    "            if i[3] in [None, \"P\", \"U-PRF\"]:\n",
    "                S1_MX = S1_BIN.sel(time=max_flood, method=\"nearest\").max(dim=\"time\")\n",
    "                err = gen_output(S1_MX, i, geopolygon, cell, aoi_m, \"maxflood\", \"max\")\n",
    "\n",
    "        # # -------------- preflood ----------------\n",
    "        if \"pre_flood\" in a_periods:\n",
    "            if i[3] in [None, \"P\", \"U-PRF\"]:\n",
    "                S1_PRF_MD = PRFS1.sel(time=pre_flood, method=\"nearest\").median(\n",
    "                    dim=\"time\"\n",
    "                )\n",
    "                S1_PRF_MN = PRFS1.sel(time=pre_flood, method=\"nearest\").mean(dim=\"time\")\n",
    "                if \"median\" in a_measures:\n",
    "                    err = gen_output(\n",
    "                        S1_PRF_MD, i, geopolygon, cell, aoi_m, \"preflood\", \"median\"\n",
    "                    )\n",
    "                    if len(err) > 0:\n",
    "                        e_log.extend(err)\n",
    "                if \"mean\" in a_measures:\n",
    "                    err = gen_output(\n",
    "                        S1_PRF_MN, i, geopolygon, cell, aoi_m, \"preflood\", \"mean\"\n",
    "                    )\n",
    "                    if len(err) > 0:\n",
    "                        e_log.extend(err)\n",
    "\n",
    "        # # --------------- flood ------------------\n",
    "        if \"flood\" in a_periods:\n",
    "            if i[3] in [None, \"P\", \"U-F\"]:\n",
    "                S1_F_MD = FS1.sel(time=flood, method=\"nearest\").median(dim=\"time\")\n",
    "                S1_F_MN = FS1.sel(time=flood, method=\"nearest\").mean(dim=\"time\")\n",
    "                if \"median\" in a_measures:\n",
    "                    err = gen_output(\n",
    "                        S1_F_MD, i, geopolygon, cell, aoi_m, \"flood\", \"median\"\n",
    "                    )\n",
    "                    if len(err) > 0:\n",
    "                        e_log.extend(err)\n",
    "                if \"mean\" in a_measures:\n",
    "                    err = gen_output(\n",
    "                        S1_F_MN, i, geopolygon, cell, aoi_m, \"flood\", \"mean\"\n",
    "                    )\n",
    "                    if len(err) > 0:\n",
    "                        e_log.extend(err)\n",
    "\n",
    "        # # ------------ flood-extents --------------\n",
    "        if \"flood_extents\" in a_periods:\n",
    "            if i[3] in [None, \"P\", \"U-FE\"]:\n",
    "                S1_FE_MD = S1_F_MD - S1_PRF_MD\n",
    "                S1_FE_MN = S1_F_MN - S1_PRF_MN\n",
    "                if \"median\" in a_measures:\n",
    "                    err = gen_output(\n",
    "                        S1_FE_MD, i, geopolygon, cell, aoi_m, \"flood_extents\", \"median\"\n",
    "                    )\n",
    "                    if len(err) > 0:\n",
    "                        e_log.extend(err)\n",
    "                if \"mean\" in a_measures:\n",
    "                    err = gen_output(\n",
    "                        S1_FE_MN, i, geopolygon, cell, aoi_m, \"flood_extents\", \"mean\"\n",
    "                    )\n",
    "                    if len(err) > 0:\n",
    "                        e_log.extend(err)\n",
    "\n",
    "        cell += 1\n",
    "        clear_output()\n",
    "\n",
    "    if len(e_log) == 0:\n",
    "        print(\n",
    "            \"\\n\\n\"\n",
    "            + \"\\033[32m\"\n",
    "            + \"GRID PROCESSED AND UPLOADED SUCCESSFULLY\"\n",
    "            + \"\\033[0m\"\n",
    "            + \"\\n\\n\"\n",
    "        )\n",
    "    else:\n",
    "        e_log = gen_elog(e_log)\n",
    "\n",
    "    # return e_log to be run again\n",
    "    return e_log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Crete the aoi-mosaic - aoi_m\n",
    "def gen_aoim(c: list, b: float, count_table) -> list:\n",
    "    \"\"\"\n",
    "    Generates the feature collection list (list of cells) using centroid coordinates and a buffer distance. Calls the main iterator for execution as well.\n",
    "\n",
    "    Parameters:\n",
    "    c: list, required\n",
    "        List of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload, if executio fails.\n",
    "    b: float, required\n",
    "        Cell half-dimension in degrees (EPSG:4326). Creates a cell by adding this distance to the centroid coordinates.\n",
    "    count_table: pd.DataFrame, optional\n",
    "        Water pixel count table.\n",
    "\n",
    "    Returns:\n",
    "    e_log: list\n",
    "        Error log list of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload.\n",
    "    \"\"\"\n",
    "    aoi_m = []\n",
    "    for i in c:\n",
    "        aoi_m.append(define_area(i[1], i[0], buffer=b))\n",
    "    # print(c, len(aoi_m))\n",
    "    e_log = iterate_grid(aoi_m, c, count_table)\n",
    "\n",
    "    # return e_log to be run again\n",
    "    return e_log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize input file\n",
    "def view_input(gdf_list: list[gpd.GeoDataFrame], grid_c: list) -> None:\n",
    "    \"\"\"\n",
    "    Visualizes cells and respective IDs  on a basemap.\n",
    "\n",
    "    Parameters:\n",
    "    gdf_list:list[gpd.GeoDataFrame], required\n",
    "        List of geodataframes to be visualized.\n",
    "    grid_c:list, required\n",
    "        List of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload, if executio fails.\n",
    "\n",
    "    Returns:\n",
    "    None\n",
    "\n",
    "    \"\"\"\n",
    "    print(\"Visualizing data...\")\n",
    "    p = gdf_list[0].dissolve()\n",
    "    center = p.centroid\n",
    "    map = folium.Map(location=[center.y, center.x], tiles=\"CartoDB Positron\")\n",
    "\n",
    "    for gdf in gdf_list:\n",
    "        folium.GeoJson(gdf, name=\"{}\".format(gdf)).add_to(map)\n",
    "\n",
    "    for c in grid_c:\n",
    "        folium.Marker(\n",
    "            location=[c[1], c[0]],\n",
    "            popup=f\"Centroid: {c[1]}, {c[0]}\",\n",
    "            icon=folium.DivIcon(\n",
    "                icon_size=(10, 10),\n",
    "                icon_anchor=(0, 0),\n",
    "                html='<div style=\"font-size: 10pt\">{}</div>'.format(c[2]),\n",
    "            ),\n",
    "        ).add_to(map)\n",
    "\n",
    "    bounds = gdf_list[0].total_bounds.tolist()\n",
    "    map.fit_bounds([bounds[:2][::-1], bounds[2:][::-1]])\n",
    "    display(map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create grid\n",
    "def create_grid(adm0: gpd.GeoDataFrame, size: float) -> gpd.GeoDataFrame:\n",
    "    \"\"\"\n",
    "    Divides adm0 AOI vectorfile into square grid based on size\n",
    "\n",
    "    Parameters:\n",
    "    adm0:gpd.GeoDataFrame, required\n",
    "        AMD0 GeoDataFrame created from ADM0 input vector file\n",
    "    size:float, required\n",
    "        Grid cell size in degrees (EPSG:4326)\n",
    "\n",
    "    Returns:\n",
    "    grid: gpd.GeoDataFrame\n",
    "        The generated grid GeoDataFrame\n",
    "    \"\"\"\n",
    "    bounds = adm0.bounds\n",
    "    minx = bounds.minx[0]  # only 1 feature at the 0th index\n",
    "    miny = bounds.miny[0]\n",
    "    maxx = bounds.maxx[0]\n",
    "    maxy = bounds.maxy[0]\n",
    "\n",
    "    grid = gpd.GeoDataFrame()\n",
    "    for x0 in np.arange(minx, maxx, size):\n",
    "        for y0 in np.arange(miny, maxy, size):\n",
    "            x1 = x0 + size\n",
    "            y1 = y0 + size\n",
    "            d = {\"geometry\": [Polygon([(x0, y0), (x1, y0), (x1, y1), (x0, y1)])]}\n",
    "            cell = gpd.GeoDataFrame(d, crs=\"EPSG:4326\")\n",
    "            flag = adm0.intersection(cell)\n",
    "            if flag[0].is_empty == False:\n",
    "                grid = pd.concat([grid, cell])\n",
    "\n",
    "    return grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check CRS and convert to 4326 if required\n",
    "def crs_check(shp: gpd.GeoDataFrame) -> gpd.GeoDataFrame:\n",
    "    \"\"\"\n",
    "    Checks input GeoDataFrame CRS and converts to EPSG 4326, if different.\n",
    "\n",
    "    Parameters:\n",
    "    shp: gpd.GeoDataFrame, required\n",
    "        Input GeoDataFrame to check.\n",
    "\n",
    "    Returns:\n",
    "    shp: gpd.GeoDataFrame\n",
    "        As is or converted GeoDataFrame.\n",
    "    \"\"\"\n",
    "    if shp.crs != \"EPSG:4326\":\n",
    "        print(\"Added ADM0 CRS is {}. Converting to EPSG:4326...\".format(shp.crs))\n",
    "        shp = shp.to_crs(\"EPSG:4326\")\n",
    "        if shp.crs == \"EPSG:4326\":\n",
    "            print(\"Done\")\n",
    "\n",
    "    return shp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_inp(x: str) -> None:\n",
    "    \"\"\"\n",
    "    Checks if input is \"y\" or \"n\"\n",
    "\n",
    "    Parameters:\n",
    "    x: str, required\n",
    "        String input to be checked\n",
    "\n",
    "    Returns:\n",
    "    None\n",
    "    \"\"\"\n",
    "    if x not in [\"y\", \"n\"]:\n",
    "        raise ValueError(\"Invalid input, must be 'y' or 'n'\")\n",
    "    elif x == \"n\":\n",
    "        raise RuntimeError(\n",
    "            \"Excecution terminated. Make necessary changes before running again\"\n",
    "        )\n",
    "\n",
    "\n",
    "def exec_checks(c, buffer, count_table) -> list:\n",
    "    \"\"\"\n",
    "    Performs checks and Run the entire application\n",
    "\n",
    "    Parameters:\n",
    "    c: list, required\n",
    "        List of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload, if executio fails.\n",
    "    buffer: float, required\n",
    "        Cell half-dimension in degrees (EPSG:4326). Creates a cell by adding this distance to the centroid coordinates.\n",
    "\n",
    "    Returns:\n",
    "    Returns:\n",
    "    e_log: list\n",
    "        Error log list of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload.\n",
    "    \"\"\"\n",
    "\n",
    "    inst = \"\"\"\n",
    "Before running the execution, ensure all requirements have been met:\n",
    "1. Create appropriate folders in Google Drive and add their Folder IDs here\n",
    "2. Check input shapefile\n",
    "3. Check grid size\n",
    "    \n",
    "    \"\"\"\n",
    "    print(inst)\n",
    "\n",
    "    x = input(\"Folder IDs verified? (y/n):\")\n",
    "    check_inp(x)\n",
    "\n",
    "    x = input(\"Input shapefile/geojson verified? (y/n):\")\n",
    "    check_inp(x)\n",
    "\n",
    "    input(\"Grid (size) verified? (y/n):\")\n",
    "    check_inp(x)\n",
    "\n",
    "    z = input(\"\\nBegin execution for entire input shapefile/geojson? (y/n):\")\n",
    "    if z not in [\"y\", \"n\"]:\n",
    "        raise ValueError(\"Invalid input, must be 'y' or 'n'\")\n",
    "    elif z == \"n\":\n",
    "        raise RuntimeError(\n",
    "            \"Excecution terminated. Make necessary changes before running again\"\n",
    "        )\n",
    "    elif z == \"y\":\n",
    "        print(\"Starting execution...\")\n",
    "        # get e_log with centroids, cell_id and error message\n",
    "        # Calling gen_aoim will run the entire iteration process\n",
    "        e_log = gen_aoim(c, buffer, count_table)\n",
    "        print(len(e_log))\n",
    "        return e_log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def del_files(path: str, ext: str) -> None:\n",
    "    \"\"\"\n",
    "    Deletes all files with specified extention at specified folder path\n",
    "\n",
    "    Parameters:\n",
    "    path:str, required\n",
    "        Folder path.\n",
    "    ext:str, required\n",
    "        File extension. \"*\" for all files.\n",
    "\n",
    "    Returns:\n",
    "    None\n",
    "    \"\"\"\n",
    "    res_files = False\n",
    "    loc = os.path.join(path, ext)\n",
    "    files = glob.glob(loc)\n",
    "    if len(files) > 0:\n",
    "        res_files = True\n",
    "        print(\"Found {} files. Deleting...\".format(len(files)))\n",
    "        for f in files:\n",
    "            os.remove(f)\n",
    "\n",
    "    return res_files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clean_dirs() -> None:\n",
    "    \"\"\"\n",
    "    Creates directories if they dont exist or deletes residual files if they exist.\n",
    "\n",
    "    Parameters:\n",
    "    None\n",
    "\n",
    "    Returns:\n",
    "    None\n",
    "    \"\"\"\n",
    "    dirs_exist = False\n",
    "    dir_dict = {\n",
    "        \"sd_flood\": \"output/flood\",\n",
    "        \"sd_preflood\": \"output/preflood\",\n",
    "        \"sd_flood_extents\": \"output/flood_extents\",\n",
    "    }\n",
    "\n",
    "    for k in dir_dict:\n",
    "        if not os.path.exists(dir_dict[k]):\n",
    "            os.makedirs(dir_dict[k])\n",
    "        else:\n",
    "            dirs_exist = True\n",
    "            r = del_files(dir_dict[k], \"*\")\n",
    "\n",
    "    if dirs_exist:\n",
    "        print(\"Output folders alredy exist.\")\n",
    "    else:\n",
    "        print(\"Output folders created.\")\n",
    "\n",
    "    if not r:\n",
    "        print(\"No residual files to delete.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_shapefile(centroids, filename, preview=False):\n",
    "    # Export grid polygon\n",
    "    o_grid = gpd.GeoDataFrame()\n",
    "    for i in centroids:\n",
    "        point = Point(i[0], i[1])  # This takes x first and then y\n",
    "        gdf = gpd.GeoDataFrame(geometry=[point])\n",
    "        buffer = cell_size / 2\n",
    "        cell = gpd.GeoDataFrame()\n",
    "        cell[\"geometry\"] = gdf.buffer(buffer, cap_style=\"square\")\n",
    "        cell.set_geometry(\"geometry\")\n",
    "        cell[\"cell_id\"] = i[2]\n",
    "        o_grid = pd.concat([o_grid, cell])\n",
    "\n",
    "    o_grid = o_grid.set_crs(\"epsg:4326\")\n",
    "    o_grid.to_file(\"output/shape/{}_grid.geojson\".format(filename), driver=\"GeoJSON\")\n",
    "    if preview:\n",
    "        view_input([o_grid], centroids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_adm(\n",
    "    shp: str,\n",
    "    boundary_buffer: float,\n",
    "    cell_size: float,\n",
    "    test_cells: list = None,\n",
    "    preview: bool = False,\n",
    ") -> list | gpd.GeoDataFrame | pd.DataFrame:\n",
    "    \"\"\"\n",
    "    Processes the input ADM file. File does not have to be ADM0 and can be any vectorfile.\n",
    "\n",
    "    Parameters:\n",
    "    shp: str, required\n",
    "        Path of input vectorfile. Will be converted to EPSG:4326.\n",
    "    boundary_buffer: float, required\n",
    "        Outer boundary buffer to be given to the input vector file in degrees (EPSG:4326)\n",
    "    cell_size: float, required\n",
    "        Grid cell size in degrees (EPSG:4326)\n",
    "    preview: bool, required\n",
    "        Show map preview, False by default\n",
    "\n",
    "    Returns:\n",
    "    c: list\n",
    "        List of [x, y, cell_id and None]. None will store the error \"P\" - Processing or \"U\" - Upload, if executio fails.\n",
    "    grid: GeoDataFrame\n",
    "        Grid file generated from input file\n",
    "    adm_df: Dataframe\n",
    "        Input file information\n",
    "    \"\"\"\n",
    "    adm0_b = gpd.read_file(shp)  # adm0 base\n",
    "    adm0_b = adm0_b.dissolve()\n",
    "    adm0_buf = adm0_b.buffer(boundary_buffer)  # adm0 with 20KM boundary buffer\n",
    "    adm0 = crs_check(adm0_buf)\n",
    "    size = cell_size  # Grid cell size 0.5 ~ 55KM\n",
    "    buffer = size / 2  # cell buffer around the centroid to create the cell\n",
    "\n",
    "    grid = create_grid(adm0, size)\n",
    "    # Calculate centroids and store in centroid list c[].\n",
    "    c = []\n",
    "    g = grid.centroid\n",
    "\n",
    "    cell_id = 1\n",
    "    for i in g:\n",
    "        c.append(\n",
    "            [round(i.x, 5), round(i.y, 5), cell_id, None]\n",
    "        )  # The array c[] has four values: x, y, cell_id and None. None will store the \"P\" or \"U\" error value\n",
    "        cell_id += 1\n",
    "\n",
    "    n = len(c)\n",
    "\n",
    "    if test_cells != None:\n",
    "        if len(test_cells) == 2:\n",
    "            c = c[test_cells[0] : test_cells[1]]\n",
    "            to_process = test_cells[1] - test_cells[0]\n",
    "        elif len(test_cells) == 0:\n",
    "            to_process = 0\n",
    "\n",
    "    if preview == True:\n",
    "        view_input([grid, adm0], c)\n",
    "\n",
    "    adm_data = {\n",
    "        \"Parameter\": [\n",
    "            \"Input File Path\",\n",
    "            \"Area\",\n",
    "            \"Area with Buffer\",\n",
    "            \"Cell Size\",\n",
    "            \"Total Cells\",\n",
    "            \"Test Cells\",\n",
    "        ],\n",
    "        \"Value\": [\n",
    "            shp,\n",
    "            \"{} KM2\".format(\n",
    "                round((adm0_b.to_crs(\"EPSG:3857\").area).iloc[0] / (10**6), 2)\n",
    "            ),\n",
    "            \"{} KM2\".format(\n",
    "                round((adm0.to_crs(\"EPSG:3857\").area).iloc[0] / (10**6), 2)\n",
    "            ),\n",
    "            \"{} deg\".format(cell_size),\n",
    "            n,\n",
    "            to_process,\n",
    "        ],\n",
    "    }\n",
    "\n",
    "    adm_df = pd.DataFrame(adm_data)\n",
    "    adm_df.style.set_caption(\"Input Data\")\n",
    "\n",
    "    return (c, grid, adm_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outputs\n",
    "## [Jump to config](#Config-Cell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load file from sandbox disc. file should be present in 'input' folder\n",
    "c, grid, adm_df = add_adm(\n",
    "    shp, buffer, cell_size, test_cells, preview=False\n",
    ")  # (shp, boundary_buffer, cell_size)\n",
    "adm_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check and Run Application"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# Init Table\n",
    "count_table = gen_table()\n",
    "\n",
    "# Calls the checklist function\n",
    "if run_checks:\n",
    "    e_log = exec_checks(c, cell_size / 2, count_table)\n",
    "else:\n",
    "    e_log = gen_aoim(c, cell_size / 2, count_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tools.misc import get_maxmonth\n",
    "\n",
    "get_maxmonth(\"output/csv/{}.csv\".format(process_name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e_log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create processed and error shapefiles as grids in output/shape/\n",
    "create_shapefile(c, \"{}_BASE\".format(process_name))\n",
    "if len(e_log) > 0:\n",
    "    create_shapefile(e_log, \"{}_ERR\".format(process_name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-run application for cells logged in e_log, if required\n",
    "if len(e_log) > 0:\n",
    "    e_log = gen_aoim(e_log, cell_size / 2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}