Add nearest way and R-tree

Author: jarodlam

Project.toml

@@ -1,7 +1,7 @@
 name = "LightOSM"
 uuid = "d1922b25-af4e-4ba3-84af-fe9bea896051"
 authors = ["Jack Chan <jchan2@deloitte.com.au>"]
-version = "0.2.5"
+version = "0.2.6"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -15,6 +15,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 QuickHeaps = "30b38841-0f52-47f8-a5f8-18d5d4064379"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpatialIndexing = "d4ead438-fe20-5cc5-a293-4fd39a41b74c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StaticGraphs = "4c8beaf5-199b-59a0-a7f2-21d17de635b6"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -30,6 +31,7 @@ NearestNeighbors = "0.4.6"
 Parameters = "0.12.1"
 QuickHeaps = "0.1.1"
 SimpleWeightedGraphs = "1.2.0"
+SpatialIndexing = "0.1.3"
 StaticArrays = "1.4.6"
 StaticGraphs = "0.3.0"
 julia = "1"

docs/src/graph_utilities.md

@@ -1,4 +1,4 @@
-# [`Graph Utilities`](@ref) Methods
+# Graph Utilities
 
 ```@docs
 index_to_node_id

docs/src/index.md

@@ -11,10 +11,11 @@ Pages = [
     "create_graph.md",
     "shortest_path.md",
     "nearest_node.md",
+    "nearest_way.md",
     "download_buildings.md",
     "create_buildings.md",
-    "geolocation.md"
-    "graph_utilities.md"
+    "geolocation.md",
+    "graph_utilities.md",
     "defaults.md"
 ]
 ```

docs/src/nearest_node.md

@@ -2,4 +2,5 @@
 
 ```@docs
 nearest_node
+nearest_nodes
 ```

docs/src/nearest_way.md

@@ -0,0 +1,7 @@
+# Nearest Node
+
+```@docs
+nearest_way
+nearest_ways
+nearest_point_on_way
+```

src/LightOSM.jl

@@ -14,11 +14,13 @@ using HTTP
 using JSON
 using LightXML
 using StaticArrays
+using SpatialIndexing
 
 export GeoLocation,
        OSMGraph,
        Node,
        Way,
+       EdgePoint,
        Restriction,
        Building,
        PathAlgorithm,
@@ -46,6 +48,9 @@ export GeoLocation,
        path_from_parents,
        nearest_node,
        nearest_nodes,
+       nearest_way,
+       nearest_ways,
+       nearest_point_on_way,
        download_osm_buildings,
        buildings_from_object,
        buildings_from_download,
@@ -75,6 +80,7 @@ include("graph_utilities.jl")
 include("traversal.jl")
 include("shortest_path.jl")
 include("nearest_node.jl")
+include("nearest_way.jl")
 include("buildings.jl")
 include("subgraph.jl")
 

src/constants.jl

@@ -228,7 +228,7 @@ optional.
     values are used instead based on the value of the `highway` way tag. If no 
     `highway` way tag is available, the value for `"other"` is used. Unit is km/h. 
     Default value: 
