<jupyter_start><jupyter_text>Prep B - Road Closure MergeIt is well known that conflict leads to road damage and closure - with various factions setting up road blocks, and munitions damaging road surfaces, sometimes in major ways. We received from the UN Logistics Cluster (via Dr. <NAME>, World Bank) the latest road closure information for Yemen, as well as 4 file versions for previous time periods. However, this information is limited to only the largest grade roads in the area (main roads, paved) whereas our own network is derived from Open Street Map which has been purposefully improved over the region. As such, we want to spatially join the information from this road shapefile over to our own network. This is what this script accomplishes.  We start by importing the usual suspects. If you have installed GOSTnets via conda, you will not need to append to the system path as we have done below.<jupyter_code>import pandas as pd
import os, sys
sys.path.append(r'C:\Users\charl\Documents\GitHub\GOST_PublicGoods\GOSTNets\GOSTNets')
sys.path.append(r'C:\Users\charl\Documents\GitHub\GOST')
import GOSTnet as gn
import importlib
import geopandas as gpd
from shapely.geometry import Point
import rasterio
from rasterio import features
from shapely.wkt import loads
import numpy as np
import networkx as nx
import time<jupyter_output>peartree version: 0.6.1 
networkx version: 2.3 
matplotlib version: 3.0.3 
osmnx version: 0.9<jupyter_text>Load Road Network We set up net_graph as the pickled network object. We are using the Salted network with time attributes.<jupyter_code>bp = r'C:\Users\charl\Documents\GOST\Yemen'
net_graph = nx.read_gpickle(os.path.join(bp, 'YEM', 'Round 3', 'G_salty_time.pickle'))<jupyter_output><empty_output><jupyter_text>Here we add a new attribute - ID - for each edge, and then generate a GeoDataFrame from the graph's edges using GOSTnets<jupyter_code>counter = 1
for u, v, data in net_graph.edges(data = True):
    data['ID'] = counter
    counter+=1
net_graph_edge_gdf = gn.edge_gdf_from_graph(net_graph)<jupyter_output><empty_output><jupyter_text>Match on Conflict information to Graph In this section, we set up for a large loop, where we iteratively read in and spatially match on the information from each shapefile on to our main graph object.<jupyter_code># Create a list of the files and their dates to iterate through
settings = [
{'date':'November1st',
'f_name':'Access_20181101.shp'},
{'date':'November8th',
'f_name':'Access_20181108.shp'},
{'date':'November14th',
  'f_name':r'Access_20181114.shp'}, 
{'date':'November25th',
'f_name':'Access_20181125.shp'},
{'date':'December17th',
'f_name':'Access_20181217.shp'},
{'date':'January24th',
'f_name':r'Access_20190124.shp'}]<jupyter_output><empty_output><jupyter_text>Below is the main loop. Pay attention to in-line commentary to see how it works.<jupyter_code># set up a spatial index using the GeoPandas .sindex method
smol = net_graph_edge_gdf
spatial_index = smol.sindex

# iterate through our settings list....
for setting in settings:
    
    # Extract from the setting dictionary the filename and the date it relates to
    f_name = setting['f_name']
    date = setting['date']

    # read in the UN logistics cluster shapefile for that date using f_name
    conflict = gpd.read_file(os.path.join(r'C:\Users\charl\Documents\GOST\Yemen\Conflicts\Road Networks (UN Logistics Cluster)',f_name))
    
    # pick out only the roads where their status is either restricted or closed. The rest of the roads are in full working operation
    conflict = conflict.loc[conflict.status.isin(['Restricted','Closed'])]
    
    # Reproject the conflict-affected roads to UTM zone 
    conflict = conflict.to_crs({'init':'epsg:32638'})
    
    # Buffer these roads by 50 meters
    conflict['buffered'] = conflict.buffer(50)
    
    # reset the GeoDataFrame's geometry attribute to the buffered road geometry
    conflict = conflict.set_geometry('buffered')
    
    # reset the projection to WGS 84
    conflict = conflict.to_crs({'init':'epsg:4326'})
    
    # split out the restricted and the closed roads into their only GeoDataFrames
    conflict_restricted = conflict.copy()
    conflict_restricted = conflict_restricted.loc[conflict_restricted.status.isin(['Restricted'])]
    conflict_closed = conflict.copy()
    conflict_closed = conflict_closed.loc[conflict_closed.status.isin(['Closed'])]
    
    # now, we iterate through each closed road, and use the spatial index to quickly identify 
    # intersecting graph edges with the buffered closed roads. This may capture a few more roads that strictly necessary, but 
    # that is a price we are willing to pay, and probably reflects the 'on the ground' situation (remember no edge is longer than 2km anyway)
    close = []
    for index, row in conflict_closed.iterrows():
        polygon = row.buffered
        possible_matches_index = list(spatial_index.intersection((polygon.bounds)))
        possible_matches = smol.iloc[possible_matches_index]
        precise_matches = possible_matches[possible_matches.intersects(polygon)]
        i = list(precise_matches.ID)
        if len(i) !=0:
            close.append(i)
    close = [item for sublist in close for item in sublist]
    
    # we perform the same operation now on restricted roads
    restrict = []
    for index, row in conflict_restricted.iterrows():
        polygon = row.buffered
        possible_matches_index = list(spatial_index.intersection((polygon.bounds)))
        possible_matches = smol.iloc[possible_matches_index]
        precise_matches = possible_matches[possible_matches.intersects(polygon)]
        i = list(precise_matches.ID)
        if len(i) !=0:
            restrict.append(i)
    restrict = [item for sublist in restrict for item in sublist]
    
    # Closed roads are a 'higher order' problem that restricted roads. As such, should the buffer nature 
    # of the process have assigned a road both to the closed AND restricted buckets, we remove it from the restricted bucket. 
    # this shouldn't happen too often, but is entirely possible where a closed road meets a restricted road. 
    for i in restrict:
        if i in close:
            restrict.remove(i)
    
    # now, we iterate through the edges, and check to see whether their IDs crop up in the list items 'close' and 'restrict' - 
    # and if so, we adjust their time accordingly. 
    for u, v, data in net_graph.edges(data = True):
        
        # if a road is on the restricted list, we double the travel time. 
        if data['ID'] in restrict:
            data['time_%s' % date] = data['time'] * 2
            data['MOD_%s' % date] = 'restricted'
            
        # if a road is closed, we make the travel time very, very large. 
        elif data['ID'] in close:
            data['time_%s' % date] = 99999999
            data['MOD_%s' % date] = 'closed'
            
        # otherwise, we pass
        else:
            data['time_%s' % date] = data['time']
            data['MOD_%s' % date] = 'normal'<jupyter_output><empty_output><jupyter_text>Save Down We save the resulting graph down as 'conflict adjusted'. Note, we did not edit the 'time' property, but set a NEW property for the travel time as of the date of the conflict info - so the original time can still be accessed from the conflcit adjusted Graph object.<jupyter_code>gn.save(net_graph, 'G_salty_time_conflict_adj', os.path.join(bp, 'YEM', 'Round 3'), nodes = False, edges = False)<jupyter_output><empty_output><jupyter_text>Snap on  At this point, we also snap on our facility file and origin files to the conflcit adjusted graph. Although we did not change any nodes, its always nice to be sure you are working with the latest and greatest file version :-)<jupyter_code>dfiles = ['HeRAMS 2018 April.csv']

