[WIP] Simplify Graph

Project.toml

@@ -4,14 +4,17 @@ authors = ["Jack Chan <[email protected]>"]
 version = "0.1.20"
 
 [deps]
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
 MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticGraphs = "4c8beaf5-199b-59a0-a7f2-21d17de635b6"

example.jl

@@ -0,0 +1,36 @@
+using LightOSM, GeoInterface, Plots, DataFrames
+using ArchGDAL: createmultilinestring
+
+g = graph_from_download(
+    :place_name, place_name="moabit, berlin germany",
+    network_type=:bike
+)
+# reverse coordinates for plotting
+reverse!.(g.node_coordinates)
+
+g_simple, weights, node_gdf, edge_gdf = simplify_graph(g)
+
+# join edges in mulitlinestring for faster plotting
+all_edges = createmultilinestring(coordinates.(edge_gdf.geom))
+
+# node validation
+
+# nodes from original graph
+plot(all_edges, color=:black, size=(1200,800))
+scatter!(first.(g.node_coordinates), last.(g.node_coordinates), color=:red)
+
+# nodes from simplified graph
+plot(all_edges, color=:black, size=(1200,800))
+scatter!(node_gdf.geom, color=:green)
+
+
+# edge validation
+
+function highway_gdf(osmg::OSMGraph)
+    function _geometrize_way(way)
+        createlinestring(map(id -> coordinates(osmg.nodes[id]), way.nodes))
+    end
+    geom = map(way -> _geometrize_way(way), values(osmg.highways))
+    return DataFrame(; id = collect(keys(osmg.highways)), geom)
+end
+

src/LightOSM.jl

@@ -4,17 +4,22 @@ using Parameters
 using DataStructures: DefaultDict, OrderedDict, MutableLinkedList, PriorityQueue, dequeue!, dequeue_pair!
 using Statistics: mean
 using SparseArrays: SparseMatrixCSC, sparse
-using Graphs: AbstractGraph, DiGraph, nv, outneighbors, weakly_connected_components, vertices
+using Graphs: AbstractGraph, DiGraph, nv, outneighbors, weakly_connected_components, vertices, all_neighbors, indegree, outdegree, add_edge!
 using StaticGraphs: StaticDiGraph
 using SimpleWeightedGraphs: SimpleWeightedDiGraph
 using MetaGraphs: MetaDiGraph
 using NearestNeighbors: KDTree, knn
 using HTTP
 using JSON
 using LightXML
+using DataFrames
+using GeoInterface
+using RecipesBase
 
 export GeoLocation,
+       AbstractOSMGraph,
        OSMGraph,
+       SimplifiedOSMGraph,
        Node,
        Way,
        Restriction,
@@ -33,7 +38,11 @@ export GeoLocation,
        download_osm_buildings,
        buildings_from_object,
        buildings_from_download,
-       buildings_from_file
+       buildings_from_file,
+       simplify_graph,
+       node_gdf,
+       edge_gdf,
+       highway_gdf
 
 include("types.jl")
 include("constants.jl")
@@ -46,5 +55,8 @@ include("traversal.jl")
 include("shortest_path.jl")
 include("nearest_node.jl")
 include("buildings.jl")
+include("simplification.jl")
+include("geodataframes.jl")
+include("plotrecipes.jl")
 
 end # module

src/geodataframes.jl

