This is a comprehensive tutorial on the usage of LightOSM.jl, see the documentation for further details on the interface and methods.
First configure the logger:
using Logging
logger = SimpleLogger(stdout)
global_logger(logger);
Install the Julia dependencies:
using Pkg
Pkg.add("LightOSM")
Pkg.add("PyCall");
Set a seed so this tutorial is reproducible:
using Random
Random.seed!(1234);
For visualisation, this tutorial uses the Python package pydeck (a Python binding for Uber's graphics library deck.gl), and is called from Julia using PyCall.jl.
The official release of pydeck does not work with PyCall, you must install this forked version to the same virtual environment that PyCall is built with (might take a couple of minutes):
git clone --single-branch --branch pydeck/julia-pycall-binding git@github.com:captchanjack/deck.gl.git
cd deck.gl/bindings/pydeck
yarn bootstrap
pip install .
Download :drive network data as an object by searching a :place_name, notice there is an option to save that object to disk, the default data format is :osm (this is basically an .xml file but can be opened and visualized with JOSM, note JOSM can only open the file if keyword argument metadata=true when downloading):
using LightOSM
data = download_osm_network(:place_name,
place_name="melbourne, australia",
network_type=:drive,
save_to_file_location="melbourne_drive_network.osm");
With the downloaded data we can instantiate an OSMGraph object, which is a container for all the parsed data, DiGraph object, and KDTree object needed for shortest path and nearest node calculations. By default the graph object is a StaticDiGraph from StaticGraphs.jl as it is more memory-efficient:
g = graph_from_object(data, weight_type=:distance); # default weight_type is travel :time
# or g = graph_from_file("melbourne_drive_network.osm")
# or g = graph_from_download(:place_name, place_name="melbourne, australia")
Mapping of all Node objects:
g.nodes # node id => Node
Let's pick two Nodes for use throughout this tutorial:
n1 = g.nodes[443298633]
n2 = g.nodes[7620884608];
Node location:
n1.location
Haversine distance between n1 and n2:
distance(n1, n2) # km
Bearing from n1 to n2:
heading(n1, n2) # degrees from North
Node tags:
n1.tags
Mapping of all Way objects:
g.highways # way id => Highway
Highway nodes:
g.highways[13752475].nodes
Highway tags:
g.highways[13752475].tags
Mapping of all restriction objects:
g.restrictions # relation id => Restriction
Restriction tags:
g.restrictions[8221605].tags
First initialise the pydeck object, view state and Mapbox token, if you don't have one you can create a free account:
using PyCall
pydeck = pyimport("pydeck")
MAPBOX_TOKEN = "pk.eyJ1IjoiY2FwdGNoYW5qYWNrIiwiYSI6ImNrMzJ1enJoZjBueWwzY245ZDV0YjJ3Z3YifQ.VAWolOVu6eDYSnj3SC4NeQ"
MAXPBOX_STYLE = "mapbox://styles/captchanjack/ckepp735v2m2719lirl762qbi"
VIEWPORT_LOCATION = GeoLocation(-37.8142176, 144.9631608) # Melbourne (lat, lon)
view_state = pydeck.ViewState(longitude=VIEWPORT_LOCATION.lon,
latitude=VIEWPORT_LOCATION.lat,
zoom=13,
min_zoom=1,
max_zoom=25,
pitch=50,
bearing=-45);
Visualise Node objects using the ScatterplotLayer:
# Transform data
nodes_pydeck_data = [
Dict(
"ID" => id,
"Type" => "Node",
"Longitude" => node.location.lon,
"Latitude" => node.location.lat,
"Coordinates" => string([node.location.lon, node.location.lat]),
"Name" => "",
"Maxspeed" => "",
"Lanes" => "",
"Oneway" => "",
) for (id, node) in g.nodes
]
# Build the layer
nodes_layer = pydeck.Layer("ScatterplotLayer",
nodes_pydeck_data,
pickable=true,
opacity=0.8,
stroked=true,
filled=true,
line_width_min_pixels=5,
line_width_max_pixels=5,
line_width_scale=1,
auto_highlight=true,
get_position=["Longitude", "Latitude"],
get_radius=1,
radius_scale=1,
get_line_width=1,
get_line_color=[255, 98, 0, 255],
get_fill_color=[255,156,93, 255])
# Define tooltip
tooltip_style = Dict(
"color" => "white",
"border-radius" => "10px",
"border-color" => "dark grey",
"background-color" => "CadetBlue",
"font-family" => "Trebuchet MS",
"z-index" => 3,
"position" => "absolute"
)
nodes_tooltip = Dict(
"html" => "<b>ID:</b> {ID}<br><b>Coordinates:</b> {Coordinates}",
"style" => tooltip_style
)
# Build the deck
r = pydeck.Deck(layers=[nodes_layer],
initial_view_state=view_state,
mapbox_key=MAPBOX_TOKEN,
map_style=MAXPBOX_STYLE,
tooltip=nodes_tooltip)
# Save to .html file and display
r.to_html("nodes.html", notebook_display=true)
Add Way objects to the same deck using the PathLayer:
# Create a mapping of number of lanes to colours (monochromatic green)
COLOUR_MAPPING = Dict(
1 => [18, 231, 114],
2 => [53, 181, 53],
3 => [0, 107, 60],
4 => [0, 86, 63],
5 => [1, 50, 32]
)
DEFAULT_GREEN = [5, 102, 68] # more than 5 lanes
delete_quotes(str) = replace(str, r"'|\"" => " ") # deals with names with quotes (single or double apostrophe)
# Transform data
ways_pydeck_data = [
Dict(
"ID" => id,
"Longitude" => "",
"Latitude" => "",
"Coordinates" => "",
"Type" => "Highway - $(way.tags["highway"])",
"color" => get(COLOUR_MAPPING, way.tags["lanes"], DEFAULT_GREEN),
"Path" => [[g.nodes[n_id].location.lon, g.nodes[n_id].location.lat] for n_id in way.nodes],
"Name" => delete_quotes(get(way.tags, "name", "")),
"Maxspeed" => way.tags["maxspeed"],
"Lanes" => way.tags["lanes"],
"Oneway" => way.tags["oneway"]
) for (id, way) in g.highways
]
# Build the layer
ways_layer = pydeck.Layer("PathLayer",
ways_pydeck_data,
pickable=true,
get_color="color",
width_scale=1,
width_min_pixels=2,
get_path="Path",
get_width=2,
auto_highlight=true,
rounded=true)
# Define tooltip
ways_tooltip = Dict(
"html" => "
<b>ID:</b> {ID}<br>
<b>Type:</b> {Type}<br>
<b>Name:</b> {Name}<br>
<b>Coordinates:</b> {Coordinates}<br>
<b>Maxspeed:</b> {Maxspeed}<br>
<b>Lanes:</b> {Lanes}<br>
<b>Oneway:</b> {Oneway}<br>
",
"style" => tooltip_style
)
# Build the deck
r = pydeck.Deck(layers=[ways_layer, nodes_layer],
initial_view_state=view_state,
mapbox_key=MAPBOX_TOKEN,
map_style=MAXPBOX_STYLE,
tooltip=ways_tooltip)
# Save to .html file and display
r.to_html("nodes_and_ways.html", notebook_display=true)
Nearest Node calculations in this package is implemented with a K-D Tree data structure for fast querying. We can pick any point on the map (either a Node or a GeoLocation) to query N nearest Nodes and the straight line distances (Euclidean) to each of these neighbours.
If we use a Node's GeoLocation, the closest node will technically be itself:
nearest_node(g, n1.location) # returns a Tuple ([[nearest nodes]], [[distances in km]])
If a Node or Node.id is given as input, itself is not considered the closest node:
nearest_node(g, n1)
We can use any [latitude, longitude] pair to query:
nearest_node(g, [-37.8142176, 144.9631608])
We can query from multiple points:
nearest_node(g, [n1, n2]) # returns ([[n1_nbrs...], [n2_nbrs...]], [[n1_nbr_dists...], [n2_nbr_dists...]])
And search for multiple neighbours:
neighbours, distances = nearest_node(g, [n1, n2], 100)
Visualise neighbours and distances in an IconLayer:
# Create a new view state
VIEWPORT_LOCATION_NBR = GeoLocation(-37.8073,144.9771)
view_state_nbr = pydeck.ViewState(longitude=VIEWPORT_LOCATION_NBR.lon,
latitude=VIEWPORT_LOCATION_NBR.lat,
zoom=14,
min_zoom=1,
max_zoom=25,
pitch=50,
bearing=225);
# Define icon data
neighbour_icon_data = Dict(
# Icon taken from https://thenounproject.com/
"url" => "data:image/svg+xml;base64,PHN2ZyBoZWlnaHQ9JzMwMHB4JyB3aWR0aD0nMzAwcHgnICBmaWxsPSIjMUExQTFBIiB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHhtbG5zOnhsaW5rPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5L3hsaW5rIiB2ZXJzaW9uPSIxLjEiIHg9IjBweCIgeT0iMHB4IiB2aWV3Qm94PSIwIDAgMTAwIDEwMCIgc3R5bGU9ImVuYWJsZS1iYWNrZ3JvdW5kOm5ldyAwIDAgMTAwIDEwMDsiIHhtbDpzcGFjZT0icHJlc2VydmUiPjxwYXRoIGQ9Ik04NS45LDM4LjVDODUuOSwxOC43LDY5LjgsMi42LDUwLDIuNlMxNC4xLDE4LjcsMTQuMSwzOC41YzAsMTYuMywxMC44LDMwLDI1LjYsMzQuNGw3LjUsMTggIGMxLjIsMi4yLDQuMywyLjIsNS41LDBsNy41LTE4Qzc1LjEsNjguNSw4NS45LDU0LjgsODUuOSwzOC41eiI+PC9wYXRoPjwvc3ZnPg==",
"width" => 242,
"height" => 242,
"anchorY" => 242,
"mask" => true # allows color to be altered with get_color kwarg
)
centroid_icon_data = Dict(
"url" => "https://cdn.iconscout.com/icon/premium/png-512-thumb/destination-flag-8-902948.png",
"width" => 242,
"height" => 242,
"anchorY" => 242,
)
# Define red-green colour scale, where 0 <= n <= 100
red(n) = Int(round(255 * n / 100))
green(n) = Int(round(255 * (100 - n) / 100))
blue(n) = 0
# Transform data
nbrs = [(neighbours...)...] # flatten
dists = [(distances...)...] # flatten
max_dist = max(dists...)
scaled_dist(d) = (d / max_dist) * 100
neighbours_pydeck_data = [
Dict(
"ID" => n,
"Type" => "Neighbour",
"Longitude" => g.nodes[n].location.lon,
"Latitude" => g.nodes[n].location.lat,
"Coordinates" => string([g.nodes[n].location.lon, g.nodes[n].location.lat]),
"Distance" => string(round(dists[i]*1000, digits=2))*"m",
"Icon" => neighbour_icon_data,
"Colour" => [red(scaled_dist(dists[i])), green(scaled_dist(dists[i])), blue(scaled_dist(dists[i]))]
) for (i, n) in enumerate(nbrs)
]
centroid_pydeck_data = [
Dict(
"ID" => n.id,
"Type" => "Centroid",
"Longitude" => n.location.lon,
"Latitude" => n.location.lat,
"Coordinates" => string([n.location.lon, n.location.lat]),
"Icon" => centroid_icon_data,
"Distance" => "0m"
) for n in [n1, n2]
]
# Build the layer
neighbours_layer = pydeck.Layer("IconLayer",
data=neighbours_pydeck_data,
get_icon="Icon",
get_size=5,
size_min_pixels=50,
size_max_pixels=30,
filled=true,
get_color="Colour",
opacity=0.5,
get_position=["Longitude", "Latitude"],
pickable=true,
auto_highlight=true)
centroid_layer = pydeck.Layer("IconLayer",
data=centroid_pydeck_data,
get_icon="Icon",
get_size=5,
size_min_pixels=60,
size_max_pixels=60,
filled=true,
get_color=[72, 209, 204],
get_position=["Longitude", "Latitude"],
pickable=true,
auto_highlight=true)
# Define tooltip
neighbours_tooltip = Dict(
"html" => "
<b>ID:</b> {ID}<br>
<b>Type:</b> {Type}<br>
<b>Coordinates:</b> {Coordinates}<br>
<b>Distance From Centroid:</b> {Distance}<br>
",
"style" => tooltip_style
)
# Build the deck
r = pydeck.Deck(layers=[neighbours_layer, centroid_layer],
initial_view_state=view_state_nbr,
mapbox_key=MAPBOX_TOKEN,
map_style=MAXPBOX_STYLE,
tooltip=neighbours_tooltip)
# Save to .html file and display
r.to_html("neighbours.html", notebook_display=true)
To calculate the shortest path between two Nodes, we can use the Dijkstra or A* algorithm. These differ to those implemented in LightGraphs and OpenStreetMapX.jl as LightOSM.jl takes into account turn restrictions.
Let's start by calculating the shortest path between n1 and n2:
path = shortest_path(g, n1, n2, algorithm=:dijkstra)
Retrieve the Edge weights for the path (:distance was selected as the :weight_type earlier):
weights = weights_from_path(g, path) # km
Cumulative edge distance:
cum_weights = cumsum(weights)
Total path distance in km is therefore the last index:
total_distance = cum_weights[end]
Edges are segments of a Way and are defined as adjacent origin-destination Node pairs:
edges = [[path[i], path[i + 1]] for i in 1:length(path) - 1]
Use Edges to find the corresponding Ways:
way_ids = [g.edge_to_highway[e] for e in edges]
Now with our Nodes, Edges and Ways we can plot the shortest path using a LineLayer:
# Define icon data
start_finish_icons = [
Dict(
"url" => "data:image/svg+xml;base64,PHN2ZyBoZWlnaHQ9JzMwMHB4JyB3aWR0aD0nMzAwcHgnICBmaWxsPSIjMUExQTFBIiB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIGRhdGEtbmFtZT0iTGF5ZXIgMSIgdmlld0JveD0iMCAwIDEwMCAxMDAiIHg9IjBweCIgeT0iMHB4Ij48dGl0bGU+VmVoaWNsZXM8L3RpdGxlPjxwYXRoIGQ9Ik0xMDAsNTcuNTZ2LS4wOGE0Ljc1LDQuNzUsMCwwLDAtLjE5LTFDOTguNTEsNTAuNDgsOTQsNDUuMDYsODgsNDMuNTJBNjQuMzksNjQuMzksMCwwLDAsNzksNDJhNy4yOSw3LjI5LDAsMCwxLTUuNzEtMy4yNmMtNS4zMy03LjYyLTEyLjc3LTEyLTIxLjg5LTEzLjY1YTc5LjYxLDc5LjYxLDAsMCwwLTI1LjgyLjA2LDE2LjU3LDE2LjU3LDAsMCwwLTExLjkzLDcuNzdjLTEuMDksMS42OS0xLjk0LDMuNTYtMy4xMiw1LjE4YTE5LjY2LDE5LjY2LDAsMCwxLTQuMSw0LjY4QTEyLjYyLDEyLjYyLDAsMCwwLC45MSw1MC44OEMuNSw1My4wOS4yMyw1NS4zMiwwLDU3LjU2SC4wNWE1LjE5LDUuMTksMCwwLDAsMCwuNjh2NWEyLDIsMCwwLDAsMiwxLjk1aDcuNWExMi4wOSwxMi4wOSwwLDAsMSwyNC4wNSwwSDY5LjM3YzAtLjE1LDAtLjMsMC0uNDVhMTIuMSwxMi4xLDAsMCwxLDI0LjE5LDBjMCwuMTUsMCwuMywwLC40NWg0LjUzYTIsMiwwLDAsMCwyLTEuOTVWNTcuODJhMi4zNSwyLjM1LDAsMCwwLDAtLjI2Wk0zOS40LDQ0LjExYy02LDAtMTEuODQuMTUtMTcuNjctLjA3LTIuODctLjEtMy41NC0xLjYyLTIuNDgtNC4zNCwzLjUyLTksOS4yLTExLjk0LDIwLjE1LTEwLjI2Wm0yNS44MSwwYy02LjUxLjEtMTMsMC0xOS44OSwwVjI5LjU5YzcuNTEtMS4zOSwxOC43OCw0LjQsMjIsMTEuMjZDNjguMyw0Myw2Ny43Myw0NC4xMiw2NS4yMSw0NC4xNloiPjwvcGF0aD48cGF0aCBkPSJNMjEuNDcsNTdBOS41MSw5LjUxLDAsMSwwLDMxLDY2LjQ3LDkuNSw5LjUsMCwwLDAsMjEuNDcsNTdabTAsMTMuNzNhNC4yMyw0LjIzLDAsMSwxLDQuMjMtNC4yMkE0LjIyLDQuMjIsMCwwLDEsMjEuNDcsNzAuNjlaIj48L3BhdGg+PHBhdGggZD0iTTkxLDY0Ljc2YTkuNTEsOS41MSwwLDEsMC05LjUsOS41MUE5LjUsOS41LDAsMCwwLDkxLDY0Ljc2Wm0tMTMuNzMsMEE0LjIzLDQuMjMsMCwxLDEsODEuNDUsNjksNC4yNCw0LjI0LDAsMCwxLDc3LjIyLDY0Ljc2WiI+PC9wYXRoPjwvc3ZnPg==",
"width" => 242,
"height" => 242,
"anchorY" => 242,
"mask" => true
),
Dict(
"url" => "data:image/svg+xml;base64,PHN2ZyBoZWlnaHQ9JzMwMHB4JyB3aWR0aD0nMzAwcHgnICBmaWxsPSIjMUExQTFBIiB4bWxuczp4PSJodHRwOi8vbnMuYWRvYmUuY29tL0V4dGVuc2liaWxpdHkvMS4wLyIgeG1sbnM6aT0iaHR0cDovL25zLmFkb2JlLmNvbS9BZG9iZUlsbHVzdHJhdG9yLzEwLjAvIiB4bWxuczpncmFwaD0iaHR0cDovL25zLmFkb2JlLmNvbS9HcmFwaHMvMS4wLyIgeG1sbnM9Imh0dHA6Ly93d3cudzMub3JnLzIwMDAvc3ZnIiB4bWxuczp4bGluaz0iaHR0cDovL3d3dy53My5vcmcvMTk5OS94bGluayIgdmVyc2lvbj0iMS4xIiB4PSIwcHgiIHk9IjBweCIgdmlld0JveD0iMCAwIDY0IDY0IiBlbmFibGUtYmFja2dyb3VuZD0ibmV3IDAgMCA2NCA2NCIgeG1sOnNwYWNlPSJwcmVzZXJ2ZSI+PG1ldGFkYXRhPjxzZncgeG1sbnM9Imh0dHA6Ly9ucy5hZG9iZS5jb20vU2F2ZUZvcldlYi8xLjAvIj48c2xpY2VzPjwvc2xpY2VzPjxzbGljZVNvdXJjZUJvdW5kcyB4PSItNzg2NSIgeT0iLTIzMCIgd2lkdGg9IjE2MzgzIiBoZWlnaHQ9IjIyOCIgYm90dG9tTGVmdE9yaWdpbj0idHJ1ZSI+PC9zbGljZVNvdXJjZUJvdW5kcz48L3Nmdz48L21ldGFkYXRhPjxwYXRoIGQ9Ik01Myw1SDM5SDI1Yy0wLjU1MiwwLTEsMC40NDgtMSwxdjEwdjEwYzAsMC41NTIsMC40NDgsMSwxLDFoMTRoMTRjMC41NTIsMCwxLTAuNDQ4LDEtMVYxNlY2QzU0LDUuNDQ4LDUzLjU1Miw1LDUzLDV6ICAgTTUyLDE1SDQwVjdoMTJWMTV6IE0yNiwxN2gxMnY4SDI2VjE3eiI+PC9wYXRoPjxwYXRoIGQ9Ik0yOC45OTcsNjAuNzg3Yy0wLjEwNC0zLjk3Ni0zLjEyNi03LjIzMS02Ljk5Ny03LjcxN1YzYzAtMC41NTItMC40NDgtMS0xLTFzLTEsMC40NDgtMSwxdjUwLjA2OSAgYy0zLjk0LDAuNDk1LTcsMy44NTktNyw3LjkzMWMwLDAuNTUyLDAuNDQ4LDEsMSwxaDE0YzAuMDA3LDAsMC4wMTMtMC4wMDEsMC4wMiwwYzAuNTUyLDAsMS0wLjQ0OCwxLTEgIEMyOS4wMiw2MC45MjcsMjkuMDEyLDYwLjg1NSwyOC45OTcsNjAuNzg3eiI+PC9wYXRoPjwvc3ZnPg==",
"width" => 242,
"height" => 242,
"anchorY" => 242,
"mask" => true,
)
]
# Define red-blue colour scale, where 0 <= n <= 100
red2(n) = Int(round(255 * n / 100))
green2(n) = 0
blue2(n) = Int(round(255 * (100 - n) / 100))
scaled_w(d) = (d / total_distance) * 100
# Transform data
labels = ["Origin Node", "Destination Node"]
start_finish_data = [
Dict(
"ID" => n.id,
"Type" => labels[i],
"Longitude" => n.location.lon,
"Latitude" => n.location.lat,
"Coordinates" => string([n.location.lon, n.location.lat]),
"Icon" => start_finish_icons[i],
"Distance" => "",
"Name" => "",
"Maxspeed" => "",
"Lanes" => "",
"Oneway" => "",
"Edge" => ""
) for (i, n) in enumerate([n1, n2])
]
shortest_path_pydeck_data = [
Dict(
"ID" => way_ids[i],
"Edge" => string([origin, destination]),
"Type" => "Highway - $(g.highways[way_ids[i]].tags["highway"])",
"start" => [g.nodes[origin].location.lon, g.nodes[origin].location.lat],
"end" => [g.nodes[destination].location.lon, g.nodes[destination].location.lat],
"Colour" => [red2(scaled_w(cum_weights[i])), green2(scaled_w(cum_weights[i])), blue2(scaled_w(cum_weights[i]))],
"Name" => delete_quotes(get(g.highways[way_ids[i]].tags, "name", "")),
"Maxspeed" => g.highways[way_ids[i]].tags["maxspeed"],
"Lanes" => g.highways[way_ids[i]].tags["lanes"],
"Oneway" => g.highways[way_ids[i]].tags["oneway"],
"Coordinates" => "",
"Distance" => string(round(cum_weights[i] * 1000, digits=2)) * "m"
) for (i, (origin, destination)) in enumerate(edges)
]
# Build the layer
start_finish_layer = pydeck.Layer("IconLayer",
data=start_finish_data,
get_icon="Icon",
get_size=5,
size_min_pixels=60,
size_max_pixels=60,
filled=true,
get_position=["Longitude", "Latitude"],
pickable=true,
highlight_color=[106, 110, 117],
auto_highlight=true)
shortest_path_layer = pydeck.Layer("LineLayer",
shortest_path_pydeck_data,
get_source_position="start",
get_target_position="end",
get_color="Colour",
get_width=7,
picking_radius=10,
auto_highlight=true,
highlight_color=[106, 110, 117],
rounded=true,
pickable=true)
# Define tooltip
shortest_path_tooltip = Dict(
"html" => "
<b>ID:</b> {ID}<br>
<b>Edge:</b> {Edge}<br>
<b>Type:</b> {Type}<br>
<b>Coordinates:</b> {Coordinates}<br>
<b>Name:</b> {Name}<br>
<b>Maxspeed:</b> {Maxspeed}<br>
<b>Lanes:</b> {Lanes}<br>
<b>Oneway:</b> {Oneway}<br>
<b>Distance From Origin:</b> {Distance}<br>
",
"style" => tooltip_style
)
# Build the deck
r = pydeck.Deck(layers=[start_finish_layer, shortest_path_layer],
initial_view_state=view_state_nbr,
mapbox_key=MAPBOX_TOKEN,
map_style=MAXPBOX_STYLE,
tooltip=shortest_path_tooltip)
# Save to .html file and display
r.to_html("shortest_path.html", notebook_display=true)
Now let's generate some random origin-destination Node pairs:
n_paths = 1000
rand_o_d_indices = rand(1:length(g.nodes), n_paths, 2)
rand_o_d_nodes = [[g.index_to_node[o], g.index_to_node[d]] for (o, d) in eachrow(rand_o_d_indices) if o != d]
Calculate shortest path between each origin-destination Node pair:
paths = []
@time for (o, d) in rand_o_d_nodes
try
p = shortest_path(g, o, d, algorithm=:dijkstra)
push!(paths, p)
catch
# Error exception will be thrown if path does not exist from origin to destination node
end
end
Retrieve total distance in km for each path:
total_distances = round.(sum.([weights_from_path(g, p) for p in paths]), digits=2)
We can plot these shortest paths using the TripsLayer, unfortunately unlike its parent library deck.gl, pydeck currently does not support dynamic animation of the TripsLayer against current_time (this requires some JavaScript methods such as requestAnimationFrame):
# Transform data
trips_pydeck_data = [
Dict(
"Coordinates" => [[g.nodes[node_id].location.lon, g.nodes[node_id].location.lat] for node_id in path],
"Timestamps" => collect(0:length(path)-1),
"Colour" => [rand(1:255), rand(1:255), rand(1:255)],
"Distance" => total_distances[i],
"Origin" => path[1],
"Destination" => path[end]
) for (i, path) in enumerate(paths)
]
# Build the layer
trips_layer = pydeck.Layer("TripsLayer",
trips_pydeck_data,
get_path="Coordinates",
get_timestamps="Timestamps",
get_color="Colour",
opacity=0.8,
width_min_pixels=5,
rounded=true,
trail_length=50,
current_time=100,
pickable=true,
auto_highlight=true)
# Define tooltip
trips_tooltip=Dict(
"html" => "
<b>Path Distance:</b> {Distance}km<br>
<b>Origin:</b> {Origin}<br>
<b>Destination:</b> {Destination}<br>
",
"style" => tooltip_style
)
# Build the deck
r = pydeck.Deck(layers=[trips_layer],
initial_view_state=view_state,
mapbox_key=MAPBOX_TOKEN,
map_style=MAXPBOX_STYLE,
tooltip=trips_tooltip)
# Save to .html file and display
r.to_html("trips.html", notebook_display=true)
Similar to downloading an OpenStreetMap network, we can also download Building polygons by searching with a :place_name, centroid :point or :bbox:
buildings_data = download_osm_buildings(:place_name,
place_name="melbourne, australia",
save_to_file_location="melbourne_buildings.osm");
We can then parse and instantiate Building objects from the data downloaded (either using the in-memory data object, a saved file or a direct download method):
buildings = buildings_from_object(buildings_data)
# or buildings = buildings_from_file("melbourne_buildings.osm")
# or buildings = buildings_from_download(:place_name, place_name="melbourne, australia")
Buildings consist of an array of Polygons, the first element is always the outer ring, followed by any optional inner rings (inner rings are holes in the outer ring), each with their own set of Nodes:
buildings[127595640].polygons[1].nodes
buildings[127595640].polygons[1].is_outer
Building metadata tags:
buildings[127595640].tags # height is in metres
Visualise Buildings with a PolygonLayer:
# Transform data
max_height = max([b.tags["height"] for (id, b) in buildings]...)
rbg(d) = 255 - Int(round((d / max_height * 255)))
buildings_pydeck_data = [
Dict(
"Polygons" => [[[node.location.lon, node.location.lat] for node in poly.nodes] for poly in b.polygons], # array of polygons (i.e. array of array of coordinates)
"Height" => b.tags["height"],
"Colour" => [rand(1:255), rand(1:255), rand(1:255)],
"Name" => delete_quotes(string(get(b.tags, "name", ""))),
"Description" => delete_quotes(string(get(b.tags, "description", ""))),
"Website" => delete_quotes(string(get(b.tags, "website", ""))),
"Colour" => [rbg(b.tags["height"]), rbg(b.tags["height"]), rbg(b.tags["height"])], # grey scale
"ID" => id
) for (id, b) in buildings
]
# Build the layer
buildings_layer = pydeck.Layer("PolygonLayer",
buildings_pydeck_data,
id="geojson",
opacity=0.2,
stroked=false,
get_polygon="Polygons",
filled=true,
extruded=true,
wireframe=true,
get_elevation="Height",
get_fill_color="Colour",
get_line_color=[255, 255, 255],
auto_highlight=true,
pickable=true,
)
# Define tooltip
buildings_tooltip = Dict(
"html" => "
<b>ID:</b> {ID}<br>
<b>Name:</b> {Name}<br>
<b>Description:</b> {Description}<br>
<b>Height:</b> {Height}m<br>
<b>Website:</b> {Website}<br>
",
"style" => tooltip_style
)
# Build the deck
r = pydeck.Deck(layers=[buildings_layer],
initial_view_state=view_state,
mapbox_key=MAPBOX_TOKEN,
map_style=MAXPBOX_STYLE,
tooltip=buildings_tooltip)
# Save to .html file and display
r.to_html("buildings.html", notebook_display=true)