RA2CE feature: Damages analysis (Huizinga)

Three types of damage curves can be introduced to this analysis: - Huizinga - OSdaMage - User defined

The Huizinga damage functions

This type of damage functions are initially developed for grid-based models, and generates damage per m2. For use in raster-based models, Huizinga (2007) developed a set of damage functions for diverse land use classes including transport infrastructure, initially for the European Union and later generalized worldwide (Huizinga et al., 2017). ‘Infrastructure’ is defined as physical damage to “roads and railroads as a result of contact with (fast flowing) water” (Huizinga et al., 2017, p.96).

This type is adapted by Van Ginkel et al. (2021) for use object-based models such as RA2CE, which can be implemented easily across the globe. The grid-based Huizinga (2007, 2017) infrastructure damage function is expressed in euros damage per inundated area (€/m2). The maximum damage for road infrastructure is 25 €/m2. To apply this function in the object-based model, they are multiplied by typical road widths (m) to obtain damage functions per unit road length

It is good to mention that this method is known to (strongly) underestimate damage to high-end road infrastructure such as motorways with street lighting (Van Ginkel et al., 2021); (strongly) overestimate damage to cheap roads in the underlying road networks, such small residential roads, unpaved roads and tracks; underestimate damage to bridges (Jongman et al., 2012).

An example for the Huizinga damage function

In this notebook, we provide an example for an event-based object-oriented Huizinga analysis.

import geopandas as gpd
import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np
import rasterio

root_dir = Path("data", "damages_analysis_huizinga")
assert root_dir.exists(), "root_dir not found."

Introducing a hazardous event

To use the flood map with RA2CE, we need to fill in the [hazard] section in the network.ini.

Specify the flood map name in the hazard_map parameter in network.ini. RA2CE expects the flood map to be located in the hazard folder. The aggregate_wl parameter in analysis.ini can be set to either ‘max’, ‘min’ or ‘mean’ to take the maximum, minimum or mean water depth per road segment when the exposure of the roads to a certain hazard (map) is determined.

Set the right CRS for the flood map in the hazard_crs parameter. This CRS can be different from the origins, destinations and road network. RA2CE will reproject the network to the CRS of the flood map and will reproject the road back to the original CRS when the CRS differs.

hazard_folder = root_dir / "static" / "hazard" # find the hazard folder where you locate your floo dmap
hazard_map = hazard_folder / "max_flood_depth.tif" # set the location of the hazard map

# Open the TIF file using rasterio
with rasterio.open(hazard_map) as src:
    # Read the TIF file as a numpy array
    tif_array = src.read(1)  # Change the band index (1) if necessary

plt.figure(figsize=(10, 10))
plt.imshow(tif_array, cmap='Blues')  # Change the colormap if desired
plt.colorbar(label='Pixel Values')
plt.title('Flood map')
plt.show()

Specifying the .ini files

Network.ini content > [project] name = beira [network] directed = False source = OSM download primary_file = None diversion_file = None file_id = rfid_c polygon = region_polygon.geojson network_type = drive road_types = motorway,motorway_link,primary,primary_link,secondary,secondary_link,tertiary,tertiary_link,residential save_gpkg = True [origins_destinations] origins = None destinations = None origins_names = None destinations_names = None id_name_origin_destination = None origin_count = None origin_out_fraction = None category = categorybr> [hazard] hazard_map = max_flood_depth.tif hazard_id = None hazard_field_name = waterdepth aggregate_wl = max hazard_crs = EPSG:32736br> [cleanup] snapping_threshold = None segmentation_length = 100 merge_lines = True merge_on_id = False cut_at_intersections = False

We now need to update our analysis initialisation files using the preferred OD-analysis (there are multiple). We will consider the damages analysis. With the aggregate_wl parameter, the user can choose which type of aggregation of the water level on the road segment (max, mean, min) the analysis should consider. For the damages analysis, the aggregate_wl=mean makes sense.

The damage_curve defines the damage curve type. ‘HZ’ to use the Huizinga damage function, ‘OSD’ to use the OSdaMage functions, and ‘MAN’ to use damage functions from manually inserted files. The event_type defines the type of the hazardous event, which is either ‘event’, or ‘return_period’. The former is an one-time event, while the later will be applied for the events witha probability of occurance.

[analysis1] name = Huizinga_damages_event analysis = damages event_type = event damage_curve = HZ save_csv = True save_gpkg =True

Set the paths to the initialization files and check if the files exist.

from ra2ce.ra2ce_handler import Ra2ceHandler
from ra2ce.network.network_config_data.network_config_data_reader import NetworkConfigDataReader
from ra2ce.analysis.analysis_config_data.analysis_config_data_reader import AnalysisConfigDataReader

# Load network data.
_network_config_data = NetworkConfigDataReader().read(root_dir.joinpath("network.ini"))

# Load analysis data.
_analysis_config_data = AnalysisConfigDataReader().read(root_dir.joinpath("analysis.ini"))

_analysis_config_data.input_path = root_dir.joinpath("input_analysis_data")

Run RA2CE.

handler = Ra2ceHandler.from_config(_network_config_data, _analysis_config_data)

handler.configure()
handler.run_analysis()

Visualising the results

analysis_output_path = root_dir / "output" / "damages"
damage_gdf = gpd.read_file(analysis_output_path / 'Huizinga_damages_event_segmented.gpkg')
damage_gdf.head() #show the origins

null_mask = damage_gdf['dam_EV1_HZ'].isnull()
damage_gdf.loc[null_mask, 'dam_EV1_HZ'] = damage_gdf.loc[null_mask, 'dam_EV1_HZ'].replace({np.nan: 0})


import folium
from branca.colormap import LinearColormap

# Create a colormap
colormap = LinearColormap(['yellow', 'orange', 'red'], vmin=0, vmax=damage_gdf['dam_EV1_HZ'].max(), caption='dam_EV1_HZ')

# Create a Folium map centered around the mean of the geometry
center = [damage_gdf['geometry'].centroid.y.mean(), damage_gdf['geometry'].centroid.x.mean()]
m = folium.Map(location=center, zoom_start=12, control_scale=True, tiles="CartoDB dark_matter")

# Create the color strip legend with white font color
colormap.caption = 'dam_EV1_HZ'
colormap.add_to(m)
colormap._repr_html_ = colormap._repr_html_().replace(
    'background-color: #000', 'background-color: #000; color: #000'
)

# Add edges to the map
for idx, row in damage_gdf.iterrows():
    color = 'darkgray' if row['dam_EV1_HZ'] == 0 else colormap(row['dam_EV1_HZ'])

    # Extracting coordinates from LineString
    coordinates = [(coord[1], coord[0]) for coord in row['geometry'].coords]

    # Create a popup with data
    popup_content = f"<strong>dam_EV1_HZ:</strong> {row['dam_EV1_HZ']}"
    popup = folium.Popup(popup_content, max_width=300)

    folium.PolyLine(
        locations=coordinates,
        color=color,
        weight=2,
        opacity=1,
        popup=popup
    ).add_to(m)

# Add colormap to the map
colormap.add_to(m)

# Display the map
m