Import Libraries and Load the Graph#

We start by importing the necessary libraries and loading the graph we created in step 2.

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 gn

GDAL is not installed - OGR functionality not available

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, 4644, {'Wkt': 'LINESTRING (-21.9985141 64.1536718, -21.9984767 64.1537859, -21.9984354 64.1538681, -21.9983764 64.1539657, -21.9982927 64.154069)', 'id': 328, 'infra_type': 'residential', 'one_way': False, 'osm_id': '5013949', 'key': 'edge_328', 'length': 0.04569859554162313, '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, 4644, {'Wkt': 'LINESTRING (-21.9985141 64.1536718, -21.9984767 64.1537859, -21.9984354 64.1538681, -21.9983764 64.1539657, -21.9982927 64.154069)', 'id': 328, 'infra_type': 'residential', 'one_way': False, 'osm_id': '5013949', 'key': 'edge_328', 'length': 45.69859554162313, 'Type': 'legitimate', 'time': 8.225747197492165, '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, 3, 4, 5, 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: 1563 seconds, or 26 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]

We can visualize this in-line by just calling it:

poly

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

# the crs of churches
churches.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# view the pandana_snap doc string
gn.pandana_snap?

We 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

churches

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:

destinations = list(set(churches.NN))

destinations

[12032, 9540, 13257, 8370, 11157, 6557]

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_gdf

origin_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].NN

2302

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

OD

array([[1335.75240603,  635.47968749,  489.87692026,  348.95959429,
        1066.23418002,  817.05710822]])

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_df

It appears nearly all of the churches less than twenty minutes away.

Now you are up to speed on the basics!