Fetch rehabilitation costs during damage calculation

src/open_gira/direct_damages.py

@@ -6,6 +6,10 @@
 import re
 from typing import Union
 
+import pandas as pd
+import snkit
+
+from open_gira import fields
 from open_gira.utils import natural_sort
 
 
@@ -214,3 +218,79 @@ def get_rp_map(name: str) -> ReturnPeriodMap:
 
     # return a concrete subclass of ReturnPeriodMap
     return prefix_class_map[prefix](name)
+
+
+def rail_rehab_cost(row: pd.Series, rehab_costs: pd.DataFrame) -> float:
+    """
+    Determine the cost of rehabilitation for a given rail segment (row).
+
+    Args:
+        row: Railway segment
+        rehab_costs: Table of rehabilitation costs for various rail types
+
+    Returns:
+        Minimum cost, maximum cost, units of cost
+    """
+
+    # bridge should be a boolean type after data cleaning step
+    if row.bridge:
+        asset_type = "bridge"
+    else:
+        asset_type = "rail"
+
+    return float(rehab_costs.loc[rehab_costs["asset_type"] == asset_type, "rehab_cost_USD_per_km"])
+
+
+def road_rehab_cost(row: pd.Series, rehab_costs: pd.DataFrame) -> float:
+    """
+    Fetch the cost of rehabilitation for a given road segment (row).
+
+    Args:
+        row: Road segment
+        rehab_costs: Table of rehabilitation costs for various road types
+
+    Returns:
+        Cost in units of USD per km per lane.
+    """
+
+    # bridge should be boolean after data cleaning step
+    if row.bridge:
+        highway_type = "bridge"
+    else:
+        highway_type = row.tag_highway
+
+    if row.paved:
+        condition = "paved"
+    else:
+        condition = "unpaved"
+
+    return float(
+        rehab_costs["rehab_cost_USD_per_km_per_lane"].loc[
+            (rehab_costs.highway == highway_type) & (rehab_costs.road_cond == condition)
+        ].squeeze()
+    )
+
+
+def annotate_rehab_cost(network: snkit.network.Network, network_type: str, rehab_cost: pd.DataFrame) -> snkit.network.Network:
+    """
+    Label a network with rehabilitation costs to be used in subsequent direct
+    damage cost estimate.
+
+    Args:
+        network: Network with edges to label.
+        network_type: Category string, currently either road or rail.
+        rehab_cost: Table containing rehabilitation cost data per km. See lookup
+            functions for more requirements.
+
+    Returns:
+        Network labelled with rehabilitation cost data.
+    """
+
+    if network_type == "road":
+        network.edges[fields.REHAB_COST] = network.edges.apply(road_rehab_cost, axis=1, args=(rehab_cost,)) * network.edges.lanes
+    elif network_type == "rail":
+        network.edges[fields.REHAB_COST] = network.edges.apply(rail_rehab_cost, axis=1, args=(rehab_cost,))
+    else:
+        raise ValueError(f"No lookup function available for {network_type=}")
+
+    return network

workflow/transport-flood/event_set_direct_damages.py