dpath = r'C:\Users\charl\Documents\GOST\Yemen\facility_files'
wpath = r'C:\Users\charl\Documents\GOST\Yemen\graphtool'

for dfile in dfiles:
    # Read in
    dest_df = pd.read_csv(os.path.join(os.path.join(dpath, dfile)), encoding = "ISO-8859-1")
    
    # Ensure coordinates are floats
    dest_df.Longitude = dest_df.Longitude.astype(float)
    dest_df.Latitude = dest_df.Latitude.astype(float)
    
    # Drop entries with no coordinates    
    dest_df2 = dest_df.copy()
    print(len(dest_df2))
    dest_df2 = dest_df2.loc[(dest_df2.Longitude != 0)]
    dest_df2 = dest_df2.loc[(dest_df2.Longitude != None)]
    dest_df2 = dest_df2.loc[(dest_df2.Longitude <= 60)]
    dest_df2 = dest_df2.loc[(dest_df2.Longitude >= 35)]
    dest_df2 = dest_df2.loc[(dest_df2.Latitude <= 30)]
    dest_df2 = dest_df2.loc[(dest_df2.Latitude >= 5)]
    print(len(dest_df2))
    dest_df = dest_df2
    
    # Generate Geometries
    dest_df['geometry'] = list(zip(dest_df.Longitude, dest_df.Latitude))
    dest_df['geometry'] = dest_df['geometry'].apply(Point)
    dest_df = gpd.GeoDataFrame(dest_df, geometry = 'geometry', crs = {'init':'espg:4326'})
    
    # Perform snap
    time.ctime()
    start = time.time()
    df = gn.pandana_snap(net_graph, dest_df, 'epsg:4326','epsg:32638', add_dist_to_node_col = True)
    
    # Save to file
    df.to_csv(os.path.join(wpath, dfile.replace('.csv', '_snapped.csv')))
    df.to_csv(os.path.join(dpath, dfile.replace('.csv', '_snapped.csv')))

    print('time elapsed: %d seconds' % (time.time() - start))
ofile = r'origins_1km.csv'
ofiles = [ofile] 
opath = r'C:\Users\charl\Documents\GOST\Yemen\origins'
for ofile in ofiles:
    
    # Read in
    dest_df = pd.read_csv(os.path.join(os.path.join(opath, ofile)), encoding = "ISO-8859-1")
    dest_df['geometry'] = dest_df['geometry'].apply(loads)
    dest_df = gpd.GeoDataFrame(dest_df, geometry = 'geometry', crs = {'init':'espg:4326'})
    
    # Perform snap
    print('Beginning snap')
    
    time.ctime()
    start = time.time()
    df = gn.pandana_snap(net_graph, dest_df, 'epsg:4326','epsg:32638', add_dist_to_node_col = True)
    
    # Save to file
    df.to_csv(os.path.join(wpath, ofile.replace('.csv', '_snapped.csv')))
    df.to_csv(os.path.join(opath, ofile.replace('.csv', '_snapped.csv')))

    print('Time elapsed: %d seconds' % (time.time() - start))<jupyter_output>Beginning snap
Time elapsed: 166 seconds