<jupyter_start><jupyter_text>Step 3: Using your GraphIn this third tutorial, we calculate an OD matrix as well as isochrones around a set of destinations. Our setting remains Reykjavik, Iceland, as we look at travel times along the network to churches.<jupyter_code>import os, sys
import geopandas as gpd
import pandas as pd
import importlib
import networkx as nx
from shapely.geometry import Point
sys.path.append(r'C:\Users\charl\Documents\GitHub\GOST_PublicGoods\GOSTNets\GOSTNets')
import GOSTnet as gn<jupyter_output>networkx version: 2.1 
osmnx version: 0.7.4<jupyter_text>First, we read in the result of the cleaning process (Step 2)<jupyter_code>pth = r'C:\Users\charl\Documents\T'
G = nx.read_gpickle(os.path.join(pth, 'tutorial_outputs', r'ISL_processed.pickle'))<jupyter_output><empty_output><jupyter_text>At this stage each edge in the network has a property called 'length'. This is added during the cleaning process. The units of this length are kilometres:<jupyter_code>gn.example_edge(G)<jupyter_output>(0, 11677, {'Wkt': 'LINESTRING (-21.7150886 64.16856079999999, -21.7150429 64.1684612, -21.7150343 64.1684424, -21.7150189 64.16841460000001, -21.7149683 64.1683788, -21.7149165 64.1683566, -21.714839 64.16832599999999, -21.7138473 64.1680006)', 'id': 5779, 'infra_type': 'residential', 'osm_id': '55759237', 'key': 'edge_5779', 'length': 0.09030580094429581, 'Type': 'legitimate'})<jupyter_text>We want to convert this to time, so that we can generate some isochrones. We do this below. We have used a factor of 1000 to convert the length back to metres, which is the assumed unit of measurement for this function call.<jupyter_code>G_time = gn.convert_network_to_time(G, distance_tag = 'length', road_col = 'infra_type', factor = 1000)<jupyter_output><empty_output><jupyter_text>We can now use the 'time' property for each edge to work out how long it takes to get from one node to another!<jupyter_code>gn.example_edge(G_time, 1)<jupyter_output>(0, 11677, {'Wkt': 'LINESTRING (-21.7150886 64.16856079999999, -21.7150429 64.1684612, -21.7150343 64.1684424, -21.7150189 64.16841460000001, -21.7149683 64.1683788, -21.7149165 64.1683566, -21.714839 64.16832599999999, -21.7138473 64.1680006)', 'id': 5779, 'infra_type': 'residential', 'osm_id': '55759237', 'key': 'edge_5779', 'length': 90.30580094429581, 'Type': 'legitimate', 'time': 16.255044169973246, 'mode': 'drive'})<jupyter_text>To do this for just one journey, we could call nx.shortest_path_length as per the below:<jupyter_code>A = 100 # node with id 100
B = 200 # node with id 200
travel_time = nx.shortest_path_length(G_time, A, B, weight = 'time')
print('The travel time between A and B is: %d seconds, or %d minutes!' % (travel_time, travel_time / 60))<jupyter_output>The travel time between A and B is: 4464 seconds, or 74 minutes!<jupyter_text>In our example, we want to use our network for Reykjavik to work out the travel time to local churches.Here, I import a shapefile for Reykjavik, and reproject it to WGS 84:<jupyter_code>rek = gpd.read_file(os.path.join(pth, 'tutorial_data', 'rek2.shp'))
rek = rek.to_crs({'init':'epsg:4326'})<jupyter_output><empty_output><jupyter_text>Next, I set a variable poly equal to just the geometry<jupyter_code>poly = rek.geometry.iloc[0]<jupyter_output><empty_output><jupyter_text>We can visualize this in-line by just calling it:<jupyter_code>poly<jupyter_output><empty_output><jupyter_text>With this in hand, I can read in a shapefile of destinations - here, the churches in Iceland. I use shapely's 'within' command to select just those in the Reykjavik area:<jupyter_code>churches = gpd.read_file(os.path.join(pth, 'tutorial_data', 'churches.shp'))
churches = churches.loc[churches.within(poly)]<jupyter_output><empty_output><jupyter_text>I want to know the closest network node to each church. For this, we use pandana snap to snap the church locations to the road network:<jupyter_code>churches = gn.pandana_snap(G_time, churches, add_dist_to_node_col = True)<jupyter_output><empty_output><jupyter_text>As we can see from the NN_dist column, our church locations are very close to a node on the network in all cases<jupyter_code>churches<jupyter_output><empty_output><jupyter_text>When calculating an OD-Matrix, we can only use the node IDs as inputs. So, we convert this column of our dataframe over to a list of unique values:<jupyter_code>destinations = list(set(list(churches.NN)))<jupyter_output><empty_output><jupyter_text>At this point, we have a choice to make. If we know our origins, we can use gn.calculate_OD() to calculate the OD matrix from each origin to all destinations. If we do not have fixed origins in mind, we could also visualize the isochrones around the network itself, with gn.make_iso_polys. Here, I generate a polygon geodataframe of the isochrones around the churches, at the 5, 10 and 15 minute travel time intervals. I have chosen to buffer each edge by 10m, and each node by 25m, to make the visualization. I am measuring in the nearest UTM zone with my choice of measure_crs. N.B. Please note - this step takes a while. If you want to speed it up, reduce the trip times options down somewhat. Here, I am generating 5 minute, 10 minute and 15 minute drive time isochrones.<jupyter_code>isopolys = gn.make_iso_polys(G_time, 
                  destinations, 
                  trip_times = [5*60, 10*60, 15*60], 
                  edge_buff=10, 
                  node_buff=25, 
                  infill=False, 
                  weight = 'time', 
                  measure_crs = {'init': 'epsg:32627'})<jupyter_output><empty_output><jupyter_text>I can cheaply visualize them here with geopandas' plot method:<jupyter_code>isopolys.plot()<jupyter_output><empty_output><jupyter_text>Or, save them down to visualize in QGIS, after resetting the crs to WGS 84 for ease of use:<jupyter_code>isopolys = isopolys.to_crs({'init': 'epsg:4326'})
