(This tutorial was also presented during JuliaCon 2020)
Any transportation network can be represented as a complex directed graph where vertices are spread an Euclidean space. The library provides a bridging functionality between real world spatial data available in the OpenStreetMap project and LightGraphs.jl and makes it possible to run real-life-sized experiment on transportation networks along with various visualizations.
A transportation system or even an entire city can be represented as a complex directed graph embedded in an Euclidean space. Such graph can model real world in 1:1 scale and be used to perform various numerical experiments. The OpenStreetMapX.jl package makes it possible to load the data from the OpenStreetMap.org project and processes such graphs with Julia. The package is using LightGraphs.jl to represent the directed graph structure object along with meta related to spatial information. For the graph vizualisation we will use Python's folium package for providing an interactive vizualization of results.
This notebook can be downloaded in Jupyter ipynb format here (right-click to download).
The source code uses a torontoF.osm file which you can also download here.
Firstly, let us start by installing the required Julia and Python packages. (If the packages are required uncomment the code below)
#using Pkg
#pkg"add OpenStreetMapX, Parameters, LightGraphs, PyCall, Conda, Colors"
#using Conda
#Conda.runconda(`install folium -c conda-forge --yes`)
Now we are ready to load all required packages
using Random
using Parameters
using OpenStreetMapX
using LightGraphs
using PyCall
using Colors
const flm = pyimport("folium");
pwd() # use to check in which foler you are, the folder should contain the torontoF.osm file
We start by loading the file. Note that we trim the map file to have a fully connected road network.
The file was downloaded from the OpenStreetMap project web page. The steps included: (1) Select some area on map; (2) Click the "Export" button at the top of the page; (3) Click "Manually select a different area" to select the central Toronto area; (4) Press the "Export" button on the left (note that sometimes the Export link does not work - in this case click one of the Overpass API link below the Export button).
m = get_map_data("torontoF.osm", use_cache=false, trim_to_connected_graph=true );
fieldnames(MapData)
The central element of the MapData object is a LightGraphs' representation of the road network
m.g
The LightGraph graph is represented by nodes, each node id can be directly mapped to OpenStreetMap.
m.n
Furhther those nodes can be presented as geographic coordinates
m.nodes
Let's take 10 cars an let them drive between randomly selected pairs of points
Random.seed!(0)
node_ids = collect(keys(m.nodes))
routes = Vector{Vector{Int}}()
for k in 1:10
a,b = rand(1:nv(m.g), 2)
route, route_time = OpenStreetMapX.shortest_route(m,m.n[a],m.n[b])
push!(routes, route)
end
Now we will plot those routes
fm = flm.Map()
colrs = distinguishable_colors(10, [RGB(1,0.6,0.5)])
for k=1:10
locs = [LLA(m.nodes[n],m.bounds) for n in routes[k]]
info = "The route of Car $k\n<BR>"*
"Length: $(length(routes[k])) nodes\n<br>" *
"From: $(routes[k][1]) $(round.((locs[1].lat, locs[1].lon),digits=4))\n<br>" *
"To: $(routes[k][end]) $(round.((locs[end].lat, locs[end].lon),digits=4))"
flm.PolyLine(
[(loc.lat, loc.lon) for loc in locs ],
popup=info,
tooltip=info,
color="#$(hex(colrs[k]))"
).add_to(fm)
end
MAP_BOUNDS = [(m.bounds.min_y,m.bounds.min_x),(m.bounds.max_y,m.bounds.max_x)]
flm.Rectangle(MAP_BOUNDS, color="black",weight=6).add_to(fm)
fm.fit_bounds(MAP_BOUNDS)
fm
Let now us try to built a simple simulation model of a pandemic in the city.
We start with a definition of an Agent
@with_kw mutable struct Agent
id::Int
last_node::Int=-1
current_node::Int
infected::Bool = false
infection_time::Int = -1
total_route_len::Int = 0
end
Once the agent is defined let us define the enviroment where they can move around.
@with_kw struct Simulation
agents::Vector{Agent}
m::OpenStreetMapX.MapData
nodes_infected::Vector{Set{Int}}
nodes_agents::Vector{Set{Int}}
infected_agents_count::Vector{Int}
end
function Simulation(m::OpenStreetMapX.MapData, N=100)
vv = size(m.g)[1] # number of vertices
nodes_infected = [ Set{Int}() for s in 1:vv ]
nodes_agents = [ Set{Int}() for s in 1:vv ]
agents = Agent[]
for i in 1:N
a = Agent(id=i, current_node=rand(1:vv))
push!(agents, a)
push!(nodes_agents[a.current_node], a.id)
end
agents[1].infected = 1 # we start with one sick agent
agents[1].infection_time = 0
push!(nodes_infected[agents[1].current_node], 1)
Simulation(agents=agents, m=m, nodes_infected=nodes_infected,
nodes_agents=nodes_agents, infected_agents_count=[1])
end
s = Simulation(m)
Now we define the plotting function which is indeed similar to the previous one
function latlon(s::Simulation,map_g_point_id::Int64)
osm_node_ix = s.m.n[map_g_point_id]
lla = LLA(s.m.nodes[osm_node_ix], s.m.bounds)
return (lla.lat, lla.lon)
end
function plot_sim(s::Simulation; tiles="Stamen Toner" )
MAP_BOUNDS = [( s.m.bounds.min_y, s.m.bounds.min_x),( s.m.bounds.max_y, s.m.bounds.max_x)]
map_size = (abs(MAP_BOUNDS[1][1]-MAP_BOUNDS[2][1]), abs(MAP_BOUNDS[1][2]-MAP_BOUNDS[2][2]))
MAP_BOUNDSx9 = [ (MAP_BOUNDS[1][1]-map_size[1], MAP_BOUNDS[1][2]-map_size[1]),
(MAP_BOUNDS[2][1]+map_size[1], MAP_BOUNDS[2][2]+map_size[1])]
m_plot = flm.Map(tiles=tiles)
for e in edges(s.m.g)
flm.PolyLine( (latlon(s,e.src), latlon(s,e.dst)),
color="red", weight=5,
opacity=1).add_to(m_plot)
end
for v in 1:nv(s.m.g)
info = "<b>Node</B>\n<br> OSM id: $(s.m.n[v])\n <br>Node: $(v) "
flm.Circle(
latlon(s,v),
popup=info,
tooltip=info,
radius=10,
color="blue",
weight=3,
fill=true,
fill_color="blue"
).add_to(m_plot)
end
jitter = 2.5e-4
for agent in s.agents
info = "Agent: $(agent.id)\n<br>Infected: " *
(agent.infected ? "<b>YES</b>" : "NO")*
"\n<br>Current node: $(agent.current_node)"*
"\n<br>Previous node: $(agent.last_node)"*
"\n<br>Total distance travelled so far $(agent.total_route_len)"
loc = latlon(s,agent.current_node) .+ jitter.*randn(2)
sizex = agent.infected ? 0.0002 : 0.00016
sizey = agent.infected ? 0.0002 : 0.00022
flm.Rectangle(
[(loc[1]-sizex, loc[2]-sizey), (loc[1]+sizex, loc[2]+sizey)],
popup=info,
tooltip=info,
color=(agent.infected ? "green" : "black"),
weight=(agent.infected ? 5 : 1.5),
fill=(agent.infected ? false : true),
fill_opacity=(agent.infected ? 0.2 : 1.0),
fill_color=(agent.infected ? "green" : "#FAFAFA"),
).add_to(m_plot)
end
MAP_BOUNDS = [( s.m.bounds.min_y, s.m.bounds.min_x),( s.m.bounds.max_y, s.m.bounds.max_x)]
flm.Rectangle(MAP_BOUNDS, color="brown",weight=4).add_to(m_plot)
m_plot.fit_bounds(MAP_BOUNDS)
#m_plot.save("mysim.html") # uncomment to save to a file
m_plot
end
plot_sim(s)
Now we need to define how the agents move around
function step!(s::Simulation)
agent = rand(s.agents)
pop!(s.nodes_agents[agent.current_node], agent.id )
agent.total_route_len += 1
agent.infected && pop!(s.nodes_infected[agent.current_node], agent.id)
ns = LightGraphs.neighbors(s.m.g, agent.current_node)
new_node = rand(ns)
agent.last_node = agent.current_node
agent.current_node = new_node
new_infections = 0
if agent.infected
for id in s.nodes_agents[new_node]
if ! s.agents[id].infected
s.agents[id].infected = true
s.agents[id].infection_time = length(s.infected_agents_count)
push!(s.nodes_infected[new_node], id)
new_infections += 1
end
end
else
if length(s.nodes_infected[new_node]) > 0
agent.infected = true
agent.infection_time = length(s.infected_agents_count)
new_infections += 1
end
end
if agent.infected
push!(s.nodes_infected[new_node], agent.id)
end
push!(s.nodes_agents[new_node], agent.id)
push!(s.infected_agents_count,s.infected_agents_count[end]+new_infections)
end
We run the simulation for 10000 steps and watch the pandemic to develop
Random.seed!(1)
s = Simulation(m)
for i in 1:5000
step!(s)
end
plot_sim(s)
using Plots
Random.seed!(1)
function randres(N=200,steps=10000)
s = Simulation(m,N)
for i in 1:steps
step!(s)
end
s.infected_agents_count ./ N
end
res_total = hcat((randres(200,15000) for i in 1:200)...);
using Statistics
mvs = mean(res_total,dims=2)
plot(mvs, legend=:bottomright, lab="mean value")
plot!(mvs.+std(res_total,dims=2), linestyle=:dash, lab="upper")
plot!(mvs.-std(res_total,dims=2), linestyle=:dash, lab="lower")
res100 = hcat((randres(100,25000) for i in 1:200)...);
res200 = hcat((randres(200,25000) for i in 1:200)...);
res500 = hcat((randres(500,25000) for i in 1:200)...);
plot(collect(1:25001) ./100 , mean(res100,dims=2), legend=:bottomright, lab="100 agents")
plot!(collect(1:25001) ./200, mean(res200,dims=2), lab="200 agents")
plot!(collect(1:25001) ./500, mean(res500,dims=2), lab="500 agents")
title!("Stay at home!")