Multi-country Analysis for Transportation Infrastructure#
In this notebook, we will perform a damage and risk assessment for transportation infrastructure, specifically focusing on roads across multiple countries. The assessment is based on combining hazard data (e.g., flood depths) with vulnerability curves to estimate the potential damage to infrastructure.
We will follow the steps outlined below to conduct the assessment:
Loading the necessary packages:
We will import the Python libraries required for data handling, analysis, and visualization.Loading the data:
The infrastructure data (e.g., roads) and hazard data (e.g., flood depths) will be loaded into the notebook.Preparing the data:
The infrastructure and hazard data will be processed and data gaps can be filled, if required.Performing the damage assessment:
We will apply vulnerability curves to estimate the potential damage based on the intensity of the hazard.Visualizing the results:
Finally, we will visualize the estimated damage using graphs and maps to illustrate the impact on transportation infrastructure.
1. Loading the Necessary Packages#
To perform the assessment, we are going to make use of several python packages.
In case you run this in Google Colab, you will have to install the packages below (remove the hashtag in front of them).
#!pip install damagescanner==0.9b14
#!pip install contextily
#!pip install exactextract
#!pip install lonboard
In this step, we will import all the required Python libraries for data manipulation, spatial analysis, and visualization.
import warnings
import xarray as xr
import numpy as np
import pandas as pd
import geopandas as gpd
import seaborn as sns
import shapely
from tqdm import tqdm
import matplotlib.pyplot as plt
import contextily as cx
import damagescanner.download as download
from damagescanner.core import DamageScanner
from damagescanner.osm import read_osm_data
from damagescanner.base_values import DICT_CIS_VULNERABILITY_FLOOD
from statistics import mode
from lonboard import viz
warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.simplefilter(action='ignore', category=RuntimeWarning) # exactextract gives a warning that is invalid
Specify the countries of interest#
Before we continue, we should specify the countries for which we want to assess the damage. We use the ISO3 code for each country to download the OpenStreetMap data.
countries_full_names = ['Togo','Benin','Ghana']
countries_iso3 = ['TGO','BEN','GHA']
2. Loading the Data#
In this step, we will load four key datasets:
Infrastructure data:
This dataset contains information on the location and type of transportation infrastructure (e.g., roads). Each asset may have attributes such as type, length, and replacement cost.Hazard data:
This dataset includes information on the hazard affecting the infrastructure (e.g., flood depth at various locations).Vulnerability curves:
Vulnerability curves define how the infrastructure’s damage increases with the intensity of the hazard. For example, flood depth-damage curves specify how much damage occurs for different flood depths.Maximum damage values:
This dataset contains the maximum possible damage (in monetary terms) that can be incurred by individual infrastructure assets. This is crucial for calculating total damage based on the vulnerability curves.
Infrastructure Data#
We will perform this example analysis for three countries. To start the analysis, we first download the OpenStreetMap data from GeoFabrik, open it, and store it in a dictionary.
%%time
features_per_country = {}
for country_iso3 in countries_iso3:
infrastructure_path = download.get_country_geofabrik(country_iso3)
features_per_country[country_iso3] = read_osm_data(infrastructure_path,asset_type='main_roads')
CPU times: total: 30.1 s
Wall time: 40.1 s
Now define a new set of sub_types, to be used in the remainder of this analysis
sub_types = features_per_country[country_iso3].object_type.unique()
Hazard Data#
For this example, we make use of the flood data provided by CDRI.
We use a 1/100 flood map to showcase the approach.
hazard_map = xr.open_dataset("https://hazards-data.unepgrid.ch/global_pc_h100glob.tif", engine="rasterio")
hazard_map
<xarray.Dataset> Size: 747GB
Dimensions: (band: 1, x: 432000, y: 216000)
Coordinates:
* band (band) int32 4B 1
* x (x) float64 3MB -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
* y (y) float64 2MB 90.0 90.0 90.0 90.0 ... -90.0 -90.0 -90.0 -90.0
spatial_ref int32 4B ...
Data variables:
band_data (band, y, x) float64 746GB ...Maximum damages#
One of the most difficult parts of the assessment is to define the maximum reconstruction cost of particular assets. Here we provide a baseline set of values, but these should be updated through national consultations.
maxdam_dict = {'trunk':2000,
'motorway' : 2000,
'primary':2000,
'secondary':1300,
'tertiary':700,
'trunk_link': 1300,
'motorway_link' : 1300,
'primary_link' : 1300,
'secondary_link' : 700,
'tertiary_link' : 700
}
To be used in our damage assessment, we convert this to a Pandas DataFrame
maxdam = pd.DataFrame.from_dict(maxdam_dict,orient='index').reset_index()
maxdam.columns = ['object_type','damage']
And check if any of the objects are missing from the dataframe.
missing = set(sub_types) - set(maxdam.object_type)
if len(missing) > 0:
print(f"Missing object types in maxdam: \033[1m{', '.join(missing)}\033[0m. Please add them before you continue.")
Vulnerability data#
Similarly to the maximum damages, specifying the vulnerability curves is complex. We generally have limited information about the quality of the assets, its level of deteriation and other asset-level characteristics. The study by Nirandjan et al. (2024) provides us with a baseline set of fragility and vulnerability curves that one can use. In the following cell, we load that data.
vulnerability_path = "https://zenodo.org/records/10203846/files/Table_D2_Multi-Hazard_Fragility_and_Vulnerability_Curves_V1.0.0.xlsx?download=1"
vul_df = pd.read_excel(vulnerability_path,sheet_name='F_Vuln_Depth')
And let’s have a look at all the available options
with pd.option_context('display.max_rows', None, 'display.max_colwidth', None):
display(vul_df.iloc[:2,:].T)
| 0 | 1 | |
|---|---|---|
| ID number | Infrastructure description | Additional characteristics |
| F1.1 | plant | Small power plants, capacity <100 MW |
| F1.2 | plant | Medium power plants, capacity 100-500 MW |
| F1.3 | plant | Large power plants, >500 MW |
| F1.4 | plant | thermal plant |
| F1.5 | plant | wind turbine |
| F1.6 | plant | wind turbine |
| F1.7 | plant | wind turbine |
| F2.1 | substation | Low Voltage Substation |
| F2.2 | substation | Medium Voltage Substation |
| F2.3 | substation | High Voltage Substation |
| F5.1 | cable | Distribution circuits buried crossings |
| F6.1 | Power (minor) line | Distribution circuits (non-crossings) |
| F6.2 | cable | Distribution circuits elevated crossings |
| F6.3 | Energy system | Generalized curve for energy assets in diked a... |
| F7.1 | Roads | NaN |
| F7.2 | Roads | NaN |
| F7.2a (lower boundary) | Roads | NaN |
| F7.2b (upper boundary) | Roads | NaN |
| F7.3 | Roads | NaN |
| F7.4 | Roads | Motorways and trunk roads - sophisticated acce... |
| F7.5 | Roads | Motorways and trunk roads - sophisticated acce... |
| F7.6 | Roads | Motorways and trunk roads - without sophistica... |
| F7.7 | Roads | Motorways and trunk roads - without sophistica... |
| F7.8 | Roads | Other roads - low flow |
| F7.9 | Roads | Other roads - high flow |
| F7.10 | Roads | NaN |
| F7.11 | Roads | NaN |
| F7.12 | Roads | NaN |
| F7.13 | Roads | NaN |
| F8.1 | Railways | Double-tracked railway |
| F8.2 | Railways | Double-tracked railway |
| F8.3 | Railways | NaN |
| F8.4 | Railways | NaN |
| F8.5 | Railways | NaN |
| F8.6 | Railways | NaN |
| F8.6a (lower boundary) | Railways | NaN |
| F8.6b (upper boundary) | Railways | NaN |
| F8.7 | Railways | NaN |
| F8.8 | Railways | Railway station |
| F9.1 | Airports | NaN |
| F9.2 | Airports | NaN |
| F9.3 | Airports | NaN |
| F10.1 | Telecommunication | Communication tower |
| F12.1 | Telecommunication | Communication system |
| F13.1 | water storage tanks | Water storage tanks at grade concrete |
| F13.2 | water storage tanks | Water storage tanks at grade steel |
| F13.3 | water storage tanks | Water storage tanks at grade wood |
| F13.4 | water storage tanks | Water storage tanks elevated |
| F13.5 | water storage tanks | Water storage tanks below grade |
| F14.1 | Water treatment plants | Small water treatment plants open/gravity - av... |
| F14.2 | Water treatment plants | Medium water treatment plants open/gravity - a... |
| F14.3 | Water treatment plants | Large water treatment plants open/gravity - av... |
| F14.4 | Water treatment plants | small/medium/large water treatment plants open... |
| F14.5 | Water treatment plants | small/medium/large water treatment plants open... |
| F14.6 | Water treatment plants | Small water treatment plants closed/pressure |
| F14.7 | Water treatment plants | Medium water treatment plants closed/pressure |
| F14.8 | Water treatment plants | Large water treatment plants closed/pressure |
| F14.9 | Water treatment plants | small/medium/large water treatment plants clos... |
| F14.10 | Water treatment plants | small/medium/large water treatment plants clos... |
| F15.1 | Water well | Wells |
| F16.1 | Transmission and distribution pipelines | Exposed transmission pipeline crossing |
| F16.2 | Transmission and distribution pipelines | Buried transmission pipeline crossing |
| F16.3 | Transmission and distribution pipelines | Pipelines (non-crossing) |
| F17.1 | Others | Pumping plants (small) below grade |
| F17.2 | Others | Pumping plants (medium/large) below grade |
| F17.3 | Others | Pumping plants (small) above grade |
| F17.4 | Others | Pumping plants (medium/large) above grade |
| F17.5 | Others | Pumping plants |
| F17.6 | Others | Control vaults and stations |
| F18.1 | Wastewater treatment plant | Small wastewater treatment plants |
| F18.2 | Wastewater treatment plant | Medium wastewater treatment plants |
| F18.3 | Wastewater treatment plant | Large wastewater treatment plants |
| F18.4 | Wastewater treatment plant | Small, medium, large Wastewater treatment plan... |
| F18.5 | Wastewater treatment plant | Small, medium, large Wastewater treatment plan... |
| F18.6 | Wastewater treatment plant | NaN |
| F19.1 | Transmission and distribution pipelines | Sewers & Interceptors: Exposed collector river... |
| F19.2 | Transmission and distribution pipelines | Sewers & Interceptors: Buried collector river ... |
| F19.3 | Transmission and distribution pipelines | Sewers & Interceptors: Pipes (non-crossings) |
| F20.1 | Others | Control vaults and control stations |
| F20.2 | Others | Lift stations: lift station (small), wet well/... |
| F20.3 | Others | Lift stations: lift station (medium large), we... |
| F20.4 | Others | Lift stations: lift station (small), submersible |
| F20.5 | Others | Lift stations: lift station (medium large), su... |
| F21.1 | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.2 | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.3 | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.3a (lower boundary) | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.3b (upper boundary) | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.4 | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.4a (lower boundary) | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.4b (upper boundary) | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.5 | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.5a (lower boundary) | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.5b (upper boundary) | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.6 | Education & Health | Generalized curve for commercial buildings, wh... |
| F21.7 | School buildings | Generalized curve for companies (incl. governm... |
| F21.8 | Health & education | Generalized curve for offices in diked areas |
| F21.9 | Education & Health | Generalized function for industry |
| F21.10 | School buildings | Generalized curve for buildings |
| F21.11 | School buildings | NaN |
| F21.12 | Health buildings | NaN |
| F21.13 | Education buildings | NaN |
And select a curve to use for each different subtype we are analysing. We select a curve from the database provided by Nirandjan et al. (2024).
ci_system = DICT_CIS_VULNERABILITY_FLOOD['main_roads']
selected_curves = []
for subtype in ci_system:
selected_curves.append(ci_system[subtype][0])
damage_curves = vul_df[['ID number']+selected_curves]
damage_curves = damage_curves.iloc[4:125,:]
damage_curves.set_index('ID number',inplace=True)
damage_curves.index = damage_curves.index.rename('Depth')
damage_curves = damage_curves.astype(np.float32)
damage_curves.columns = list(ci_system.keys())
damage_curves = damage_curves.fillna(method='ffill')
Make sure we set the index of the damage curves (the inundation depth) in the same metric as the hazard data (e.g. meters or centimeters). In our case, the hazard data is in centimeters, so we have to multiply our inundation values of the damage curves by 100.
damage_curves.index = damage_curves.index*100
damage_curves
| motorway | motorway_link | trunk | trunk_link | primary | primary_link | secondary | secondary_link | tertiary | tertiary_link | |
|---|---|---|---|---|---|---|---|---|---|---|
| Depth | ||||||||||
| 0 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 5.0 | 0.00200 | 0.00200 | 0.00200 | 0.00200 | 0.00200 | 0.00200 | 0.012000 | 0.012000 | 0.012000 | 0.012000 |
| 10.0 | 0.00400 | 0.00400 | 0.00400 | 0.00400 | 0.00400 | 0.00400 | 0.024000 | 0.024000 | 0.024000 | 0.024000 |
| 15.0 | 0.00600 | 0.00600 | 0.00600 | 0.00600 | 0.00600 | 0.00600 | 0.036000 | 0.036000 | 0.036000 | 0.036000 |
| 20.0 | 0.00800 | 0.00800 | 0.00800 | 0.00800 | 0.00800 | 0.00800 | 0.048000 | 0.048000 | 0.048000 | 0.048000 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 580.0 | 0.21500 | 0.21500 | 0.21500 | 0.21500 | 0.21500 | 0.21500 | 0.346500 | 0.346500 | 0.346500 | 0.346500 |
| 585.0 | 0.21625 | 0.21625 | 0.21625 | 0.21625 | 0.21625 | 0.21625 | 0.347375 | 0.347375 | 0.347375 | 0.347375 |
| 590.0 | 0.21750 | 0.21750 | 0.21750 | 0.21750 | 0.21750 | 0.21750 | 0.348250 | 0.348250 | 0.348250 | 0.348250 |
| 595.0 | 0.21875 | 0.21875 | 0.21875 | 0.21875 | 0.21875 | 0.21875 | 0.349125 | 0.349125 | 0.349125 | 0.349125 |
| 600.0 | 0.22000 | 0.22000 | 0.22000 | 0.22000 | 0.22000 | 0.22000 | 0.350000 | 0.350000 | 0.350000 | 0.350000 |
121 rows × 10 columns
Ancilliary data for processing#
world = gpd.read_file("https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp")
world_plot = world.to_crs(3857)
admin1 = gpd.read_file("https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_1_states_provinces.shp")
3. Performing the Damage Assessment#
We will use the DamageScanner approach. This is a fully optimised exposure, vulnerability damage and risk calculation method, that can capture a wide range of inputs to perform such an assessment.
total_damage_results = {}
for country_iso3 in countries_iso3:
country_bounds = world.loc[world.ADM0_ISO == country_iso3].bounds
country_geom = world.loc[world.ADM0_ISO == country_iso3].geometry
hazard_country = hazard_map.rio.clip_box(minx=country_bounds.minx.values[0],
miny=country_bounds.miny.values[0],
maxx=country_bounds.maxx.values[0],
maxy=country_bounds.maxy.values[0]
)
features = features_per_country[country_iso3]
total_damage_results[country_iso3] = DamageScanner(hazard_country, features, damage_curves, maxdam).calculate()
Overlay raster with vector: 100%|██████████████████████████████████████████████████████| 44/44 [00:18<00:00, 2.39it/s]
convert coverage to meters: 100%|██████████████████████████████████████████████████| 325/325 [00:00<00:00, 1109.28it/s]
Calculating damage: 100%|█████████████████████████████████████████████████████████| 325/325 [00:00<00:00, 14310.97it/s]
Overlay raster with vector: 100%|██████████████████████████████████████████████████████| 91/91 [00:24<00:00, 3.72it/s]
convert coverage to meters: 100%|███████████████████████████████████████████████████| 806/806 [00:01<00:00, 635.80it/s]
Calculating damage: 100%|██████████████████████████████████████████████████████████| 806/806 [00:00<00:00, 8582.73it/s]
Overlay raster with vector: 100%|████████████████████████████████████████████████████| 117/117 [00:56<00:00, 2.06it/s]
convert coverage to meters: 100%|█████████████████████████████████████████████████| 1032/1032 [00:01<00:00, 880.22it/s]
Calculating damage: 100%|████████████████████████████████████████████████████████| 1032/1032 [00:00<00:00, 9463.15it/s]
damage_results_countries = pd.concat(total_damage_results)
4. Save the Results#
For further analysis, we can save the results in their full detail, or save summary estimates per subnational unit. We show how to do this using the damage results.
hazard = 'river_flood'
return_period = '1_100'
for country_iso3,country_full_name in zip(countries_iso3,countries_full_names):
total_damage_results[country_iso3].to_file(f'Road_Damage_{country_full_name}_{hazard}_{return_period}.gpkg')
admin1_country = admin1.loc[admin1.sov_a3.isin(countries_iso3)]
damage_results_countries = damage_results_countries.sjoin(admin1_country[['adm1_code','name','geometry']])
admin1_damage = admin1_country.merge(damage_results_countries[['name','damage']].groupby('name').sum(),
left_on='name',
right_on='name',
how='outer')[['name','adm1_code','geometry','damage']]
admin1_damage.to_file(f'Admin1_Road_Damage_{hazard}_{return_period}.gpkg')
5. Visualizing the Results#
The results of the damage assessment can be visualized using charts and maps. This will provide a clear representation of which infrastructure is most affected by the hazard and the expected damage levels.
Plot all damages
fig, ax = plt.subplots(1,1,figsize=(10, 5))
damage_results_countries.to_crs(3857).plot(ax=ax,column='damage',cmap='Reds',legend=True,
legend_kwds={'shrink': 0.5,'label':'Absolute Damage'},zorder=5)
world_plot.loc[world.SOV_A3.isin(countries_iso3)].plot(ax=ax,facecolor="none",edgecolor='black')
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron,alpha=0.5)
ax.set_axis_off()
viz(damage_results_countries)
9. Performing the Risk Assessment#
To do so, we need to select the return periods we want to include, and create a dictionary as input. We will create this below.
return_periods = [2,5,10,50,100,200,500,1000]
total_risk_results = {}
for country_iso3 in countries_iso3:
country_bounds = world.loc[world.ADM0_ISO == country_iso3].bounds
country_geom = world.loc[world.ADM0_ISO == country_iso3].geometry
hazard_country = hazard_map.rio.clip_box(minx=country_bounds.minx.values[0],
miny=country_bounds.miny.values[0],
maxx=country_bounds.maxx.values[0],
maxy=country_bounds.maxy.values[0]
)
features = features_per_country[country_iso3]
hazard_dict = {}
for return_period in return_periods:
hazard_map = xr.open_dataset(f"https://hazards-data.unepgrid.ch/global_pc_h{return_period}glob.tif", engine="rasterio")
hazard_dict[return_period] = hazard_map.rio.clip_box(minx=country_bounds.minx.values[0],
miny=country_bounds.miny.values[0],
maxx=country_bounds.maxx.values[0],
maxy=country_bounds.maxy.values[0]
)
total_risk_results[country_iso3] = DamageScanner(hazard_country, features, damage_curves, maxdam).risk(hazard_dict)
10. Visualize the Risk Assessment#
To understand the relative shares of each component of the system, we create a pie chart with the percentages of risk to each component.
risk_results = pd.concat(total_risk_results)
risk_per_object = risk_results[['object_type','risk']].groupby('object_type').sum()
risk_per_object
| risk | |
|---|---|
| object_type | |
| primary | 4558.844763 |
| primary_link | 2.141788 |
| secondary | 21127.809544 |
| secondary_link | 0.989814 |
| tertiary | 52417.900611 |
| tertiary_link | 466.775363 |
| trunk | 2989.609524 |
| trunk_link | 4.575081 |
And create some figures. We start with a figure of all countries:
colors = ['#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00']
labels = [x.capitalize() for x in list(risk_per_object.index)]
sizes = risk_per_object.risk.values
pie = plt.pie(sizes,autopct='%1.1f%%', labeldistance=1.15, wedgeprops = { 'linewidth' : 1, 'edgecolor' : 'white' }, colors=colors);
plt.axis('equal')
plt.legend(loc = 'right', labels=labels,bbox_to_anchor=(1.25, 0.5),)
plt.title(f'Road Transport River Flood Risk for {', '.join(countries_full_names)}',fontweight='bold')
Text(0.5, 1.0, 'Road Transport River Flood Risk for Togo, Benin, Ghana')
And here an overview per country, in a single plot
country_risk_object = risk_results.reset_index()[['level_0','object_type','risk']].groupby(['level_0','object_type']).sum()
fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharey=True)
fig.suptitle('Road Transport River Flood Risk for Selected Countries', fontweight='bold')
all_labels = list(country_risk_object.index.get_level_values(1).unique())
inputs_for_legend = []
for ax, (country, risks) in zip(axes, country_risk_object.groupby('level_0')):
labels = [x.capitalize() for x in list(risks.index.get_level_values(1))]
if len(labels) == len(all_labels):
inputs_for_legend.append(labels)
sizes = risks.risk.values
pie = ax.pie(sizes, autopct='%1.1f%%', labeldistance=1.15,
wedgeprops={'linewidth': 1, 'edgecolor': 'white'}, colors=colors)
if len(pie[0]) == len(all_labels):
inputs_for_legend.append(pie[0])
ax.set_title(f'Risk for {country}', fontweight='bold')
ax.axis('equal')
# Add a legend to the right of the figure
fig.legend(inputs_for_legend[1], inputs_for_legend[0], loc='center right', bbox_to_anchor=(1.05, 0.5))
plt.tight_layout(rect=[0, 0, 0.9, 1]) # Adjust layout to fit legend