isopolys.to_csv(os.path.join(pth, 'tutorial_outputs', 'church_isopolys.csv'))<jupyter_output><empty_output><jupyter_text>Note, we could, with a little extra work transform the pandas dataframe into a geopandas GeoDataFrame, and then send this to file as a geoJSON or shapefile. Next, I decide I would like to make an OD-matrix where the origin is the cottage I am renting in the city, and the destinations are the churches in Reykjavik. This will help me to wok out how many churches I can see today!. First, I need to create my origin. It has coordinates: 64.152215, -22.002099, so I make a point of this:<jupyter_code>my_house = [Point(64.152215, -22.002099)]<jupyter_output><empty_output><jupyter_text>Next, I load it into a geodataframe and snap it to the network:<jupyter_code>mini_gdf = gpd.GeoDataFrame({'point':my_house}, crs = {'init':'epsg:4326'}, geometry = 'point', index = [1])
mini_gdf
origin_gdf = gn.pandana_snap(G_time, mini_gdf)
origin_gdf<jupyter_output><empty_output><jupyter_text>Now, I can calcuate my OD-matrix using GN.calculate_OD(). Bear in mind it takes list objects as inputs:<jupyter_code>origin = [2129]
OD = gn.calculate_OD(G_time, origin, destinations, fail_value = 9999999)<jupyter_output><empty_output><jupyter_text>I can convert this nicely into a pandas Dataframe, using minutes as the measure as so:<jupyter_code>OD = OD / 60
OD_df = pd.DataFrame(OD, columns = destinations, index = origin)
OD_df<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09809782608695652, 'fraction_numerical': 0.07567934782608696, 'mean_word_length': 5.423211169284468, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 77, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 5, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '257474', 'license': 'MIT', 'n_tokens_mistral': 2940, 'n_tokens_neox': 2474, 'n_words': 822, 'path': 'Tutorial/Step 3 - Using your Graph.ipynb', 'repo_name': 'd3netxer/GOSTnets', 'size': '7360'}