<jupyter_start><jupyter_text>Data preparation and infrastructure exposure to floodingThis notebook forms the basis of "Hands-On 5" in the CCG course.1. Extract infrastructure data from OpenStreetMap2. Extract flood hazard data from Aqueduct3. Intersect floods with roads to calculate exposure4. Open QGIS to look at the data  Activity 1: Extract infrastructure data Step 1) On your desktop, create a folder called `ghana_tutorial` Step 2) Create a variable to store the folder locationIn the cell below, type in the path to your desktop, by changing NAME to match your username as shown on your computer<jupyter_code># edit this if using a Mac (otherwise delete)
data_folder = "/Users/NAME/Desktop/ghana_tutorial"

# edit this if using Windows (otherwise delete)
data_folder = "C:\\Users\\NAME\\Desktop\\ghana_tutorial"

# delete this line
data_folder = "../data"<jupyter_output><empty_output><jupyter_text>Step 3) Load Python libraries<jupyter_code># 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

# 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 
# 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<jupyter_output><empty_output><jupyter_text>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`  Step 6) Load the road dataset you've just downloaded<jupyter_code>roads = gpd.read_file(
    os.path.join(data_folder, 'ghana-latest-free.shp', 'gis_osm_roads_free_1.shp'))<jupyter_output><empty_output><jupyter_text>Step 7) To take a look at the data and the attribute table fill in and run the next two cells<jupyter_code>roads
roads.fclass.unique()<jupyter_output><empty_output><jupyter_text>Step 8) Next we want to make a couple of changes to the data Filter out minor and residential roads, tracks and paths.<jupyter_code># 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',
    })<jupyter_output><empty_output><jupyter_text>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<jupyter_code>road_network = snkit.Network(edges=roads)
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<jupyter_output><empty_output><jupyter_text>Calculate the length of each road segment in meters<jupyter_code>geod = Geod(ellps='WGS84')
roads['length_m'] = roads.geometry.apply(geod.geometry_length)
roads.tail()
roads.crs = {'init': 'epsg:4326'}
road_nodes.crs = {'init': 'epsg:4326'}<jupyter_output><empty_output><jupyter_text>Step 9) Save the pre-processed dataset<jupyter_code>roads.to_file(
    os.path.join(data_folder, 'GHA_OSM_roads.gpkg'),
    layer='edges', 
    driver="GPKG")
road_nodes.to_file(
    os.path.join(data_folder, 'GHA_OSM_roads.gpkg'),
    layer='nodes', 
    driver="GPKG")<jupyter_output><empty_output><jupyter_text>Activity 2: Extract and polygonise hazard data  Step 1) Download flood hazard data from Aqueduct The full [Aqueduct dataset](https://www.wri.org/resources/data-sets/aqueduct-floods-hazard-maps) is available to download openly. There are some scripts and summary of the data you may find useful at [nismod/aqueduct](https://github.com/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.For later tutorials, we provide a preprocessed set of hazard polygons for the Ghana example. The next steps show how to clip a region out of the global dataset and polygonise it, in case you wish to reproduce this analysis in another part of the world.For now, we suggest downloading [inunriver_historical_000000000WATCH_1980_rp00100.tif](http://wri-projects.s3.amazonaws.com/AqueductFloodTool/download/v2/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.  Step 2) Run the code below to polygonise the tif filesThis converts the flood maps from *tiff files (raster data)* into *shape files (vector data)*. It will take a little time to run.<jupyter_code>xmin = "-3.262509"
ymin = "4.737128"
xmax = "1.187968"
ymax = "11.162937"

for root, dirs, files in os.walk(data_folder, 'flood_layer'): 
    print("Looking in", root) 
    for file in sorted(files): 
        if file.endswith(".tif") and not file.endswith("p.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, f"{stem}_clip.tif")
            try:
                os.remove(clip_file)
            except FileNotFoundError:
                pass
            p = subprocess.run([
                "gdalwarp", "-te", xmin, ymin, xmax, ymax, input_file, clip_file],
                capture_output=True)
            print(p.stdout.decode('utf8'))
            print(p.stderr.decode('utf8'))
            print(clip_file)

            # Create vector outline of raster areas
            # note that this rounds the floating-point values of flood depth from 
            # the raster to the nearest integer in the vector outlines
            polygons_file = os.path.join(root, f"{stem}.gpkg") 
            try:
                os.remove(polygons_file)
            except FileNotFoundError:
                pass
            p = subprocess.run([
                "gdal_polygonize.py", clip_file,'-q','-f', 'GPKG', polygons_file],
                capture_output=True)
            print(p.stdout.decode('utf8'))
            print(p.stderr.decode('utf8'))  
            print(polygons_file)<jupyter_output><empty_output><jupyter_text>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<jupyter_code>flood_path = os.path.join(
    data_folder, 
    'flood_layer', 
    'inunriver_historical_000000000WATCH_1980_rp00100.gpkg')

output_path = os.path.join(
    data_folder, 
    'results', 
    'inunriver_historical_000000000WATCH_1980_rp00100_exposure.gpkg')

flood = gpd.read_file(flood_path) \
    .rename(columns={'DN': 'depth_m'})
flood = flood[flood.depth_m > 0]<jupyter_output><empty_output><jupyter_text>Step 2) Run the intersection<jupyter_code>flood_intersections = gpd.overlay(GHA_OSM_roads, flood, how='intersection')
flood_intersections<jupyter_output><empty_output><jupyter_text>Calculate the exposed length<jupyter_code>geod = Geod(ellps='WGS84')
flood_intersections['flood_length_m'] = flood_intersections.geometry.apply(geod.geometry_length)
flood_intersections.tail()<jupyter_output><empty_output><jupyter_text>Calculate the proportion of roads in our dataset which are exposed to >=1m flood depths in this scenario<jupyter_code>exposed_length = flood_intersections.flood_length_m.sum()
exposed_length
all_roads_in_dataset_length = GHA_OSM_roads.length_m.sum()
all_roads_in_dataset_length
proportion = exposed_length / all_roads_in_dataset_length
proportion
f"{proportion:.0%} of roads in this dataset are exposed to flood depths of >= 1m in a historical 1-in-100 year flood"<jupyter_output><empty_output><jupyter_text>Save to file (with spatial data)<jupyter_code>flood_intersections.to_file(output_path, driver="GPKG")
flood_intersections<jupyter_output><empty_output><jupyter_text>Save to CSV (without spatial data)<jupyter_code>flood_intersections.drop(columns='geometry').to_csv(output_path.replace(".gpkg", ".csv"))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09499117387466902, 'fraction_numerical': 0.017541924095322153, 'mean_word_length': 4.326086956521739, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 64, 'https://': 8, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '521141', 'license': 'MIT', 'n_tokens_mistral': 3072, 'n_tokens_neox': 2821, 'n_words': 817, 'path': 'tutorials/01-data-preparation-ghana.ipynb', 'repo_name': 'nismod/snail', 'size': '9064'}