@@ -0,0 +1,38 @@
+coordinates(node::Node) = [node.location.lon, node.location.lat]
+
+function node_gdf(g::OSMGraph)
+    ids = collect(keys(g.nodes))
+    geom = map(ids) do id
+       coordinates(g.nodes[id])
+    end
+    return DataFrame(;id=ids, geom=Point.(geom))
+end
+
+function highway_gdf(g::OSMGraph)
+    ids = collect(keys(g.highways))
+    _way_coordinates(way) = map(way.nodes) do id
+        coordinates(g.nodes[id])
+    end
+    geom = map(id -> _way_coordinates(g.highways[id]), ids)
+    return DataFrame(;id=ids, geom=LineString.(geom))
+end
+
+function node_gdf(sg::SimplifiedOSMGraph)
+    ids = collect(keys(sg.node_to_index))
+    geom = map(ids) do id
+        coordinates(sg.parent.nodes[id])
+    end
+    return DataFrame(;id=ids, geom=Point.(geom))
+end
+
+highway_gdf(sg::SimplifiedOSMGraph) = highway_gdf(sg.parent)
+
+function edge_gdf(sg::SimplifiedOSMGraph)
+    edge_ids = collect(keys(sg.edges))
+    geom = map(edge_ids) do edge
+        path = sg.edges[edge]
+        reverse.(sg.parent.node_coordinates[path])
+    end
+    u, v, key = map(i -> getindex.(edge_ids, i), 1:3)
+    return DataFrame(;u, v, key, geom=LineString.(geom))
+end

src/plotrecipes.jl

@@ -0,0 +1,13 @@
+function aspect_ratio(g::OSMGraph) 
+    max_y, min_y = extrema(first, g.node_coordinates)
+    mid_y = (max_y + min_y)/2
+    return 1/cos(mid_y * pi/180)
+end
+aspect_ratio(sg::SimplifiedOSMGraph) = aspect_ratio(sg.parent)
+
+
+RecipesBase.@recipe function f(g::AbstractOSMGraph)
+    color --> :black
+    aspect_ratio --> aspect_ratio(g)
+    MultiLineString(GeoInterface.coordinates.(highway_gdf(g).geom))
+end

src/simplification.jl

@@ -0,0 +1,122 @@
+
+#adapted from osmnx: https://github.com/gboeing/osmnx/blob/main/osmnx/simplification.py
+
+"""
+Predicate wether v is an edge endpoint in the simplified version of g
+"""
+function is_endpoint(g::AbstractGraph, v)
+    neighbors = all_neighbors(g, v)
+    if v in neighbors # has self loop
+        return true
+    elseif outdegree(g, v) == 0 || indegree(g, v) == 0 # sink or source
+        return true
+    elseif length(neighbors) != 2 || indegree(g, v) != outdegree(g, v) # change to/from one way
+        return true
+    end
+    return false
+end
+
+"""
+iterator over all endpoints in g
+"""
+endpoints(g::AbstractGraph) = (v for v in vertices(g) if is_endpoint(g, v))
+
+"""
+iterator over all paths in g which can be contracted
+"""
+function paths_to_reduce(g::AbstractGraph)
+    (path_to_endpoint(g, (u, v)) for u in endpoints(g) for v in outneighbors(g, u))
+end
+
+"""
+path to the next endpoint starting in edge (ep, ep_succ)
+"""
+function path_to_endpoint(g::AbstractGraph, (ep, ep_succ)::Tuple{T,T}) where {T<:Integer}
+    path = [ep, ep_succ]
+    head = ep_succ
+    # ep_succ not in endpoints -> has 2 neighbors and degree 2 or 4
+    while !is_endpoint(g, head)
+        neighbors = [n for n in outneighbors(g, head) if n != path[end-1]]
+        @assert length(neighbors) == 1 "found unmarked endpoint!"
+        head, = neighbors
+        push!(path, head)
+        (head == ep) && return path # self loop
+    end
+    return path
+end
+
+"""
+Return the total weight of a path given as a Vector of Ids.
+"""
+function total_weight(g::OSMGraph, path::Vector{<:Integer})
+    sum((g.weights[path[[i, i+1]]...] for i in 1:length(path)-1))
+end
+
+function ways_in_path(g::OSMGraph, path::Vector{<:Integer})
+    ways = Set{Int}()
+    for i in 1:(length(path)-1)
+        edge = [g.index_to_node[path[i]], g.index_to_node[path[i+1]]]
+        push!(ways, g.edge_to_way[edge])
+    end
+    return collect(ways)
+end
+
+"""
+Build a new graph which simplifies the topology of osmg.graph.
+The resulting graph only contains intersections and dead ends from the original graph.
+The geometry of the contracted nodes is kept in the edge_gdf DataFrame
+"""
+function simplify_graph(osmg::OSMGraph{U, T, W}) where {U, T, W}
+    g = osmg.graph
+    relevant_nodes = collect(endpoints(g))
+    n_relevant = length(relevant_nodes)
+    graph = DiGraph(n_relevant) 
+    weights = similar(osmg.weights, (n_relevant, n_relevant))
+    node_coordinates = Vector{Vector{W}}(undef, n_relevant)
+    node_to_index = OrderedDict{T,U}()
+    index_to_node = OrderedDict{U,T}()
+
+    index_mapping = Dict{U,U}()
+    for (new_i, old_i) in enumerate(relevant_nodes)
+        index_mapping[old_i] = new_i
+        node_coordinates[new_i] = osmg.node_coordinates[old_i]
+        node = osmg.index_to_node[old_i]
+        index_to_node[new_i] = node
+        node_to_index[node] = new_i
+    end
+
+    edges = Dict{NTuple{3,U}, Vector{U}}()
+    edge_count = Dict{Tuple{U,U}, Int}()
+    for path in paths_to_reduce(g)
+        u = index_mapping[first(path)]
+        v = index_mapping[last(path)]
+        path_weight = total_weight(osmg, path)
+        if add_edge!(graph, (u, v))
+            key = 0
+            weights[u, v] = path_weight
+            edge_count[u,v] = 1
+        else # parallel edge
+            key = edge_count[u,v]
+            edge_count[u,v] += 1
+            weights[u, v] = min(path_weight, weights[u, v])
+        end
+        edges[u,v,key] = path
+    end
+
+    edge_to_way = Dict{NTuple{3,U}, Vector{T}}()
+    for (edge, path) in edges
+        edge_to_way[edge] = ways_in_path(osmg, path)
+    end
+
+    return SimplifiedOSMGraph(
+                osmg,
+                node_coordinates,
+                node_to_index,
+                index_to_node,
+                edge_to_way,
+                graph,
+                edges,
+                weights,
+                nothing
+            )
+end

