In step 3 of this tutorial, we use our cleaned graph to create an Origin-Destination matrix (OD). Our setting remains Reykjavik, Iceland, as we look at travel times along the network to churches.
Import Libraries and Load the Graph¶
We start by importing the necessary libraries and loading the graph we created in step 2.
Import Libraries¶
We start by importing the necessary libraries:
osto set the working directorynetworkxto work with the graphgeopandasto work with the geospatial datapandasto work with the datapickleto load the graphshapelyto work with the geometryGOSTnetsto apply custom functions to the network
import os
import networkx as nx
import geopandas as gpd
import pandas as pd
import pickle as pkl
from shapely.geometry import Point
# import GOSTnet library
import GOSTnets as gnRead in the Graph¶
We set the path and read in the graph from the result of the cleaning process (Step 2)
pth = "./" # change this path to your working folder
G = pkl.load(
open(os.path.join(pth, "tutorial_outputs", r"iceland_network_clean.pickle"), "rb")
)Calculate Travel Times For Each Edge of the Graph¶
At this stage each edge in the network has a property called ‘length’. This was actually computed during Step 1 when the generateRoadsGDF function was run. The units of this length are in kilometres.
gn.example_edge(G)(1, 989, {'infra_type': 'footway', 'one_way': None, 'osm_id': 'b12d37b5-8a83-4b90-bd1d-8b37a86cb3ac', 'key': 'edge_365', 'length': 0.00253020699858194, 'Wkt': 'LINESTRING (-21.8314556 64.0944827, -21.8314875 64.0945006)', 'Type': 'legitimate'})
We want to convert length to time, so that we can conduct analysis on how long it takes to reach certain destinations.
We do this using the convert_network_to_time function.
We have used a factor of 1000, because the function is expecting meters, so we need to convert the units of kilometers to meters.
The convert_network_to_time function uses a default speed dictionary that assigns speed limits to OSM highway types.
However, it is possible to specify your own speed dictionary.
G_time = gn.convert_network_to_time(
G, distance_tag="length", road_col="infra_type", factor=1000
)We can now use the ‘time’ property for each edge to work out how long it takes to get from one node to another!
gn.example_edge(G_time, 1)(1, 989, {'infra_type': 'footway', 'one_way': None, 'osm_id': 'b12d37b5-8a83-4b90-bd1d-8b37a86cb3ac', 'key': 'edge_365', 'length': 2.53020699858194, 'Wkt': 'LINESTRING (-21.8314556 64.0944827, -21.8314875 64.0945006)', 'Type': 'legitimate', 'time': 0.45543725974474913, 'mode': 'drive'})
To do this for just one journey, we could call nx.shortest_path_length on any given origin or destination node.
Let’s list 10 of our nodes using this networkX function:
list(G_time.nodes)[:10][1, 2, 3, 4, 6, 7, 8, 9, 10, 11]Now we will calculate the travel time between the first and tenth node in the list.
A = list(G_time.nodes)[0] # first node in list
B = list(G_time.nodes)[10] # 10th node in list
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)
)The travel time between A and B is: 61 seconds, or 1 minutes!
Calculate Travel Times to Churches¶
In our example, we want to use our network for Reykjavik to work out the travel time to local churches.
Here, we import a shapefile for Reykjavik, and reproject it to WGS 84:
rek = gpd.read_file(os.path.join(pth, "tutorial_data", "rek2.shp"))
rek = rek.to_crs("epsg:4326")Next, we set a variable poly equal to just the geometry
poly = rek.geometry.iloc[0]With this in hand, we can read in a shapefile of destinations - here, the churches in Iceland. We use Shapely’s ‘within’ command to select just those in the Reykjavik area:
churches = gpd.read_file(os.path.join(pth, "tutorial_data", "churches.shp"))
churches = churches.loc[churches.within(poly)]In order to perform network analysis we want to know the closest network node to each church.
For this, we use the GOSTnets pandana_snap function to snap the church locations to the road network:
churches# view the pandana_snap doc string
gn.pandana_snap?Signature:
gn.pandana_snap(
G,
point_gdf,
source_crs='epsg:4326',
target_crs='epsg:4326',
add_dist_to_node_col=True,
time_it=False,
)
Docstring:
snaps points to a graph at very high speed
Parameters
----------
G : nx.Graph
a graph object, or the node geodataframe of a graph
point_gdf : gpd.GeoDataFrame
a geodataframe of points, in the same source crs as the geometry of the graph object
source_crs : str
The crs for the input G and input point_gdf in format 'epsg:32638'
target_crs : str
The measure crs how distances between points are calculated. The returned point GeoDataFrame's CRS does not get modified. The crs object in format 'epsg:32638'
add_dist_to_node_col : bool
return distance to nearest node in the units of the target_crs
Returns
-------
GeoDataFrame
returns a GeoDataFrame that is the same as the input point_gdf but adds a column containing the id of the nearest node in the graph, and the distance if add_dist_to_node_col is True
File: c:\wbg\anaconda3\envs\gostnets\lib\site-packages\gostnets\core.py
Type: functionWe want the nearest node distance (NN_dist) to be measured in meters, so that is why we include the target_crs parameter specifying the correct UTM zone.
churches = gn.pandana_snap_c(
G_time,
churches,
source_crs="epsg:4326",
target_crs="epsg:32627",
add_dist_to_node_col=True,
)As we can see from the NN_dist column, our church locations are very close to a node on the network in all cases
churchesWhen 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:
destinations = list(set(churches.NN))destinations[1289, 275, 517]Further Analysis¶
We would like to make an OD matrix where the origin is the cottage we are renting in the city, and the destinations are the churches in Reykjavik.
This will help us work out how many churches we can see today!
First, we need to create the origin, it has coordinates: 64.152215, -22.002099 (Lat,Lon), so we make a shapely.geometry.Point:
# A list with a single Shapely Point object is created with (x,y)
my_house = [Point(-22.002099, 64.152215)]Next, we load it into a GeoDataFrame and snap it to the network:
mini_gdf = gpd.GeoDataFrame(
{"geometry": my_house}, crs="epsg:4326", geometry="geometry", index=[1]
)mini_gdforigin_gdf = gn.pandana_snap_c(
G_time, mini_gdf, source_crs="epsg:4326", target_crs="epsg:32627"
)origin_gdf# This is the nearest node (NN)
origin_gdf.iloc[0].NNnp.int64(993)Now, we can calculate the OD matrix using the GOSTNets calculate_OD function.
Bear in mind it takes list objects as inputs:
origin = [origin_gdf.iloc[0].NN]OD = gn.calculate_OD(G_time, origin, destinations, fail_value=9999999)The OD matrix displays the time in seconds to reach each church
ODarray([[376.84044847, 572.71218581, 105.30502261]])We can use minutes as the measure by dividing every value in the OD Matrix by 60. Then we can convert the array nicely into a pandas Dataframe,
OD = OD / 60
OD_df = pd.DataFrame(OD, columns=destinations, index=origin)OD_dfIt appears nearly all of the churches less than twenty minutes away.
Now you are up to speed on the basics!