Data preparation and infrastructure exposure to flooding

This notebook forms the basis of “Hands-On 5” in the CCG course.

1. Extract infrastructure data from OpenStreetMap

2. Extract flood hazard data from Aqueduct

3. Intersect floods with roads to calculate exposure

4. Open QGIS to look at the data

# The os and subprocess modules are built into Python
# see https://docs.python.org/3/library/os.html
import os

# see https://docs.python.org/3/library/subprocess.html
import subprocess

# see https://docs.python.org/3/library/time.html
import time

# see https://docs.python.org/3/library/pathlib.html
from pathlib import Path

Activity 1: Extract infrastructure data

Step 1) In the parent directory of this package, create a folder called ghana_tutorial

Step 2) Create a variable to store the folder location

dir = (
    Path(os.getcwd()).resolve().parents[3]
)  # get parent directory of snail package
new_folder = os.path.join(dir, "ghana_tutorial")
if not os.path.exists(new_folder):
    os.makedirs(new_folder)

# edit this if using a Mac (otherwise delete)
data_folder = Path(new_folder)  # Path("YOUR_PATH/ghana_tutorial")

# edit this if using Windows (otherwise delete)
# data_folder = Path("C:YOUR_PATH/ghana_tutorial")

# delete this line
# data_folder = Path("../data")

Step 3) Load Python libraries

# Pandas and GeoPandas are libraries for working with datasets
# see https://geopandas.org/
import geopandas as gpd

gpd._compat.USE_PYGEOS = False
# see https://pandas.pydata.org/
import pandas as pd

# This package interacts with a risk data extract service, also accessible at
# https://global.infrastructureresilience.org/downloads
import irv_autopkg_client

# We'll use snail to intersect roads with flooding
import snail.intersection
import snail.io

# snkit helps generate connected networks from lines and nodes
# see https://snkit.readthedocs.io/
import snkit
import snkit.network

# PyPROJ is a library for working with geographic projections
# see https://pyproj4.github.io/
from pyproj import Geod

from matplotlib import pyplot as plt

from urllib.request import urlretrieve
import zipfile

Step 4) and 5) Download and save data

Download the ghana-latest-free.shp.zip dataset from http://download.geofabrik.de/africa/ghana.html, extract the zip folder and save the extracted folder within your new folder ghana_tutorial

download_url = "https://download.geofabrik.de/africa/ghana-latest-free.shp.zip"
file_path = os.path.join(data_folder, "ghana-osm.zip")
# check if directory already exists
if not os.path.exists(Path(os.path.splitext(file_path)[0])):
    # check if data already exists in this directory
    # if len(os.listdir(Path(os.path.splitext(file_path)[0]))) == 0:
    urlretrieve(download_url, file_path)  # note: this can take a few minutes
    with zipfile.ZipFile(file_path, "r") as zip_ref:
        zip_ref.extractall(os.path.splitext(file_path)[0])
else:
    print("data already exists")

Step 6) Load the road dataset you’ve just downloaded

roads = gpd.read_file(
    os.path.join(os.path.splitext(file_path)[0], "gis_osm_roads_free_1.shp")
)

Step 7) To take a look at the data and the attribute table fill in and run the next two cells

roads.head(5)

roads.fclass.unique()

Step 8) Next we want to make a couple of changes to the data

Filter out minor and residential roads, tracks and paths.

# Keep only the specified columns
roads = roads[["osm_id", "fclass", "name", "geometry"]]
# Keep only the roads whose "fclass" is in the list
roads = roads[
    roads.fclass.isin(
        [
            "motorway",
            "motorway_link",
            "trunk",
            "trunk_link",
            "primary",
            "primary_link",
            "secondary",
            "secondary_link",
            "tertiary",
            "tertiary_link",
        ]
    )
]
# Rename some columns
roads = roads.rename(
    columns={
        "fclass": "road_type",
    }
)

Create topological network information - this adds information that will let us find routes over the road network. - add nodes at the start and end of each road segment - split roads at junctions, so each segment goes from junction to junction - add ids to each node and edge, and add from_id and to_id to each edge

road_network = snkit.Network(edges=roads)

type(road_network)

with_endpoints = snkit.network.add_endpoints(road_network)
split_edges = snkit.network.split_edges_at_nodes(with_endpoints)
with_ids = snkit.network.add_ids(
    split_edges, id_col="id", edge_prefix="roade", node_prefix="roadn"
)
connected = snkit.network.add_topology(with_ids)
roads = connected.edges
road_nodes = connected.nodes

Calculate the length of each road segment in meters

geod = Geod(ellps="WGS84")
roads["length_m"] = roads.geometry.apply(geod.geometry_length)

roads.tail(5)

roads.set_crs(4326, inplace=True)
road_nodes.set_crs(4326, inplace=True)
road_nodes.crs

main_roads = roads[
    roads["road_type"].isin(
        [
            "trunk",
            "secondary",
        ]
    )
]

f, ax = plt.subplots()

main_roads.plot(
    ax=ax,
    alpha=1,
    linewidth=0.5,
)