@@ -3,17 +3,16 @@
 fraction for exposed assets.
 """
 
-from collections import defaultdict
 import logging
 import sys
 import warnings
 
 import geopandas as gpd
 import pandas as pd
 from scipy.interpolate import interp1d
-from scipy.integrate import simpson
 
 from open_gira import fields
+from open_gira.direct_damages import annotate_rehab_cost
 from open_gira.io import write_empty_frames, read_damage_curves
 from open_gira.utils import natural_sort
 
@@ -26,6 +25,7 @@
         damage_fraction_path: str = snakemake.output["damage_fraction"]
         damage_cost_path: str = snakemake.output["damage_cost"]
         damage_curves_dir: str = snakemake.config["direct_damages"]["curves_dir"]
+        rehabilitation_costs_path = snakemake.config["transport"]["rehabilitation_costs_path"]
         network_type: str = snakemake.params["network_type"]
         hazard_type: str = snakemake.params["hazard_type"]
         asset_types: set[str] = set(snakemake.config["direct_damages"]["asset_types"])
@@ -60,6 +60,10 @@
             write_empty_frames(path)
         sys.exit(0)  # exit gracefully so snakemake will continue
 
+    logging.info("Annotate network with rehabilitation costs")
+    rehab_cost = pd.read_excel(rehabilitation_costs_path, sheet_name=network_type)
+    exposure = annotate_rehab_cost(exposure, network_type, rehab_cost)
+
     # column groupings for data selection
     hazard_columns = [col for col in exposure.columns if col.startswith(fields.HAZARD_PREFIX)]
     non_hazard_columns = list(set(exposure.columns) - set(hazard_columns))

workflow/transport-flood/return_period_direct_damages.py

@@ -29,6 +29,7 @@
         expected_annual_damages_path: str = snakemake.output["expected_annual_damages"]
         return_period_and_ead_path: str = snakemake.output["return_period_and_ead"]
         damage_curves_dir: str = snakemake.config["direct_damages"]["curves_dir"]
+        rehabilitation_costs_path = snakemake.config["transport"]["rehabilitation_costs_path"]
         network_type: str = snakemake.params["network_type"]
         hazard_type: str = snakemake.params["hazard_type"]
         asset_types: set[str] = set(snakemake.config["direct_damages"]["asset_types"])
@@ -65,6 +66,10 @@
             write_empty_frames(path)
         sys.exit(0)  # exit gracefully so snakemake will continue
 
+    logging.info("Annotate network with rehabilitation costs")
+    rehab_cost = pd.read_excel(rehabilitation_costs_path, sheet_name=network_type)
+    exposure = annotate_rehab_cost(exposure, network_type, rehab_cost)
+
     # column groupings for data selection
     hazard_columns = [col for col in exposure.columns if col.startswith(fields.HAZARD_PREFIX)]
     non_hazard_columns = list(set(exposure.columns) - set(hazard_columns))

workflow/transport/create_rail_network.py

@@ -9,8 +9,6 @@
 import warnings
 
 import geopandas as gpd
-import pandas as pd
-import snkit
 
 from utils import annotate_country, get_administrative_data
 from open_gira.assets import RailAssets
@@ -19,71 +17,14 @@
 from open_gira.utils import str_to_bool
 
 
-def annotate_rehabilitation_costs(
-    network: snkit.network.Network, rehab_costs: pd.DataFrame
-) -> snkit.network.Network:
-
-    def get_rehab_costs(row: pd.Series, rehab_costs: pd.DataFrame) -> float:
-        """
-        Determine the cost of rehabilitation for a given rail segment (row).
-
-        Args:
-            row (pd.Series): Road segment
-            rehab_costs: (pd.DataFrame): Table of rehabilitation costs for various rail types
-
-        Returns:
-            Tuple[float, float, str]: Minimum cost, maximum cost, units of cost
-        """
-
-        # bridge should be a boolean type after data cleaning step
-        if row.bridge:
-            asset_type = "bridge"
-        else:
-            asset_type = "rail"
-
-        return float(rehab_costs.loc[rehab_costs["asset_type"] == asset_type, "rehab_cost_USD_per_km"])
-
-    # lookup costs
-    network.edges["rehab_cost_USD_per_km"] = network.edges.apply(
-        get_rehab_costs, axis=1, args=(rehab_costs,)
-    )
-
-    return network
-
-
 if __name__ == "__main__":
 
-    try:
-        osm_edges_path = snakemake.input["edges"]  # type: ignore
-        osm_nodes_path = snakemake.input["nodes"]  # type: ignore
-        administrative_data_path = snakemake.input["admin"]  # type: ignore
-        nodes_output_path = snakemake.output["nodes"]  # type: ignore
-        edges_output_path = snakemake.output["edges"]  # type: ignore
-        slice_number = snakemake.params["slice_number"]  # type: ignore
-        rehabilitation_costs_path = snakemake.config["transport"]["rehabilitation_costs_path"]  # type: ignore
-
-    except NameError:
-        # If "snakemake" doesn't exist then must be running from the
-        # command line.
-        (
-            osm_edges_path,
-            osm_nodes_path,
-            administrative_data_path,
-            nodes_output_path,
-            edges_output_path,
-            slice_number,
-            rehabilitation_costs_path,
-            transport_costs_path,
-            flow_cost_time_factor,
-        ) = sys.argv[1:]
-
-        # osm_edges_path = ../../results/geoparquet/tanzania-latest_filter-rail/slice-0_edges.geoparquet
-        # osm_nodes_path = ../../results/geoparquet/tanzania-latest_filter-rail/slice-0_nodes.geoparquet
-        # administrative_data_path = ../../results/input/admin-boundaries/gadm36_levels.gpkg
-        # nodes_output_path = ../../results/geoparquet/tanzania-latest_filter-rail/slice-0_nodes.geoparquet
-        # edges_output_path = ../../results/geoparquet/tanzania-latest_filter-rail/slice-0_edges.geoparquet
-        # slice_number = 0
-        # rehabilitation_costs_path = ../../bundled_data/rehabilitation.xlsx
+    osm_edges_path = snakemake.input["edges"]  # type: ignore
+    osm_nodes_path = snakemake.input["nodes"]  # type: ignore
+    administrative_data_path = snakemake.input["admin"]  # type: ignore
+    nodes_output_path = snakemake.output["nodes"]  # type: ignore
+    edges_output_path = snakemake.output["edges"]  # type: ignore
+    slice_number = int(snakemake.params["slice_number"])  # type: ignore
 
     slice_number = int(slice_number)
     osm_epsg = 4326
@@ -144,12 +85,6 @@ def get_rehab_costs(row: pd.Series, rehab_costs: pd.DataFrame) -> float:
         get_administrative_data(administrative_data_path, to_epsg=osm_epsg),
     )
 
-    logging.info("Annotate network with rehabilitation costs")
-    network = annotate_rehabilitation_costs(
-        network,
-        pd.read_excel(rehabilitation_costs_path, sheet_name="rail")
-    )
-
     logging.info("Writing network to disk")
     network.edges.to_parquet(edges_output_path)
     network.nodes.to_parquet(nodes_output_path)