diff --git a/config/config.yaml b/config/config.yaml index 0361825f..065f954a 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -1,6 +1,6 @@ -################################## -### TRANSPORT DAMAGES WORKFLOW ### -################################## +######################### +### TRANSPORT DAMAGES ### +######################### # OSM datasets in PBF format, principally from: https://download.geofabrik.de/ # infrastructure_datasets: @@ -66,9 +66,9 @@ keep_tags: # Number of slices to cut dataset into -- must be a square number slice_count: 64 -##################################### -### TRANSPORT DISRUPTION WORKFLOW ### -##################################### +#################################### +### MULTI-MODAL NETWORK CREATION ### +#################################### # country for which trade OD has been prepared, and we are to route land trade flows study_country_iso_a3: "THA" @@ -90,6 +90,10 @@ intermodal_cost_USD_t: maritime_road: 4 maritime_rail: 5 +############################ +### TRANSPORT DISRUPTION ### +############################ + # drop trade flows with less volume than this (accelerate flow allocation) # N.B. 50t threshold preserves 91% of total volume and 88% of total value # for combined agriculture & manufacturing flows between Thailand road nodes -> partner countries @@ -98,9 +102,9 @@ minimum_flow_volume_t: 50 # if disrupting a network, remove edges experiencing hazard values in excess of this edge_failure_threshold: 0.5 -######################### -### FLOODING WORKFLOW ### -######################### +########################## +### TRANSPORT FLOODING ### +########################## hazard_datasets: # hazard rasters to retrieve @@ -131,9 +135,9 @@ direct_damages: ] -############################### -### STORM / ENERGY WORKFLOW ### -############################### +################################## +### CYCLONE / ELECTRICITY GRID ### +################################## # sets of storm ids to process for potentially many country networks storm_sets: diff --git a/docs/src/user-guide/usage/network-creation/composite-networks.md b/docs/src/user-guide/usage/network-creation/composite-networks.md new file mode 100644 index 00000000..11fc4827 --- /dev/null +++ b/docs/src/user-guide/usage/network-creation/composite-networks.md @@ -0,0 +1,63 @@ +# Composite network generation + +Create road or rail networks with varying asset filters (to primary routes in +one region, tertiary in others, etc.). + +## Description + +For each `infrastructure_dataset` and `network_filter`: +1. Download or copy `.osm.pbf` file to input data directory. +1. Filter OSM file with `osmium tags-filter` for elements matching the filter file (see configuration, below). +1. Define a bounding box for the OSM file, write to disk as JSON. +1. Define a grid partitioning the bounding box into a square number of 'slices' as a series of JSON files, one for each slice. +1. Cut the OSM file into these slices. +1. Convert the sliced OSM files into geoparquet, retaining the `keep_tags` as configured. +1. Clean and annotate features in the geoparquet files (joining additional data such as country, rehabiliation costs, etc.). +1. Join sliced network components together. +1. Join all datasets together. + +## Configuration + +- Review and amend `config/composite_networks/*.csv` + - Filter files should include two columns, infrastructure_dataset and + network_filter. + For example: + ```bash + infrastructure_dataset,network_filter + thailand-latest,road-secondary + laos-latest,road-primary + cambodia-latest,road-primary + myanmar-latest,road-primary + ``` + These then map to the infrastructure datasets and network filters + specified in the `config/config.yaml` file. +- Review and amend `config/config.yaml`: + - The `composite_network` mapping is from an identifier key to a filter file + path. + - The `infrastructure_datasets` map should contain a key pointing to an `.osm.pbf` + file URL for the desired areas. + There are currently entries for the planet, for (some definition of) + continents and several countries. We use the + [geofabrik](http://download.geofabrik.de/) service for continent and + country-level OSM extracts. + - Check the OSM filter file pointed to by `network_filters.road`. + This file specifies which [elements](https://wiki.openstreetmap.org/wiki/Elements) + (nodes, ways or relations) to keep (or reject) from the multitude of data + in an OSM file. See the filter expressions section + [here](https://docs.osmcode.org/osmium/latest/osmium-tags-filter.html) + for more information on the syntax of these files. + - Check and amend `keep_tags.road` and/or `keep_tags.rail`. This list of + strings specifies which `tags` (attributes) to retain on the filtered + elements we extract from the `.osm.pbf` file. + - Review `slice_count`. This controls the degree of parallelism possible. + With it set to 1, there is no spatial slicing (we create the network in + a single chunk). To speed network creation for large domains, it can be + set to a larger square number. The first square number greater than your + number of available CPUs is a good heuristic. + +## Creation + +Here's an example creation command: +```bash +snakemake --cores all results/composite_network/south-east-asia-road/edges.gpq +``` diff --git a/docs/src/user-guide/usage/network-creation/multi-modal.md b/docs/src/user-guide/usage/network-creation/multi-modal.md new file mode 100644 index 00000000..986d6609 --- /dev/null +++ b/docs/src/user-guide/usage/network-creation/multi-modal.md @@ -0,0 +1,86 @@ +# Multi-modal network integration + +Combine different transport network modes (road, rail, maritime) together into a +connected global network, centered on a single country. + +## Description + +1. We require input land networks to connect. These can be created with other workflows in open-gira or you may bring your own. + +Using open-gira you can create a road network for Thailand with the following: +```bash +snakemake --cores all results/thailand-latest_filter-road-primary/edges.gpq +``` +See more details in [road.md](./road.md) and [rail.md](./rail.md). + +Or make a composite network from +```bash +snakemake --cores all results/composite_network/south-east-asia-road/edges.gpq +``` +For more details, see [composite-networks.md](./composite-networks.md). + +However they're created, input land networks should be placed in the following locations: +```python +road_network_nodes = "{OUTPUT_DIR}/input/networks/road/{PROJECT_SLUG}/nodes.gpq", +road_network_edges = "{OUTPUT_DIR}/input/networks/road/{PROJECT_SLUG}/edges.gpq", +rail_network_nodes = "{OUTPUT_DIR}/input/networks/rail/{PROJECT_SLUG}/nodes.gpq", +rail_network_edges = "{OUTPUT_DIR}/input/networks/rail/{PROJECT_SLUG}/edges.gpq", +``` + +1. We also require input maritime network data: + +Network nodes, including ports for import, export and trans-shipment. +```python +nodes = "{OUTPUT_DIR}/input/networks/maritime/nodes.gpq", +``` +Sample node data: +``` + id infra name iso3 Continent_Code geometry + port3 port Aberdeen_United Kingdom GBR EU POINT (-2.07662 57.13993) + port84 port Avonmouth_United Kingdom GBR EU POINT (-2.70877 51.49812) +port121 port Barry_United Kingdom GBR EU POINT (-3.24712 51.40467) +``` + +Network edges with costs (used for routing over). The referenced port nodes are +appended with `in` or `out` as the network is bi-directional. +```python +edges_no_geom = "{OUTPUT_DIR}/input/networks/maritime/edges_by_cargo/maritime_base_network_general_cargo.pq", +``` +Sample edge data: +``` + from_id to_id from_infra to_infra mode from_iso3 to_iso3 distance_km cost_USD_t_km from_port to_port + port0_out port1099_in port port maritime AUS ZAF 14296.116346 0.002385 port0 port1099 +port1001_out port1099_in port port maritime BRA ZAF 9016.208416 0.003095 port1001 port1099 +port1005_out port1099_in port port maritime ITA ZAF 10973.998701 0.004198 port1005 port1099 +``` + +Network edges without costs, but with a geometry that can be used for visualising flows over. +```python +edges_visualisation = "{OUTPUT_DIR}/input/networks/maritime/edges.gpq", +``` +Sample visualisation edge data: +``` + from_id to_id from_infra to_infra id geometry + maritime0 maritime1265 maritime maritime maritimeroute_0 LINESTRING (11.00000 -19.00002, 5.51647 -19.58... +maritime213 maritime1265 maritime maritime maritimeroute_1 LINESTRING (9.99999 -30.00001, 4.79550 -25.083... +maritime964 maritime1265 maritime maritime maritimeroute_2 LINESTRING (-10.00000 -10.00002, -7.58272 -12.... +``` + +## Configuration + +- Review and amend `config/config.yaml`: + - `study_country_iso_a3`: the ISO 3 letter code of the country in question + - `road_cost_USD_t_km`: cost in USD per tonne per km of road freight + - `road_cost_USD_t_h`: cost in USD per tonne per hour of road freight + - `road_default_speed_limit_km_h`: speed limit in km per hour if no other is available + - `rail_cost_USD_t_km`: cost in USD per tonne per km of rail freight + - `rail_cost_USD_t_h`: cost in USD per tonne per hour of rail freight + - `rail_average_freight_speed_km_h`: speed assumption for rail freight in km per hour + - `intermodal_cost_USD_t`: cost in USD per tonne of changing from one network mode to another (e.g. road to rail) + +## Creation + +Here's an example command to create a multi-modal network, bringing together land networks under `project-thailand` and joining them to a maritime network: +```bash +snakemake --cores all results/multi-modal_network/project-thailand/edges.gpq +``` \ No newline at end of file diff --git a/docs/src/user-guide/usage/network-creation/rail.md b/docs/src/user-guide/usage/network-creation/rail.md index 9e91b72b..346d1fe1 100644 --- a/docs/src/user-guide/usage/network-creation/rail.md +++ b/docs/src/user-guide/usage/network-creation/rail.md @@ -43,8 +43,6 @@ To specify a desired network: a single chunk). To speed network creation for large domains, it can be set to a larger square number. The first square number greater than your number of available CPUs is a good heuristic. - - Check and amend the values of `transport.rail`, which provide some - defaults for OSM data gap-filling. ## Creation diff --git a/docs/src/user-guide/usage/network-creation/road.md b/docs/src/user-guide/usage/network-creation/road.md index 5a8cf995..3d62872f 100644 --- a/docs/src/user-guide/usage/network-creation/road.md +++ b/docs/src/user-guide/usage/network-creation/road.md @@ -42,8 +42,6 @@ To specify a desired network: a single chunk). To speed network creation for large domains, it can be set to a larger square number. The first square number greater than your number of available CPUs is a good heuristic. - - Check and amend the values of `transport.road`, which provide some - defaults for OSM data gap-filling. ## Creation diff --git a/docs/src/user-guide/usage/risk-analysis/flow-allocation.md b/docs/src/user-guide/usage/risk-analysis/flow-allocation.md new file mode 100644 index 00000000..5b4028f9 --- /dev/null +++ b/docs/src/user-guide/usage/risk-analysis/flow-allocation.md @@ -0,0 +1,53 @@ +# Flow allocation + +This workflow can allocate transport flows from origin-destination matrices over +networks. These networks may be intact or degraded. + +## Description + +1. We require an input trade origin-destination matrix (OD). + +This should have the following format: +``` + id partner_GID_0 value_kusd volume_tons +thailand-latest_14_10 ABW 4.536817e-07 5.172909e-07 +thailand-latest_14_10 AFG 3.524201e-06 8.233424e-07 +thailand-latest_14_10 AGO 4.582354e-04 1.108041e-03 +``` + +Where: +- `id` are node ids in the network to be routed over +- `partner_GID_0` is a partner country ISO 3-letter code +- `value_kusd` is the trade value in thousands of USD per year +- `volume_tons` is the trade mass (a.k.a. volume) in tonnes per year + +This file should be located here: +```python +od = "{OUTPUT_DIR}/input/trade_matrix/{PROJECT_SLUG}/trade_nodes_total.parquet", +``` + +1. This workflow can route over intact or disrupted networks. If disrupting a +network, the raster used to represent the hazard must be available at: +```python +raster = "{OUTPUT_DIR}/hazard/{HAZARD_SLUG}.tif", +``` + +## Configuration + +Review and amend `config/config.yaml`: +- `minimum_flow_volume_t`: discard flows in the OD with a `volume_tons` value less than this +- `edge_failure_threshold`: when failing a network subject to a hazard raster, +remove edges experiencing intensities greater than this + + +## Creation + +Here's an example creation command where `PROJECT_SLUG` is `project-thailand`: +```bash +snakemake --cores all results/flow_allocation/project-thailand/routes_with_costs.pq", +``` + +And for the disrupted case, where `HAZARD_SLUG` is `hazard-thai-floods-2011`: +```bash +snakemake --cores all results/flow_allocation/project-thailand/hazard-thai-floods-2011/routes_with_costs.pq", +``` diff --git a/src/open_gira/disruption.py b/src/open_gira/disruption.py new file mode 100644 index 00000000..c9a65bd0 --- /dev/null +++ b/src/open_gira/disruption.py @@ -0,0 +1,67 @@ +import geopandas as gpd +import pandas as pd +import rasterio + +import snail.intersection + + +def filter_edges_by_raster( + edges: gpd.GeoDataFrame, + raster_path: str, + failure_threshold: float, + band: int = 1 +) -> gpd.GeoDataFrame: + """ + Remove edges from a network that are exposed to gridded values in excess of + a given threshold. + + Args: + edges: Network edges to consider. Must contain geometry column containing linestrings. + raster_path: Path to raster file on disk, openable by rasterio. + failure_threshold: Edges experiencing a raster value in excess of this will be + removed from the network. + band: Band of raster to read data from. + + Returns: + Network without edges experiencing raster values in excess of threshold. + """ + + # we will create a new edge_id column, fail if we would overwrite + assert "edge_id" not in edges.columns + + # split out edges without geometry as snail will not handle them gracefully + print("Parition edges by existence of geometry...") + no_geom_edges = edges.loc[edges.geometry.type.isin({None}), :].copy() + with_geom_edges = edges.loc[~edges.geometry.type.isin({None}), :].copy() + + # we need an id to select edges by + with_geom_edges["edge_id"] = range(len(with_geom_edges)) + + print("Read raster transform...") + grid = snail.intersection.GridDefinition.from_raster(raster_path) + + print("Prepare linestrings...") + edges = snail.intersection.prepare_linestrings(with_geom_edges) + + print("Split linestrings...") + splits = snail.intersection.split_linestrings(edges, grid) + + print("Lookup raster indices...") + splits_with_indices = snail.intersection.apply_indices(splits, grid) + + print("Read raster data into memory...") + with rasterio.open(raster_path) as dataset: + raster = dataset.read(band) + + print("Lookup raster value for splits...") + values = snail.intersection.get_raster_values_for_splits(splits_with_indices, raster) + splits_with_values = splits_with_indices.copy() + splits_with_values["raster_value"] = values + + print(f"Filter out edges with splits experiencing values in excess of {failure_threshold} threshold...") + failed_splits_mask = splits_with_values.raster_value > failure_threshold + failed_edge_ids = set(splits_with_values[failed_splits_mask].edge_id.unique()) + surviving_edges_with_geom = edges.loc[~edges.edge_id.isin(failed_edge_ids), :] + + print("Done filtering edges...") + return pd.concat([no_geom_edges, surviving_edges_with_geom]).sort_index().drop(columns=["edge_id"]) \ No newline at end of file diff --git a/src/open_gira/plot/utils.py b/src/open_gira/plot/utils.py index b73c1ecd..c8197d97 100644 --- a/src/open_gira/plot/utils.py +++ b/src/open_gira/plot/utils.py @@ -3,6 +3,60 @@ """ +import geopandas as gpd +import numpy as np +import pandas as pd +import shapely +from shapely.ops import split + + +def chop_at_antimeridian(gdf: gpd.GeoDataFrame, drop_null_geometry: bool = False) -> gpd.GeoDataFrame: + """ + Cut LineStrings either side of antimeridian, then drop the fragments that + intersect with antimeridian. + + Warning: Will create new rows (split geometries) with duplicate indices. + + Args: + gdf: Table with geometry to chop at antimeridian + drop_null_geometry: If true, drop any null geometry rows before plotting + + Returns: + Table, potentially with new rows. No rows in the table should have + geometries that cross the antimeridian. + """ + if drop_null_geometry: + gdf = gdf.loc[~gdf.geometry.isna(), :] + + assert set(gdf.geometry.type) == {'LineString'} + + def split_on_meridian(gdf: gpd.GeoDataFrame, meridian: shapely.geometry.LineString) -> gpd.GeoDataFrame: + return gdf.assign(geometry=gdf.apply(lambda row: split(row.geometry, meridian), axis=1)).explode(index_parts=False) + + xlim = 179.9 + ylim = 90 + + split_e = split_on_meridian(gdf, shapely.geometry.LineString([(xlim, ylim), (xlim, -ylim)])) + split_e_and_w = split_on_meridian(split_e, shapely.geometry.LineString([(-xlim, ylim), (-xlim, -ylim)])) + + def crosses_antimeridian(row: pd.Series) -> bool: + """ + Check if there are longitudes in a geometry that are near the antimeridian + (i.e. -180) and both sides of it. If so, return true. + """ + x, _ = row.geometry.coords.xy + longitudes_near_antimeridian = np.array(x)[np.argwhere(np.abs(np.abs(x) - 180) < xlim).ravel()] + if len(longitudes_near_antimeridian) == 0: + return False + hemispheres = np.unique(np.sign(longitudes_near_antimeridian)) + if (-1 in hemispheres) and (1 in hemispheres): + return True + else: + return False + + return split_e_and_w[~split_e_and_w.apply(crosses_antimeridian, axis=1)] + + def figure_size( min_x: float, min_y: float, diff --git a/src/open_gira/routing.py b/src/open_gira/routing.py new file mode 100644 index 00000000..740a54e9 --- /dev/null +++ b/src/open_gira/routing.py @@ -0,0 +1,222 @@ +""" +Allocate flows (value and volume) from an origin-destination (OD) file across a network of edges. +""" + +import multiprocessing +import os +import tempfile +import time + +import igraph as ig +import geopandas as gpd +import pandas as pd +from tqdm import tqdm + + +# dict containing: +# 'value_kusd' -> float +# 'volume_tons' -> float +# 'edge_indicies' -> list[indicies] +FlowResult = dict[str, float | list[int]] + +# dict with FlowResult values +# source node, destination country -> FlowResult dict +# e.g. +# ('thailand_4_123', 'GID_0_GBR') -> FlowResult dict +RouteResult = dict[tuple[str, str], FlowResult] + + +# We do not do routing within partner countries, instead we terminate at a port +# in the destination country. Destination country nodes are connected to ports by +# special edges. We give these edges a very high value, so that least cost route +# finding will take them as a last resort (no other lower cost route is +# available). + +# When calculating the total cost of a route, we want to be able to identify +# these imaginary links and remove them. +DESTINATION_LINK_COST_USD_T: float = 1E6 + + +def init_worker(graph_filepath: str, od_filepath: str) -> None: + """ + Create global variables referencing graph and OD to persist through worker lifetime. + + Args: + graph_filepath: Filepath of pickled igraph.Graph to route over. + od_filepath: Filepath to table of flows from origin node 'id' to + destination country 'partner_GID_0', should also contain 'value_kusd' + and 'volume_tons'. + """ + print(f"Process {os.getpid()} initialising...") + global graph + graph = ig.Graph.Read_Pickle(graph_filepath) + global od + od = pd.read_parquet(od_filepath) + return + + +def route_from_node(from_node: str) -> RouteResult: + """ + Route flows from single 'from_node' to destinations across graph. Record value and + volume flowing across each edge. + + Args: + from_node: Node ID of source node. + + Returns: + Mapping from (source node, destination country node) key, to value of + flow, volume of flow and list of edge ids of route. + """ + print(f"Process {os.getpid()} routing {from_node}...") + + from_node_od = od[od.id == from_node] + destination_nodes: list[str] = [f"GID_0_{iso_a3}" for iso_a3 in from_node_od.partner_GID_0.unique()] + + routes_edge_list = [] + try: + routes_edge_list: list[list[int]] = graph.get_shortest_paths( + f"road_{from_node}", + destination_nodes, + weights="cost_USD_t", + output="epath" + ) + except ValueError as error: + if "no such vertex" in str(error): + print(f"{error}... skipping destination") + pass + else: + raise error + + routes: RouteResult = {} + + if routes_edge_list: + assert len(routes_edge_list) == len(destination_nodes) + else: + return routes + + # lookup trade value and volume for each pairing of from_node and partner country + for i, destination_node in enumerate(destination_nodes): + # "GID_0_GBR" -> "GBR" + iso_a3 = destination_node.split("_")[-1] + route = from_node_od[ + (from_node_od.id == from_node) & (from_node_od.partner_GID_0 == iso_a3) + ] + value_kusd, = route.value_kusd + volume_tons, = route.volume_tons + + routes[(from_node, destination_node)] = { + "value_kusd": value_kusd, + "volume_tons": volume_tons, + "edge_indices": routes_edge_list[i] + } + + print(f"Process {os.getpid()} finished routing {from_node}...") + return routes + + +def route_from_all_nodes(od: pd.DataFrame, edges: gpd.GeoDataFrame, n_cpu: int) -> RouteResult: + """ + Route flows from origins to destinations across graph. + + Args: + od: Table of flows from origin node 'id' to destination country + 'partner_GID_0', should also contain 'value_kusd' and 'volume_tons'. + edges: Table of edges to construct graph from. First column should be + source node id and second destination node id. + n_cpu: Number of CPUs to use for routing. + + Returns: + Mapping from source node, to destination country node, to flow in value + and volume along this route and list of edge indices constituting + the route. + """ + + print("Creating graph...") + # cannot add vertices as edges reference port493_out, port281_in, etc. which are missing from nodes file + # use_vids=False as edges.from_id and edges_to_id are not integers + graph = ig.Graph.DataFrame(edges, directed=True, use_vids=False) + + temp_dir = tempfile.TemporaryDirectory() + + print("Writing graph to disk...") + graph_filepath = os.path.join(temp_dir.name, "graph.pickle") + graph.write_pickle(graph_filepath) + + print("Writing OD to disk...") + od_filepath = os.path.join(temp_dir.name, "od.pq") + od.to_parquet(od_filepath) + + print("Routing...") + start = time.time() + from_nodes = od.id.unique() + args = ((from_node,) for from_node in from_nodes) + # as each process is created, it will load the graph and od from disk in + # init_worker and then persist these in memory as globals between chunks + with multiprocessing.Pool( + processes=n_cpu, + initializer=init_worker, + initargs=(graph_filepath, od_filepath), + ) as pool: + routes: list[RouteResult] = pool.starmap(route_from_node, args) + + print("\n") + print(f"Routing completed in {time.time() - start:.2f}s") + + temp_dir.cleanup() + + # flatten our list of RouteResult dicts into one dict + return {k: v for item in routes for (k, v) in item.items()} + + +def lookup_route_costs( + routes_path: str, + edges_path: str, + destination_link_cost_USD_t: float = DESTINATION_LINK_COST_USD_T +) -> pd.DataFrame: + """ + For each route (source -> destination pair), lookup the edges + of the least cost route (the route taken) and sum those costs. + Store alongside value and volume of route. + + Args: + routes_path: Path to routes table, should have multi-index: (source node, + destination node) and include value_kusd, volume_tons and edge_indices + columns + edges_path: Path to edges table, should have cost_USD_t column which we + will positional index into with edge_indices from the routes table. + destination_link_cost_USD_t: Cost of traversing 'destination' links, to + partner entities. There should only be one of these links in any given + route. + + Returns: + Routes appended with their total cost in USD t-1 + """ + routes_with_edge_indices: pd.DataFrame = pd.read_parquet(routes_path) + edges: gpd.GeoDataFrame = gpd.read_parquet(edges_path) + cost_col_id = edges.columns.get_loc("cost_USD_t") + routes = [] + for index, route_data in tqdm(routes_with_edge_indices.iterrows(), total=len(routes_with_edge_indices)): + source_node, destination_node = index + cost_including_destination_link_USD_t = edges.iloc[route_data.edge_indices, cost_col_id].sum() + + cost_USD_t: float = cost_including_destination_link_USD_t % destination_link_cost_USD_t + + if int(cost_including_destination_link_USD_t // destination_link_cost_USD_t) != 1: + # must have exactly 1 destination link, otherwise not a valid route + continue + + if cost_USD_t != 0: + routes.append( + ( + source_node, + destination_node.split("_")[-1], + route_data.value_kusd, + route_data.volume_tons, + cost_USD_t + ) + ) + + return pd.DataFrame( + routes, + columns=["source_node", "destination_node", "value_kusd", "volume_tons", "cost_USD_t"] + ) diff --git a/tests/config/config.yaml b/tests/config/config.yaml index 1256fa9b..8813e3e3 100644 --- a/tests/config/config.yaml +++ b/tests/config/config.yaml @@ -55,6 +55,39 @@ direct_damages: ] +##################################### +### TRANSPORT DISRUPTION WORKFLOW ### +##################################### + +# country for which trade OD has been prepared, and we are to route land trade flows +study_country_iso_a3: "THA" + +# transport cost information +# road +road_cost_USD_t_km: 0.05 +road_cost_USD_t_h: 0.48 +road_default_speed_limit_km_h: 80 +# rail +rail_cost_USD_t_km: 0.05 +rail_cost_USD_t_h: 0.38 +rail_average_freight_speed_km_h: 40 + +# cost of changing transport mode in USD per tonne +# from mistral/ccg-critical-minerals/processed_data/transport_costs/intermodal.xlsx, 20240611 +intermodal_cost_USD_t: + road_rail: 5 + maritime_road: 4 + maritime_rail: 5 + +# drop trade flows with less volume than this (accelerate flow allocation) +# N.B. 50t threshold preserves 91% of total volume and 88% of total value +# for combined agriculture & manufacturing flows between Thailand road nodes -> partner countries +minimum_flow_volume_t: 50 + +# if disrupting a network, remove edges experiencing hazard values in excess of this +edge_failure_threshold: 0.5 + + ############################### ### STORM / ENERGY WORKFLOW ### ############################### diff --git a/workflow/Snakefile b/workflow/Snakefile index 616db8c0..0b31eb4d 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -134,6 +134,7 @@ include: "transport/create_overall_bbox.smk" include: "transport/join_data.smk" include: "transport/maritime/maritime.smk" include: "transport/multi-modal/multi_modal.smk" +include: "transport/flow_allocation/allocate.smk" include: "flood/aqueduct.smk" include: "flood/trim_hazard_data.smk" diff --git a/workflow/transport/flow_allocation/allocate.py b/workflow/transport/flow_allocation/allocate.py new file mode 100644 index 00000000..a5f3d113 --- /dev/null +++ b/workflow/transport/flow_allocation/allocate.py @@ -0,0 +1,48 @@ + +import geopandas as gpd +import pandas as pd +from tqdm import tqdm + +from open_gira.routing import route_from_all_nodes, RouteResult + + +if __name__ == "__main__": + + print("Reading network...") + # read in global multi-modal transport network + edges = gpd.read_parquet(snakemake.input.edges) + available_destinations = edges[edges["mode"] == "imaginary"].to_id.unique() + available_country_destinations = [d.split("_")[-1] for d in available_destinations if d.startswith("GID_")] + + print("Reading OD matrix...") + # read in trade OD matrix + od = pd.read_parquet(snakemake.input.od) + print(f"OD has {len(od):,d} flows") + + # 5t threshold drops THL road -> GID_0 OD from ~21M -> ~2M + minimum_flow_volume_tons = snakemake.config["minimum_flow_volume_t"] + od = od[od.volume_tons > minimum_flow_volume_tons] + print(f"After dropping flows with volume < {minimum_flow_volume_tons}t, OD has {len(od):,d} flows") + + # drop any flows we can't find a route to + od = od[od.partner_GID_0.isin(available_country_destinations)] + print(f"After dropping unrouteable destination countries, OD has {len(od):,d} flows") + + routes: RouteResult = route_from_all_nodes(od, edges, snakemake.threads) + + print("Writing routes to disk as parquet...") + pd.DataFrame(routes).T.to_parquet(snakemake.output.routes) + + print("Assigning route flows to edges...") + edges["value_kusd"] = 0 + edges["volume_tons"] = 0 + value_col_id = edges.columns.get_loc("value_kusd") + volume_col_id = edges.columns.get_loc("volume_tons") + for (from_node, destination_country), route_data in tqdm(routes.items()): + edges.iloc[route_data["edge_indices"], value_col_id] += route_data["value_kusd"] + edges.iloc[route_data["edge_indices"], volume_col_id] += route_data["volume_tons"] + + print("Writing edge flows to disk as geoparquet...") + edges.to_parquet(snakemake.output.edges_with_flows) + + print("Done") \ No newline at end of file diff --git a/workflow/transport/flow_allocation/allocate.smk b/workflow/transport/flow_allocation/allocate.smk new file mode 100644 index 00000000..065b3d0a --- /dev/null +++ b/workflow/transport/flow_allocation/allocate.smk @@ -0,0 +1,81 @@ +rule allocate_intact_network: + """ + Allocate a trade OD matrix across a multi-modal transport network + + Test with: + snakemake -c1 -- results/flow_allocation/project-thailand/edges.gpq + """ + input: + edges = "{OUTPUT_DIR}/multi-modal_network/{PROJECT_SLUG}/edges.gpq", + od = "{OUTPUT_DIR}/input/trade_matrix/{PROJECT_SLUG}/trade_nodes_total.parquet", + threads: workflow.cores + params: + # if this changes, we want to trigger a re-run + minimum_flow_volume_t = config["minimum_flow_volume_t"], + output: + routes = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/routes.pq", + edges_with_flows = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/edges.gpq", + script: + "./allocate.py" + + +rule allocate_degraded_network: + """ + Allocate a trade OD matrix across a multi-modal transport network which has + lost edges as a result of intersection with a hazard map. + + Test with: + snakemake -c1 -- results/flow_allocation/project-thailand/hazard-thai-floods-2011-JBA/edges.gpq + """ + input: + edges = "{OUTPUT_DIR}/multi-modal_network/{PROJECT_SLUG}/{HAZARD_SLUG}/edges.gpq", + od = "{OUTPUT_DIR}/input/trade_matrix/{PROJECT_SLUG}/trade_nodes_total.parquet", + threads: workflow.cores + params: + # if this changes, we want to trigger a re-run + minimum_flow_volume_t = config["minimum_flow_volume_t"], + output: + routes = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/{HAZARD_SLUG}/routes.pq", + edges_with_flows = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/{HAZARD_SLUG}/edges.gpq", + script: + "./allocate.py" + + +rule accumulate_route_costs_intact: + """ + For each route in the OD (source -> destination pair), lookup the edges of + the least cost route (the route taken) and sum those costs. Store alongside + value and volume of route. + + Test with: + snakemake -c1 -- results/flow_allocation/project-thailand/routes_with_costs.pq + """ + input: + routes = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/routes.pq", + edges_with_flows = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/edges.gpq", + output: + routes_with_costs = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/routes_with_costs.pq", + run: + from open_gira.routing import lookup_route_costs + + lookup_route_costs(input.routes, input.edges_with_flows).to_parquet(output.routes_with_costs) + + +rule accumulate_route_costs_degraded: + """ + For each route in the OD (source -> destination pair), lookup the edges of + the least cost route (the route taken) and sum those costs. Store alongside + value and volume of route. + + Test with: + snakemake -c1 -- results/flow_allocation/project-thailand/hazard-thai-floods-2011-JBA/routes_with_costs.pq + """ + input: + routes = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/{HAZARD_SLUG}/routes.pq", + edges_with_flows = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/{HAZARD_SLUG}/edges.gpq", + output: + routes_with_costs = "{OUTPUT_DIR}/flow_allocation/{PROJECT_SLUG}/{HAZARD_SLUG}/routes_with_costs.pq", + run: + from open_gira.routing import lookup_route_costs + + lookup_route_costs(input.routes, input.edges_with_flows).to_parquet(output.routes_with_costs) diff --git a/workflow/transport/multi-modal/multi_modal.smk b/workflow/transport/multi-modal/multi_modal.smk index ff3061de..d101298e 100644 --- a/workflow/transport/multi-modal/multi_modal.smk +++ b/workflow/transport/multi-modal/multi_modal.smk @@ -29,7 +29,7 @@ rule remove_edges_in_excess_of_threshold: edges = "{OUTPUT_DIR}/multi-modal_network/{PROJECT_SLUG}/edges.gpq", raster = "{OUTPUT_DIR}/hazard/{HAZARD_SLUG}.tif", output: - edges = "{OUTPUT_DIR}/multi-modal_network/{PROJECT_SLUG}/{HAZARD_SLUG}/chunks/{CHUNK_SLUG}/edges.gpq" + edges = temp("{OUTPUT_DIR}/multi-modal_network/{PROJECT_SLUG}/{HAZARD_SLUG}/chunks/{CHUNK_SLUG}/edges.gpq") run: import geopandas as gpd import numpy as np