ax.grid()
ax.set_title("Main roads of Ghana")
ax.set_xlabel("Longitude [deg]")
ax.set_ylabel("Latitude [deg]")

Step 9) Save the pre-processed dataset

roads.to_file(
    data_folder / "GHA_OSM_roads.gpkg",
    layer="edges",
    driver="GPKG",
)
road_nodes.to_file(
    data_folder / "GHA_OSM_roads.gpkg",
    layer="nodes",
    driver="GPKG",
)

Activity 2: Extract hazard data

The full Aqueduct dataset is available to download openly.

Country-level extracts are available through the Global Systemic Risk Assessment Tool (G-SRAT). This section uses that service to download an extract for Ghana.

Alternative: download flood hazard data from Aqueduct

The full Aqueduct dataset is available to download. There are some scripts and summary of the data you may find useful at nismod/aqueduct.

There are almost 700 files in the full Aqueduct dataset, of up to around 100MB each, so we don’t recommend downloading all of them unless you intend to do further analysis.

The next steps show how to clip a region out of the global dataset, in case you prefer to work from the original global Aqueduct files.

To follow this step, we suggest downloading inunriver_historical_000000000WATCH_1980_rp00100.tif to work through the next steps. Save the downloaded file in a new folder titled flood_layer under your data_folder.

country_iso = "gha"

Create a client to connect to the data API:

client = irv_autopkg_client.Client()

# client.dataset_list()

# client.dataset("wri_aqueduct.version_2")

client.extract_download(
    country_iso,
    data_folder / "flood_layer",
    # there may be other datasets available, but only download the following
    dataset_filter=["wri_aqueduct.version_2"],
    overwrite=True,
)

xmin = "-3.262509"
ymin = "4.737128"
xmax = "1.187968"
ymax = "11.162937"

for root, dirs, files in os.walk(os.path.join(data_folder, "flood_layer")):
    print("Looking in", root)
    for file_ in sorted(files):
        if file_.endswith(".tif") and not file_.endswith(
            f"-{country_iso}.tif"
        ):
            print("Found tif file", file_)
            stem = file_[:-4]
            input_file = os.path.join(root, file_)

            # Clip file to bounds
            clip_file = os.path.join(
                root,
                "gha",
                "wri_aqueduct_version_2",
                f"{stem}-{country_iso}.tif",
            )
            try:
                os.remove(clip_file)
            except FileNotFoundError:
                pass
            cmd = [
                "gdalwarp",
                "-te",
                xmin,
                ymin,
                xmax,
                ymax,
                input_file,
                clip_file,
            ]
            print(cmd)
            p = subprocess.run(cmd, capture_output=True)
            print(p.stdout.decode("utf8"))
            print(p.stderr.decode("utf8"))
            print(clip_file)

Activity 3: Intersect hazard

Let us now intersect the hazard and the roads, starting with one hazard initially so we save time.

Step 1) Specify your input and output path as well as the name of the intersection

flood_path = Path(
    data_folder,
    "flood_layer",
    "gha",
    "wri_aqueduct_version_2",
    "wri_aqueduct-version_2-inunriver_historical_000000000WATCH_1980_rp00100-gha.tif",
)

output_path = Path(
    data_folder,
    "results",
    "inunriver_historical_000000000WATCH_1980_rp00100__roads_exposure.gpkg",
)

Read in pre-processed road edges, as created earlier.

roads = gpd.read_file(data_folder / "GHA_OSM_roads.gpkg", layer="edges")

Step 2) Run the intersection

grid, bands = snail.io.read_raster_metadata(flood_path)

prepared = snail.intersection.prepare_linestrings(roads)
flood_intersections = snail.intersection.split_linestrings(prepared, grid)
flood_intersections = snail.intersection.apply_indices(
    flood_intersections, grid
)
flood_data = snail.io.read_raster_band_data(flood_path)
flood_intersections["inunriver__epoch_historical__rcp_baseline__rp_100"] = (
    snail.intersection.get_raster_values_for_splits(
        flood_intersections, flood_data
    )
)

Calculate the exposed length

geod = Geod(ellps="WGS84")
flood_intersections["flood_length_m"] = flood_intersections.geometry.apply(
    geod.geometry_length
)

flood_intersections.tail(2)

Calculate the proportion of roads in our dataset which are exposed to >=1m flood depths in this scenario

exposed_1m = flood_intersections[
    flood_intersections.inunriver__epoch_historical__rcp_baseline__rp_100 >= 1
]
exposed_length_km = exposed_1m.flood_length_m.sum() * 1e-3
exposed_length_km

all_roads_in_dataset_length_km = roads.length_m.sum() * 1e-3
all_roads_in_dataset_length_km

proportion = exposed_length_km / all_roads_in_dataset_length_km
proportion

f"{proportion:.1%} of roads in this dataset are exposed to flood depths of >= 1m in a historical 1-in-100 year flood"

output_path.parent.mkdir(parents=True, exist_ok=True)

Save to file (with spatial data)

flood_intersections.to_file(output_path, driver="GPKG")

Save to CSV (without spatial data)

flood_intersections.drop(columns="geometry").to_csv(
    output_path.parent / output_path.name.replace(".gpkg", ".csv")
)