Merge pull request #187 from nismod/feature/trade_flow_routing

Trade flow routing on multi-modal networks

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:

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
+```

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
+```

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
 

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
 

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",
+```

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"])