{ "cells": [ { "cell_type": "markdown", "id": "professional-display", "metadata": {}, "source": [ "## Testing multiple points of failure\n", "\n", "This notebook forms the basis of \"Hands-On 7\" in the CCG course.\n", "\n", "1. Show asset exposure based on flood return period maps, select all exposed assets in an area\n", " - note that if historic flood outlines are available, these could be used to select multiple assets instead\n", "2. Demonstrate how the network effects of multiple asset failures would be disproportionately worse if there is no redundancy or potential for substitution. \n", "3. Understand the risk of multiple failures across the network - target the combinations of failures of greatest consequence\n", "\n", "By the end of this tutorial you should be able to:\n", "* Assess direct and some indirect impacts of multiple asset failures\n", "* Compare flooding within regions as source of multiple failure scenarios\n", "* Understand approaches to stress-testing the system under multiple failures" ] }, { "cell_type": "code", "execution_count": null, "id": "rolled-wireless", "metadata": {}, "outputs": [], "source": [ "# Imports from Python standard library\n", "import os\n", "import warnings\n", "from glob import glob\n", "from math import factorial\n", "from pathlib import Path\n", "\n", "# Imports from other Python packages\n", "import contextily as cx\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 pyproj import Geod\n", "from tqdm.notebook import tqdm" ] }, { "cell_type": "markdown", "id": "ordinary-council", "metadata": {}, "source": [ "Change this to point to your data folder as in the previous tutorial:" ] }, { "cell_type": "code", "execution_count": null, "id": "international-saint", "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\")" ] }, { "cell_type": "code", "execution_count": null, "id": "53a0c038", "metadata": {}, "outputs": [], "source": [ "# TODO add a note / auto download the geoboundaries dataset\n", "# https://data.humdata.org/dataset/geoboundaries-admin-boundaries-for-ghana/resource/e324acb9-fd6b-4f5b-96ee-e2bb11917942" ] }, { "cell_type": "markdown", "id": "6ad32f28", "metadata": {}, "source": [ "### 0.5 Preliminary step\n", "\n", "Go to the geoboundatries at: https://www.geoboundaries.org/countryDownloads.html\n", "\n", "Search for *Ghana* in the **Name** field and download the data set from **Source** Open Street Map. It will show up as a folder in your downloads directory. \n", "\n", "Put this folder in the data directory of this tutorial. The one you specify in the `data_folder` variable. # TODO add these steps to do this automatically - https://www.geoboundaries.org/api.html#api" ] }, { "cell_type": "code", "execution_count": null, "id": "bc9bf1c8", "metadata": {}, "outputs": [], "source": [ "# https://www.geoboundaries.org/api/current/gbOpen/GHA/ADM1/\n", "\n", "import requests\n", "import json\n", "from urllib.request import urlretrieve\n", "import zipfile\n", "\n", "url = \"https://www.geoboundaries.org/api/current/gbOpen/GHA/ADM1/\"\n", "response = requests.post(url)\n", "test = json.loads(response.content)\n", "download_url = test[\"staticDownloadLink\"]\n", "file_path = os.path.join(data_folder, \"geoBoundaries-GHA-ADM1-all.zip\")\n", "if not os.path.exists(Path(os.path.splitext(file_path)[0])):\n", " # check if data already exists in this directory\n", " # if len(os.listdir(Path(os.path.splitext(file_path)[0]))) == 0:\n", " urlretrieve(download_url, file_path) # note: this can take a few minutes\n", " with zipfile.ZipFile(file_path, \"r\") as zip_ref:\n", " zip_ref.extractall(os.path.splitext(file_path)[0])\n", "else:\n", " print(\"data already exists\")" ] }, { "cell_type": "markdown", "id": "optimum-clear", "metadata": {}, "source": [ "### 1. Map exposure" ] }, { "cell_type": "code", "execution_count": null, "id": "middle-mercy", "metadata": {}, "outputs": [], "source": [ "roads = gpd.read_file(data_folder + \"/GHA_OSM_roads.gpkg\", layer=\"edges\")\n", "road_nodes = gpd.read_file(data_folder + \"/GHA_OSM_roads.gpkg\", layer=\"nodes\")" ] }, { "cell_type": "code", "execution_count": null, "id": "structured-hurricane", "metadata": {}, "outputs": [], "source": [ "regions = gpd.read_file(\n", " os.path.join(\n", " data_folder,\n", " \"geoBoundaries-GHA-ADM1-all\",\n", " \"geoBoundaries-GHA-ADM1.shp\",\n", " )\n", ")[[\"shapeName\", \"shapeISO\", \"geometry\"]]" ] }, { "cell_type": "code", "execution_count": null, "id": "20abd8e0", "metadata": {}, "outputs": [], "source": [ "regions.head(5)" ] }, { "cell_type": "code", "execution_count": null, "id": "worse-million", "metadata": {}, "outputs": [], "source": [ "roads = gpd.sjoin(roads, regions).drop(columns=\"index_right\")" ] }, { "cell_type": "markdown", "id": "amateur-hometown", "metadata": {}, "source": [ "Filter roads by region name to find all roads in Greater Accra:" ] }, { "cell_type": "code", "execution_count": null, "id": "alone-diameter", "metadata": {}, "outputs": [], "source": [ "accra_roads = roads[roads.shapeName == \"Greater Accra Region\"]\n", "accra_roads.head(5)" ] }, { "cell_type": "code", "execution_count": null, "id": "legitimate-metallic", "metadata": {}, "outputs": [], "source": [ "exposure = gpd.read_parquet(\n", " data_folder + \"/results\" + \"/GHA_OSM_roads_edges___exposure.geoparquet\"\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "raised-lafayette", "metadata": {}, "outputs": [], "source": [ "exposure = gpd.sjoin(exposure, regions).drop(columns=\"index_right\")" ] }, { "cell_type": "code", "execution_count": null, "id": "16b6fc9b", "metadata": {}, "outputs": [], "source": [ "exposure.head(5)" ] }, { "cell_type": "markdown", "id": "usual-vacuum", "metadata": {}, "source": [ "Filter exposure by region, RCP and return period to find all roads exposed to a historical 100-year flood in Greater Accra:" ] }, { "cell_type": "code", "execution_count": null, "id": "presidential-writing", "metadata": {}, "outputs": [], "source": [ "accra_exposure = exposure[(exposure.shapeName == \"Greater Accra Region\")]" ] }, { "cell_type": "code", "execution_count": null, "id": "3ca9f390", "metadata": {}, "outputs": [], "source": [ "accra_exposure.head(5)" ] }, { "cell_type": "code", "execution_count": null, "id": "5e78168b", "metadata": {}, "outputs": [], "source": [ "[col for col in accra_exposure.columns if \"1980\" in col]" ] }, { "cell_type": "code", "execution_count": null, "id": "proof-prague", "metadata": {}, "outputs": [], "source": [ "accra_exposure.plot(\n", " column=\"wri_aqueduct-version_2-inunriver_historical_000000000WATCH_1980_rp00100-gha\"\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "2326d851", "metadata": {}, "outputs": [], "source": [ "flood_col = \"wri_aqueduct-version_2-inunriver_historical_000000000WATCH_1980_rp00100-gha\"\n", "accra_exposure_100yr = accra_exposure[accra_exposure[flood_col] > 0.5].copy()\n", "accra_exposure_100yr[[\"id\", \"road_type\", \"name\", \"length_m\", flood_col]].head()" ] }, { "cell_type": "markdown", "id": "portuguese-welding", "metadata": {}, "source": [ "### 2. Multiple failures" ] }, { "cell_type": "markdown", "id": "joint-legislation", "metadata": {}, "source": [ "Direct damage can be summed directly, if we assume that all roads are damaged in the same event:" ] }, { "cell_type": "code", "execution_count": null, "id": "proprietary-drinking", "metadata": {}, "outputs": [], "source": [ "(\n", " \"Total direct exposure, \"\n", " \"in Accra under a historical 100-year flood, is estimated to be \"\n", " f\"{int(accra_exposure_100yr.length_m.sum() // 1e3)}km (of {int(accra_exposure.length_m.sum() // 1e3)}km total roads).\"\n", ")" ] }, { "cell_type": "markdown", "id": "durable-distributor", "metadata": {}, "source": [ "Indirect damage can be assessed in different ways, some beyond the scope of this notebook. In this section, we look at the effects of disruption on a single route across the Greater Accra region. In a fuller analysis, we could extend this to look at many trips made within the region, and calculate the number of passengers or value of freight disrupted, along with the effects on transport time and cost.\n", "\n", "Start by creating a networkx graph from the roads, using `from_id`, `to_id` and `length_m`:" ] }, { "cell_type": "code", "execution_count": null, "id": "superb-american", "metadata": {}, "outputs": [], "source": [ "G = nx.Graph()\n", "G.add_edges_from(\n", " (r.from_id, r.to_id, {\"id\": r.id, \"weight\": r.length_m})\n", " for r in roads.itertuples()\n", ")" ] }, { "cell_type": "markdown", "id": "round-difference", "metadata": {}, "source": [ "Then find the shortest path from one node to another:" ] }, { "cell_type": "code", "execution_count": null, "id": "pleasant-saturn", "metadata": {}, "outputs": [], "source": [ "route_nodes = nx.algorithms.shortest_path(\n", " G, \"roadn_6700\", \"roadn_1011\", weight=\"weight\"\n", ")" ] }, { "cell_type": "markdown", "id": "ranging-broadway", "metadata": {}, "source": [ "Then find the edges in the shortest path, and sum over their lengths to find the length of the route:" ] }, { "cell_type": "code", "execution_count": null, "id": "sustained-irrigation", "metadata": {}, "outputs": [], "source": [ "def edge_ids_from_nodes(G, route_nodes):\n", " next_nodes = iter(route_nodes)\n", " next(next_nodes)\n", " return [G.edges[u, v][\"id\"] for u, v in zip(route_nodes, next_nodes)]\n", "\n", "\n", "route_edge_ids = edge_ids_from_nodes(G, route_nodes)" ] }, { "cell_type": "code", "execution_count": null, "id": "smoking-stockholm", "metadata": {}, "outputs": [], "source": [ "route = roads[roads.id.isin(route_edge_ids)]" ] }, { "cell_type": "code", "execution_count": null, "id": "guided-region", "metadata": {}, "outputs": [], "source": [ "f\"Best route: {round(route.length_m.sum() / 1e3, 2)}km\"" ] }, { "cell_type": "code", "execution_count": null, "id": "talented-contemporary", "metadata": {}, "outputs": [], "source": [ "ax = route.plot()\n", "ax.set_title(\"Direct route\")\n", "cx.add_basemap(ax, crs=route.crs, source=cx.providers.CartoDB.Positron)\n", "ax" ] }, { "cell_type": "markdown", "id": "mental-programming", "metadata": {}, "source": [ "Save figure to file:" ] }, { "cell_type": "code", "execution_count": null, "id": "horizontal-prerequisite", "metadata": {}, "outputs": [], "source": [ "fig = ax.get_figure()\n", "fig.savefig(os.path.join(data_folder, \"results\", \"direct_route.png\"))" ] }, { "cell_type": "markdown", "id": "individual-council", "metadata": {}, "source": [ "Define a function which runs the process we went through above:\n", "- build a graph\n", "- reweight all failed edges to have infinite cost - failed edges need to be provided as a list of `(\"from_id\", \"to_id\")` tuples.\n", "- find the shortest route from source to target\n", "- return a dataframe with the route's road segments" ] }, { "cell_type": "code", "execution_count": null, "id": "nasty-camera", "metadata": {}, "outputs": [], "source": [ "def calc_route(roads, failures, source, target):\n", " G = nx.Graph()\n", " G.add_edges_from(\n", " (r.from_id, r.to_id, {\"id\": r.id, \"weight\": r.length_m})\n", " for r in roads.itertuples()\n", " )\n", "\n", " reweight = {}\n", " for from_id, to_id in failures:\n", " reweight[(from_id, to_id)] = float(\"inf\")\n", " nx.set_edge_attributes(G, reweight, \"weight\")\n", "\n", " route_nodes = nx.algorithms.shortest_path(\n", " G, source, target, weight=\"weight\"\n", " )\n", " route_edge_ids = edge_ids_from_nodes(G, route_nodes)\n", " route = roads[roads.id.isin(route_edge_ids)]\n", "\n", " return route" ] }, { "cell_type": "code", "execution_count": null, "id": "193fff97", "metadata": {}, "outputs": [], "source": [ "def plot_route(df):\n", " ax = df.plot()\n", " cx.add_basemap(ax, crs=route.crs, source=cx.providers.CartoDB.Positron)\n", " return ax" ] }, { "cell_type": "markdown", "id": "suspended-prague", "metadata": {}, "source": [ "Test a single road failure to find if disruption makes a difference to the overall route: " ] }, { "cell_type": "code", "execution_count": null, "id": "chicken-session", "metadata": {}, "outputs": [], "source": [ "single_failures = [(\"roadn_8900\", \"roadn_9227\")]\n", "single_fail_route = calc_route(\n", " roads, single_failures, \"roadn_6700\", \"roadn_1011\"\n", ")\n", "print(f\"Best route: {round(single_fail_route.length_m.sum() / 1e3, 2)}km\")\n", "plot_route(single_fail_route)" ] }, { "cell_type": "markdown", "id": "small-cherry", "metadata": {}, "source": [ "The single road failure above has almost no effect. In our dataset, the lanes of this road are represented separately, so the routing algorithm finds a route which goes around the failed link by switching to the other lane, and the whole journey is only about 10m longer. \n", "\n", "Let's see the effect of both lanes flooded at the same time, which may be more realistic:" ] }, { "cell_type": "code", "execution_count": null, "id": "decimal-cloud", "metadata": {}, "outputs": [], "source": [ "both_lanes_failures = [\n", " (\"roadn_8900\", \"roadn_9227\"),\n", " (\"roadn_9226\", \"roadn_8899\"),\n", "]\n", "both_lanes_fail_route = calc_route(\n", " roads, both_lanes_failures, \"roadn_6700\", \"roadn_1011\"\n", ")\n", "print(f\"Best route: {round(both_lanes_fail_route.length_m.sum() / 1e3, 1)}km\")\n", "plot_route(both_lanes_fail_route)" ] }, { "cell_type": "markdown", "id": "useful-argument", "metadata": {}, "source": [ "This results in a much longer route around the flooded link.\n", "\n", "What if more than one road is disrupted at the same time? Let's test what happens if we assume that all roads exposed to 1-in-100 year flood events anywhere in Greater Accra are impassible." ] }, { "cell_type": "code", "execution_count": null, "id": "republican-payment", "metadata": {}, "outputs": [], "source": [ "multi_failures = [\n", " (road.from_id, road.to_id) for road in accra_exposure.itertuples()\n", "]\n", "multi_fail_route = calc_route(\n", " roads, multi_failures, \"roadn_6700\", \"roadn_1011\"\n", ")\n", "print(f\"Best route: {round(multi_fail_route.length_m.sum() / 1e3, 1)}km\")\n", "plot_route(multi_fail_route)" ] }, { "cell_type": "markdown", "id": "abroad-moscow", "metadata": {}, "source": [ "This gives a longer route again.\n", "\n", "This is a quick way of coming up with a hypothetical flood event, but it is not a rigorous method of analysis. With historic flood outlines, we could test and validate this simple model against past events. With an event set output from a hydrological model (rather than just the return-period hazard map that we've been using), we could test an ensemble of potential events.\n", "\n", "The next section looks at testing all possible combinations of failures, which doesn't require any additional data or modelling. " ] }, { "cell_type": "markdown", "id": "further-palace", "metadata": {}, "source": [ "### 3. Test combinations " ] }, { "cell_type": "markdown", "id": "thrown-breach", "metadata": {}, "source": [ "We can calculate the number of possible combinations of failures, and it gets very large quite quickly. \n", "\n", "For example, with three roads, {A, B, C}, there are three possible single failures ({only A}, {only B} or {only C}), three possible double failures ({A and B}, {B and C} or {A and C}), and one possible triple failure ({A, B and C}).\n", "\n", "More formally, if a set has $n$ elements, the number of ways of picking $k$ elements from it can be shown to be:\n", "\n", "$$ {\\binom {n}{k}}={\\frac {n(n-1)\\dotsb (n-k+1)}{k(k-1)\\dotsb 1}}={\\frac {n!}{k!(n-k)!}} $$\n", "\n", "and is zero when $ k > n $.\n", "\n", "The function `n_choose_k` calculates this:" ] }, { "cell_type": "code", "execution_count": null, "id": "cognitive-absolute", "metadata": {}, "outputs": [], "source": [ "def n_choose_k(n, k):\n", " if k > n:\n", " return 0\n", " return int(factorial(n) / (factorial(k) * factorial(n - k)))" ] }, { "cell_type": "markdown", "id": "macro-repeat", "metadata": {}, "source": [ "Try out a few values to see how the function behaves:" ] }, { "cell_type": "code", "execution_count": null, "id": "breeding-pointer", "metadata": {}, "outputs": [], "source": [ "n_choose_k(3, 2)" ] }, { "cell_type": "code", "execution_count": null, "id": "naughty-builder", "metadata": {}, "outputs": [], "source": [ "n_choose_k(200, 2)" ] }, { "cell_type": "code", "execution_count": null, "id": "figured-examination", "metadata": {}, "outputs": [], "source": [ "n_choose_k(200, 3)" ] }, { "cell_type": "markdown", "id": "liquid-family", "metadata": {}, "source": [ "Calculate some of the numbers of possible failure combinations within our road network:" ] }, { "cell_type": "code", "execution_count": null, "id": "compatible-surrey", "metadata": {}, "outputs": [], "source": [ "n = len(roads)\n", "print(f\"With {n} roads\")\n", "for k in range(4):\n", " print(\n", " f\"there are {n_choose_k(n, k):,} total possible combinations of {k} roads failing\"\n", " )" ] }, { "cell_type": "markdown", "id": "departmental-designation", "metadata": {}, "source": [ "Use the `np.random.choice` to sample failure combinations at random from all roads (regardless of whether they intersect with any hazard):" ] }, { "cell_type": "code", "execution_count": null, "id": "prerequisite-proof", "metadata": {}, "outputs": [], "source": [ "k = 500\n", "ids = np.random.choice(roads.id, size=k, replace=False)\n", "failed_roads = roads[roads.id.isin(ids)]\n", "failed_roads.head(2)" ] }, { "cell_type": "code", "execution_count": null, "id": "underlying-popularity", "metadata": {}, "outputs": [], "source": [ "random_failures = [\n", " (road.from_id, road.to_id) for road in failed_roads.itertuples()\n", "]\n", "random_fail_route = calc_route(\n", " roads, random_failures, \"roadn_6700\", \"roadn_1011\"\n", ")\n", "print(f\"Best route: {round(random_fail_route.length_m.sum() / 1e3, 1)}km\")\n", "plot_route(random_fail_route)" ] }, { "cell_type": "markdown", "id": "electoral-begin", "metadata": {}, "source": [ "Sample 100 different sets of 500 failures to test how the best route length (for this arbitrarily-chosen route) changes under random failure conditions:" ] }, { "cell_type": "code", "execution_count": null, "id": "compound-lender", "metadata": {}, "outputs": [], "source": [ "k = 500\n", "n_samples = 100\n", "lengths = []\n", "for _ in tqdm(range(n_samples)):\n", " ids = np.random.choice(roads.id, size=k, replace=False)\n", " failed_roads = roads[roads.id.isin(ids)]\n", " random_failures = [\n", " (road.from_id, road.to_id) for road in failed_roads.itertuples()\n", " ]\n", " random_fail_route = calc_route(\n", " roads, random_failures, \"roadn_6700\", \"roadn_1011\"\n", " )\n", " length = round(random_fail_route.length_m.sum() / 1e3, 1)\n", " lengths.append(length)" ] }, { "cell_type": "code", "execution_count": null, "id": "finished-passenger", "metadata": {}, "outputs": [], "source": [ "sampled_failures = pd.DataFrame({\"length_km\": lengths})" ] }, { "cell_type": "markdown", "id": "mounted-engineer", "metadata": {}, "source": [ "Calculate basic summary statistics from this sample:" ] }, { "cell_type": "code", "execution_count": null, "id": "engaging-conditions", "metadata": {}, "outputs": [], "source": [ "sampled_failures.describe()" ] }, { "cell_type": "markdown", "id": "statistical-hawaiian", "metadata": {}, "source": [ "Plot all the route lengths as a scatter plot, to get some visual idea of the distribution:" ] }, { "cell_type": "code", "execution_count": null, "id": "sharing-button", "metadata": {}, "outputs": [], "source": [ "sampled_failures.reset_index().plot.scatter(x=\"index\", y=\"length_km\")" ] }, { "cell_type": "markdown", "id": "exceptional-clark", "metadata": {}, "source": [ "Plot the empirical cumulative distribution function to summarise the distribution in another way. \n", "\n", "In the samples we've taken while testing this tutorial, it shows most of the time 500 random failures in the road network has little effect on the route length (around 55-60km), but some combinations of failures see a route length of up to around 180km." ] }, { "cell_type": "code", "execution_count": null, "id": "swiss-singapore", "metadata": {}, "outputs": [], "source": [ "sns.kdeplot(sampled_failures.length_km, cumulative=True)" ] } ], "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 }