src/types.jl

@@ -99,7 +99,9 @@ Container for storing OpenStreetMap node, way, relation and graph related obejct
 - `kdtree::Union{KDTree,Nothing}`: KDTree used to calculate nearest nodes.
 - `weight_type::Union{Symbol,Nothing}`: Either `:distance`, `:time` or `:lane_efficiency`.
 """
-@with_kw mutable struct OSMGraph{U <: Integer,T <: Integer,W <: Real}
+
+abstract type AbstractOSMGraph end
+@with_kw mutable struct OSMGraph{U <: Integer,T <: Integer,W <: Real} <: AbstractOSMGraph
     nodes::Dict{T,Node{T}} = Dict{T,Node{T}}()
     node_coordinates::Vector{Vector{W}} = Vector{Vector{W}}() # needed for astar heuristic
     ways::Dict{T,Way{T}} = Dict{T,Way{T}}()
@@ -116,6 +118,18 @@ Container for storing OpenStreetMap node, way, relation and graph related obejct
     weight_type::Union{Symbol,Nothing} = nothing
 end
 
+struct SimplifiedOSMGraph{U <: Integer,T <: Integer,W <: Real} <: AbstractOSMGraph
+    parent::OSMGraph{U,T,W}
+    node_coordinates::Vector{Vector{W}} # needed for astar heuristic
+    node_to_index::OrderedDict{T,U}
+    index_to_node::OrderedDict{U,T}
+    edge_to_way::Dict{NTuple{3,U},Vector{T}}
+    graph::Union{AbstractGraph,Nothing}
+    edges::Dict{NTuple{3, U}, Vector{U}}
+    weights::Union{SparseMatrixCSC{W,U},Nothing}
+    dijkstra_states::Union{Vector{Vector{U}},Nothing}
+end
+
 function Base.getproperty(g::OSMGraph, field::Symbol)
     # Ensure renaming of "highways" to "ways" is backwards compatible
     if field === :highways