-    ```julia
+    ```
     Dict(
         "motorway" => 100,
         "trunk" => 100,
@@ -244,7 +244,7 @@ optional.
     values are used instead based on the value of the `highway` way tag. If no 
     `highway` way tag is available, the value for `"other"` is used. 
     Default value: 
-    ```julia
+    ```
     Dict(
         "motorway" => 3,
         "trunk" => 3,
@@ -259,7 +259,7 @@ optional.
 - `lane_efficiency::AbstractDict{<:Integer,<:Real}`: Gives the lane efficiency based on 
     number of lanes. `1.0` is used for any number of lanes not specified here.
     Default value:
-    ```julia
+    ```
     LANE_EFFICIENCY = Dict(
         1 => 0.7,
         2 => 0.8,
@@ -270,12 +270,12 @@ optional.
 - `building_height_per_level::Integer`: If the `height` building tag is not available, 
     it is calculated by multiplying this value by the number of levels from the 
     `building:levels` tag. Unit is metres. Default value:
-    ```julia
+    ```
     4
     ```
 - `max_building_levels::Integer`: If the `building:levels` tag is not available, a number 
     is randomly chosen between 1 and this value. Default value:
-    ```julia
+    ```
     3
     ```
 """

src/geometry.jl

@@ -213,3 +213,51 @@ function bounding_box_from_point(point::GeoLocation, radius::Number)::NamedTuple
     bottom_left, top_right = calculate_location([point, point], [225, 45], [radius, radius])
     return (minlat = bottom_left.lat, minlon = bottom_left.lon, maxlat = top_right.lat, maxlon = top_right.lon)
 end
+
+"""
+    nearest_point_on_line(x1::T, 
+                          y1::T, 
+                          x2::T, 
+                          y2::T, 
+                          x::T, 
+                          y::T
+                          )::Tuple{T,T,T} where {T <: AbstractFloat}
+
+Finds the nearest position along a straight line to a given point.
+
+# Arguments
+- `x1::T`, `y1::T`: Starting point of the line.
+- `x2::T`, `y2::T`: Ending point of the line.
+- `x::T`, `y::T`: Point to nearest position to.
+
+# Returns
+- `::Tuple`:
+  - `::T`: x-coordinate of nearest position.
+  - `::T`: y-coordinate of nearest position.
+  - `::T`: Position along the line, from 0 to 1.
+"""
+function nearest_point_on_line(x1::T, 
+                               y1::T, 
+                               x2::T, 
+                               y2::T, 
+                               x::T, 
+                               y::T
+                               )::Tuple{T,T,T} where {T <: AbstractFloat}
+    A = x - x1
+    B = y - y1
+    C = x2 - x1
+    D = y2 - y1
+    dot = A * C + B * D
+    len_sq = C * C + D * D
+    param = -one(T)
+    if len_sq != 0 # in case of 0 length line
+        param = dot / len_sq
+    end
+    if param < 0.0
+        return (x1, y1, zero(T))
+    elseif param > 1.0
+        return (x2, y2, one(T))
+    else
+        return (x1 + param * C, y1 + param * D, param)
+    end
+end

src/graph.jl

@@ -55,7 +55,7 @@ function graph_from_object(osm_data_object::Union{XMLDocument,Dict};
         g.dijkstra_states = Vector{Vector{U}}(undef, length(g.nodes))
     end
 
-    add_kdtree!(g)
+    add_kdtree_and_rtree!(g)
     @info "Created OSMGraph object with kwargs: network_type=$network_type, weight_type=$weight_type, graph_type=$graph_type, precompute_dijkstra_states=$precompute_dijkstra_states, largest_connected_component=$largest_connected_component"
     return g
 end
@@ -454,16 +454,76 @@ function add_dijkstra_states!(g::OSMGraph{U,T,W}) where {U <: Integer,T <: Integ
     set_dijkstra_state!(g, collect(vertices(g.graph)))
 end
 
+"""
+    get_cartesian_locations(g::OSMGraph)
+
+Calculates the Cartesian location of all nodes in the graph.
+
+Returns a 3-by-n matrix where each column is the `xyz` coordinates of a node. Column indices 
+correspond to the `g.graph` vertex indices.
+"""
+function get_cartesian_locations(g::OSMGraph)
+    node_locations = [index_to_node(g, index).location for index in 1:nv(g.graph)]
+    return to_cartesian(node_locations)
+end
+
+"""
+    add_kdtree_and_rtree!(g::OSMGraph)
+
+Adds k-d tree and R-tree to `OSMGraph` for finding nearest nodes and ways.
+"""
+function add_kdtree_and_rtree!(g::OSMGraph)
+    cartesian_locations = get_cartesian_locations(g)
+    add_kdtree!(g, cartesian_locations)
+    add_rtree!(g, cartesian_locations)
+end
+
 """
     add_kdtree!(g::OSMGraph)
+    add_kdtree!(g::OSMGraph, cartesian_locations::Matrix{Float64})
 
 Adds KDTree to `OSMGraph` for finding nearest neighbours.
 """
-function add_kdtree!(g::OSMGraph)
-    node_locations = [node.location for (id, node) in g.nodes]  # node locations must have the same order as node indices
-    cartesian_locations = to_cartesian(node_locations)
+function add_kdtree!(g::OSMGraph, cartesian_locations::Matrix{Float64})
     g.kdtree = KDTree(cartesian_locations)
 end
+function add_kdtree!(g::OSMGraph)
+    cartesian_locations = get_cartesian_locations(g)
+    add_kdtree!(g, cartesian_locations)
+end
+
+""" 
+    add_rtree!(g::OSMGraph)
+    add_rtree!(g::OSMGraph, cartesian_locations::Matrix{Float64})
+
+Adds an R-tree to `OSMGraph` for finding nearest ways.
+
+# Warning
+Make sure to suppress outputs! 
+Behaviour as of SpatialIndexing.jl 0.1.3 will print a line for every single OSM way, which 
+will flood the terminal if not suppressed. Use with caution for now.
+"""
+function add_rtree!(g::OSMGraph{U,T,W}, cartesian_locations::Matrix{Float64}) where {U,T,W}
+    # Get bounding box for every way ID
+    way_ids = collect(keys(g.ways))
+    data = map(way_ids) do way_id
+        node_indices = node_id_to_index(g, g.ways[way_id].nodes)
+        x = [cartesian_locations[1,i] for i in node_indices]
+        y = [cartesian_locations[2,i] for i in node_indices]
+        z = [cartesian_locations[3,i] for i in node_indices]
+        min_pt = (minimum(x), minimum(y), minimum(z))
+        max_pt = (maximum(x), maximum(y), maximum(z))
+        return SpatialElem(SpatialIndexing.Rect(min_pt, max_pt), way_id, nothing)
+    end
+
+    tree = RTree{Float64,3}(T, Nothing)
+    SpatialIndexing.load!(tree, data)
+    g.rtree = tree
+end
+function add_rtree!(g::OSMGraph)
+    cartesian_locations = get_cartesian_locations(g)
+    add_rtree!(g, cartesian_locations)
+end
 
 """
     get_graph_type(g::OSMGraph)