{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "global-length",
   "metadata": {},
   "source": [
    "### Evaluating infrastructure adaptation options\n",
    "\n",
    "This notebook forms the basis of \"Hands-On 8\" in the CCG course.\n",
    "\n",
    "1. Take the risk results for the Ghana road damage and disruption analysis from previous hands-on sessions\n",
    "2. Assume some adaptation options - explain what this means - and show their costs\n",
    "3. Explain cost-benefit analysis (CBA) and show how to calculate Net Present Values for benefits (avoided risks) and costs\n",
    "\n",
    "By the end of this tutorial you should be able to:\n",
    "* Quantify the potential risk reduction of adaptation options\n",
    "* Prioritise assets based on cost-benefit analysis for different adaptation options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rocky-continent",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports from Python standard library\n",
    "import math\n",
    "import os\n",
    "import warnings\n",
    "from glob import glob\n",
    "from pathlib import Path\n",
    "\n",
    "# Imports from other Python packages\n",
    "import geopandas as gpd\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from tqdm.notebook import tqdm\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "found-equilibrium",
   "metadata": {},
   "source": [
    "Change this to point to your data folder as in the previous tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "exciting-portal",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir = (\n",
    "    Path(os.getcwd()).resolve().parents[3]\n",
    ")  # get parent directory of snail package\n",
    "data_folder = os.path.join(dir, \"ghana_tutorial\")\n",
    "# data_folder = Path(\"YOUR_PATH/ghana_tutorial\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "hindu-orleans",
   "metadata": {},
   "source": [
    "### 1. Load risk results"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "professional-plain",
   "metadata": {},
   "source": [
    "Read in regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "essential-translation",
   "metadata": {},
   "outputs": [],
   "source": [
    "regions = gpd.read_file(\n",
    "    data_folder + \"/geoBoundaries-GHA-ADM1-all\" + \"/geoBoundaries-GHA-ADM1.shp\"\n",
    ")[[\"shapeName\", \"shapeISO\", \"geometry\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "184e66bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, 1)\n",
    "\n",
    "regions.plot(\n",
    "    ax=ax, alpha=1, linewidth=0.5, column=\"shapeName\", edgecolor=\"black\"\n",
    ")\n",
    "\n",
    "ax.set_title(\"Admin boundaries of Ghana\")\n",
    "ax.set_xlabel(\"Longitude [deg]\")\n",
    "ax.set_ylabel(\"Latitude [deg]\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "specific-failure",
   "metadata": {},
   "source": [
    "Read in roads, join regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "wired-solomon",
   "metadata": {},
   "outputs": [],
   "source": [
    "roads = gpd.read_file(\n",
    "    data_folder + \"/GHA_OSM_roads.gpkg\", layer=\"edges\"\n",
    ").rename(columns={\"id\": \"road_id\"})\n",
    "roads = gpd.sjoin(roads, regions).drop(columns=\"index_right\")\n",
    "roads.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "surprising-vision",
   "metadata": {},
   "source": [
    "Read in risk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "economic-technical",
   "metadata": {},
   "outputs": [],
   "source": [
    "risk = pd.read_csv(data_folder + \"/results\" + \"/inunriver_damages_ead.csv\")[\n",
    "    [\"id\", \"rcp\", \"gcm\", \"epoch\", \"ead_usd\"]\n",
    "].rename(columns={\"id\": \"road_id\"})\n",
    "risk.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "composed-objective",
   "metadata": {},
   "outputs": [],
   "source": [
    "exposed_roads = roads[roads.road_id.isin(risk.road_id.unique())]\n",
    "exposed_roads.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sophisticated-pickup",
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure = pd.read_csv(data_folder + \"/results\" + \"/inunriver_damages_rp.csv\")[\n",
    "    [\"id\", \"length_m\", \"rcp\", \"gcm\", \"epoch\", \"rp\"]\n",
    "].rename(columns={\"id\": \"road_id\", \"length_m\": \"flood_length_m\"})\n",
    "\n",
    "# sum over any segments exposed within the same return period\n",
    "exposure = exposure.groupby([\"road_id\", \"rcp\", \"gcm\", \"epoch\", \"rp\"]).sum()\n",
    "\n",
    "# pick max length exposed over all return periods\n",
    "exposure = (\n",
    "    exposure.reset_index()\n",
    "    .groupby([\"road_id\", \"rcp\", \"gcm\", \"epoch\"])\n",
    "    .max()\n",
    "    .reset_index()\n",
    ")\n",
    "\n",
    "exposure.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "italian-color",
   "metadata": {},
   "outputs": [],
   "source": [
    "roads_with_risk = exposed_roads.merge(risk, on=\"road_id\", how=\"inner\").merge(\n",
    "    exposure, on=[\"road_id\", \"rcp\", \"gcm\", \"epoch\"]\n",
    ")\n",
    "roads_with_risk.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "prime-perception",
   "metadata": {},
   "source": [
    "### 2. Introduce adaptation options"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "encouraging-poverty",
   "metadata": {},
   "source": [
    "Introduce costs of road upgrade options.\n",
    "\n",
    "These costs are taken purely as an example, and further research is required to make reasonable estimates. They are intended to represent upgrade to a bituminous or concrete road design, with a single-lane design for currently-unpaved roads. The routine maintenance costs are estimated for rehabilitation and routine maintenance that should take place every year. The periodic maintenance costs are estimated for resurfacing and surface treatment that may take place approximately every five years.\n",
    "\n",
    "As before with cost estimates, the analysis is likely to be highly sensitive to these assumptions, which should be replaced by better estimates if available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "furnished-closer",
   "metadata": {},
   "outputs": [],
   "source": [
    "options = pd.DataFrame(\n",
    "    {\n",
    "        \"kind\": [\"four_lane\", \"two_lane\", \"single_lane\"],\n",
    "        \"initial_cost_usd_per_km\": [1_000_000, 500_000, 125_000],\n",
    "        \"routine_usd_per_km\": [20_000, 10_000, 5_000],\n",
    "        \"periodic_usd_per_km\": [100_000, 50_000, 25_000],\n",
    "    }\n",
    ")\n",
    "options"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "elder-prince",
   "metadata": {},
   "source": [
    "Set a discount rate. This will be used to discount the cost of annual and periodic maintenance, as well as the present value of future expected annual damages.\n",
    "\n",
    "This is another sensitive parameter which will affect the net present value calculations for both costs and benefits. As an exercise, try re-running the remainder of the analysis with different values here. What economic or financial justification could there be for assuming different discount rates?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rapid-award",
   "metadata": {},
   "outputs": [],
   "source": [
    "discount_rate_percentage = 3"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "refined-success",
   "metadata": {},
   "source": [
    "Given initial and routine costs and a discount rate, we can calculate the net present value for each adaptation option.\n",
    "\n",
    "- start by calculating the normalised discount rate for each year over the time horizon\n",
    "- add the initial costs for each option\n",
    "- calculate the discounted routine costs for each option (assumed to be incurred each year)\n",
    "- calculate the discounted periodic costs for each option (assumed to be incurred every five years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "indoor-digit",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up a costs dataframe\n",
    "costs = pd.DataFrame()\n",
    "\n",
    "# create a row per year over the time-horizon of interest\n",
    "costs[\"year\"] = np.arange(2020, 2081)\n",
    "costs[\"year_from_start\"] = costs.year - 2020\n",
    "\n",
    "# calculate the normalised discount rate\n",
    "discount_rate = 1 + discount_rate_percentage / 100\n",
    "costs[\"discount_rate_norm\"] = costs.year_from_start.apply(\n",
    "    lambda y: 1.0 / math.pow(discount_rate, y)\n",
    ")\n",
    "# calculate the sum over normalised discount rates for the time horizon\n",
    "# this will be useful later, to calculate NPV of expected damages\n",
    "discount_rate_norm = costs.discount_rate_norm.sum()\n",
    "\n",
    "# link each of the options, so we have a row per-option, per-year\n",
    "costs[\"link\"] = 1\n",
    "options[\"link\"] = 1\n",
    "costs = costs.merge(options, on=\"link\").drop(columns=\"link\")\n",
    "\n",
    "# set initial costs to zero in all years except start year\n",
    "costs.loc[costs.year_from_start > 0, \"initial_cost_usd_per_km\"] = 0\n",
    "\n",
    "# discount routine and periodic maintenance costs\n",
    "costs.routine_usd_per_km = costs.discount_rate_norm * costs.routine_usd_per_km\n",
    "costs.periodic_usd_per_km = (\n",
    "    costs.discount_rate_norm * costs.periodic_usd_per_km\n",
    ")\n",
    "# set periodic costs to zero except for every five years\n",
    "costs.loc[costs.year_from_start == 0, \"periodic_usd_per_km\"] = 0\n",
    "costs.loc[costs.year_from_start % 5 != 0, \"periodic_usd_per_km\"] = 0\n",
    "costs.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "missing-disclaimer",
   "metadata": {},
   "source": [
    "This table can then be summarised by summing over all years in the time horizon, to calculate the net present value of all that future investment in maintenance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "intermediate-mouth",
   "metadata": {},
   "outputs": [],
   "source": [
    "npv_costs = (\n",
    "    costs[\n",
    "        [\n",
    "            \"kind\",\n",
    "            \"initial_cost_usd_per_km\",\n",
    "            \"routine_usd_per_km\",\n",
    "            \"periodic_usd_per_km\",\n",
    "        ]\n",
    "    ]\n",
    "    .groupby(\"kind\")\n",
    "    .sum()\n",
    "    .reset_index()\n",
    ")\n",
    "npv_costs[\"total_cost_usd_per_km\"] = (\n",
    "    npv_costs.initial_cost_usd_per_km\n",
    "    + npv_costs.routine_usd_per_km\n",
    "    + npv_costs.periodic_usd_per_km\n",
    ")\n",
    "npv_costs"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "grand-depth",
   "metadata": {},
   "source": [
    "### 3. Estimate costs and benefits"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "honest-relation",
   "metadata": {},
   "source": [
    "Apply road kind assumptions for adaptation upgrades:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tracked-wagner",
   "metadata": {},
   "outputs": [],
   "source": [
    "def kind(road_type):\n",
    "    if road_type in (\"trunk\", \"trunk_link\", \"motorway\"):\n",
    "        return \"four_lane\"\n",
    "    elif road_type in (\"primary\", \"primary_link\", \"secondary\"):\n",
    "        return \"two_lane\"\n",
    "    else:\n",
    "        return \"single_lane\"\n",
    "\n",
    "\n",
    "roads_with_risk[\"kind\"] = roads_with_risk.road_type.apply(kind)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "compound-spine",
   "metadata": {},
   "source": [
    "Join adaptation cost estimates (per km)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "initial-independence",
   "metadata": {},
   "outputs": [],
   "source": [
    "roads_with_costs = roads_with_risk.merge(\n",
    "    npv_costs[[\"kind\", \"total_cost_usd_per_km\"]], on=\"kind\"\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "foreign-quilt",
   "metadata": {},
   "source": [
    "Calculate total cost estimate for length of roads exposed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "floppy-crowd",
   "metadata": {},
   "outputs": [],
   "source": [
    "roads_with_costs[\"total_adaptation_cost_usd\"] = (\n",
    "    roads_with_costs.total_cost_usd_per_km\n",
    "    / 1e3\n",
    "    * roads_with_costs.flood_length_m\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "simple-secretariat",
   "metadata": {},
   "source": [
    "Calculate net present value of avoided damages over the time horizon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "judicial-rehabilitation",
   "metadata": {},
   "outputs": [],
   "source": [
    "roads_with_costs[\"total_adaptation_benefit_usd\"] = (\n",
    "    roads_with_costs.ead_usd * discount_rate_norm\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "preliminary-plenty",
   "metadata": {},
   "outputs": [],
   "source": [
    "discount_rate_norm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "criminal-burning",
   "metadata": {},
   "source": [
    "Calculate benefit-cost ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accepted-charger",
   "metadata": {},
   "outputs": [],
   "source": [
    "roads_with_costs[\"bcr\"] = (\n",
    "    roads_with_costs.total_adaptation_benefit_usd\n",
    "    / roads_with_costs.total_adaptation_cost_usd\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "noble-newark",
   "metadata": {},
   "source": [
    "Filter to pull out just the historical climate scenario:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "affecting-piano",
   "metadata": {},
   "outputs": [],
   "source": [
    "historical = roads_with_costs[roads_with_costs.rcp == \"historical\"]\n",
    "historical.describe()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "rubber-destruction",
   "metadata": {},
   "source": [
    "Filter to find cost-beneficial adaptation options under historic flood scenarios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "treated-average",
   "metadata": {},
   "outputs": [],
   "source": [
    "candidates = historical[historical.bcr > 1]\n",
    "candidates.head(5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "caring-trader",
   "metadata": {},
   "source": [
    "Summarise by region to explore where cost-beneficial adaptation options might be located.\n",
    "\n",
    "We need to sum over exposed lengths of road, costs and benefits, while finding the mean benefit-cost ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "negative-liquid",
   "metadata": {},
   "outputs": [],
   "source": [
    "candidates.groupby(\"shapeName\").agg(\n",
    "    {\n",
    "        \"flood_length_m\": \"sum\",\n",
    "        \"total_adaptation_benefit_usd\": \"sum\",\n",
    "        \"total_adaptation_cost_usd\": \"sum\",\n",
    "        \"bcr\": \"mean\",\n",
    "    }\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "helpful-negative",
   "metadata": {},
   "source": [
    "Given the aggregation, filtering and plotting you've seen throughout these tutorials, what other statistics would be interesting to explore from these results?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "357dfd1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional: delete .zip files downloaded in the process\n",
    "# for item in os.listdir(data_folder):\n",
    "#     if item.endswith(\".zip\"):\n",
    "#         os.remove(os.path.join(data_folder, item))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}