Merge pull request #61 from tedsteiner/geodesy

Move geodesy code to an independent package

REQUIRE

@@ -1,6 +1,8 @@
-julia 0.3-
+julia 0.3
 Compat
-LightXML
+Reexport
+Geodesy
+Graphs
 LibExpat
+LightXML
 Winston
-Graphs

docs/examples.rst

@@ -16,7 +16,7 @@ Define map boundary and crop:
 
 .. code-block:: python
 
-    bounds = OpenStreetMap.Bounds(42.365, 42.3675, -71.1, -71.094)
+    bounds = Bounds(42.365, 42.3675, -71.1, -71.094)
 
     cropMap!(nodes, bounds, highways=hwys, buildings=builds, features=feats, delete_nodes=false)
 
@@ -42,9 +42,9 @@ Convert map nodes to ENU coordinates:
 
 .. code-block:: python
 
-    reference = OpenStreetMap.centerBounds(bounds)
-    nodesENU = lla2enu(nodes, reference)
-    boundsENU = lla2enu(bounds, reference)
+    reference = center(bounds)
+    nodesENU = ENU(nodes, reference)
+    boundsENU = ENU(bounds, reference)
 
 Create transportation network:
 
@@ -58,8 +58,8 @@ Route planning:
 
 .. code-block:: python
 
-    loc_start = OpenStreetMap.ENU(-5000, 5500, 0)
-    loc_end = OpenStreetMap.ENU(5500, -4000, 0)
+    loc_start = ENU(-5000, 5500, 0)
+    loc_end = ENU(5500, -4000, 0)
 
     node0 = nearestNode(nodesENU, loc_start, network.v_inv)
     node1 = nearestNode(nodesENU, loc_end, network.v_inv)

docs/workers.rst

@@ -46,10 +46,10 @@ Converting Map Coordinate Systems
 
 OpenStreetMap.jl is capable of converting map data between LLA, ECEF, and ENU coordinates (see "Data Types") for definitions of these standard coordinates. Because point location data is ONLY stored in the ``nodes`` dictionary (type ``Dict{Int,Point-Type}``), only this object needs to be converted. Note that Bounds objects also need to be converted, although they don't technically store map data. The following functions can be used to convert between coordinate systems:
 
-* ``lla2ecef(nodes::Dict{Int,LLA})``
-* ``ecef2lla(nodes::Dict{Int,ECEF})``
-* ``ecef2enu(nodes::Dict{Int,ECEF}, reference::LLA)``
-* ``lla2enu(nodes::Dict{Int,LLA}, reference::LLA)``
+* ``ECEF(nodes::Dict{Int,LLA})``
+* ``LLA(nodes::Dict{Int,ECEF})``
+* ``ENU(nodes::Dict{Int,ECEF}, reference::LLA)``
+* ``ENU(nodes::Dict{Int,LLA}, reference::LLA)``
 
 East-North-Up coordinates require an additional input parameter, ``reference``, which gives the origin of the ENU coordinate system. LLA and ECEF coordinates both have their origins fixed at the center of the earth.
 

src/OpenStreetMap.jl

@@ -6,18 +6,19 @@
 
 module OpenStreetMap
 
-import LightXML
-import LibExpat
-import Winston
-import Graphs
-import Compat
+using Reexport
+@reexport using Geodesy
+using LightXML
+using LibExpat
+using Winston
+using Graphs
+using Compat
 
 export parseMapXML, getOSMData, getBounds
 export plotMap, cropMap!
 export findIntersections, nearestNode, segmentHighways, highwaySegments
-export lla2enu, lla2ecef, ecef2lla, ecef2enu
 export roadways, walkways, cycleways, classify
-export createGraph, shortestRoute, fastestRoute, distance, routeEdges
+export createGraph, shortestRoute, fastestRoute, routeEdges
 export nodesWithinRange, nodesWithinDrivingDistance, nodesWithinDrivingTime
 export findHighwaySets, findIntersectionClusters, replaceHighwayNodes!
 export simCityGrid
@@ -42,4 +43,6 @@ include("routing.jl")
 
 include("simulate.jl")
 
+include("deprecated.jl")
+
 end # module OpenStreetMap

src/crop.jl

@@ -105,53 +105,6 @@ function crop!(nodes::Dict, bounds::Bounds, features::Dict{Int,Feature})
     return nothing
 end
 
-### Check whether a location is within bounds ###
-function inBounds{T<:Union(LLA,ENU)}(loc::T, bounds::Bounds{T})
-    x, y = getX(loc), getY(loc)
-
-    bounds.min_x <= x <= bounds.max_x &&
-    bounds.min_y <= y <= bounds.max_y
-end
-
-function onBounds{T<:Union(LLA,ENU)}(loc::T, bounds::Bounds{T})
-    x, y = getX(loc), getY(loc)
-
-    x == bounds.min_x || x == bounds.max_x ||
-    y == bounds.min_y || y == bounds.max_y
-end
-
-function boundaryPoint{T<:Union(LLA,ENU)}(p1::T, p2::T, bounds::Bounds)
-    x1, y1 = getX(p1), getY(p1)
-    x2, y2 = getX(p2), getY(p2)
-
-    x, y = x1, y1
-
-    # checks assume inBounds(p1) != inBounds(p2)
-    if x1 < bounds.min_x < x2 || x1 > bounds.min_x > x2
-        x = bounds.min_x
-        y = y1 + (y2 - y1) * (bounds.min_x - x1) / (x2 - x1)
-    elseif x1 < bounds.max_x < x2 || x1 > bounds.max_x > x2
-        x = bounds.max_x
-        y = y1 + (y2 - y1) * (bounds.max_x - x1) / (x2 - x1)
-    end
-
-    p3 = T(XY(x, y))
-    inBounds(p3, bounds) && return p3
-
-    if y1 < bounds.min_y < y2 || y1 > bounds.min_y > y2
-        x = x1 + (x2 - x1) * (bounds.min_y - y1) / (y2 - y1)
-        y = bounds.min_y
-    elseif y1 < bounds.max_y < y2 || y1 > bounds.max_y > y2
-        x = x1 + (x2 - x1) * (bounds.max_y - y1) / (y2 - y1)
-        y = bounds.max_y
-    end
-
-    p3 = T(XY(x, y))
-    inBounds(p3, bounds) && return p3
-
-    error("Failed to find boundary point.")
-end
-
 function cropHighway!(nodes::Dict, bounds::Bounds, highway::Highway, valids::BitVector)
     prev_id, prev_valid = highway.nodes[1], valids[1]
     ni = 1
@@ -164,8 +117,9 @@ function cropHighway!(nodes::Dict, bounds::Bounds, highway::Highway, valids::Bit
         end
         if valid != prev_valid
             prev_node, node = nodes[prev_id], nodes[id]
-            if !(onBounds(prev_node, bounds) || onBounds(node, bounds))
-                new_node = boundaryPoint(prev_node, node, bounds)
+            if !(Geodesy.onBounds(prev_node, bounds) ||
+                 Geodesy.onBounds(node, bounds))
+                new_node = Geodesy.boundaryPoint(prev_node, node, bounds)
                 new_id = addNewNode(nodes, new_node)
                 insert!(highway.nodes, ni + !valid, new_id)
                 ni += 1

src/deprecated.jl

@@ -0,0 +1,5 @@
+Base.@deprecate lla2ecef ECEF
+Base.@deprecate ecef2lla LLA
+Base.@deprecate ecef2enu ENU
+Base.@deprecate lla2enu ENU
+Base.@deprecate centerBounds(bounds::Bounds) center(bounds)

src/parseMap.jl

@@ -181,7 +181,7 @@ function collectValues(handler::LibExpat.XPStreamHandler, name::String)
             if attr.oneway_reverse
                 reverse!(attr.way_nodes)
             end
-            osm.highways[attr.id] = Highway(attr.class, attr.lanes, 
+            osm.highways[attr.id] = Highway(attr.class, attr.lanes,
                                             (attr.oneway && !attr.oneway_override),
                                             attr.sidewalk, attr.cycleway, attr.bicycle,
                                             attr.name, copy(attr.way_nodes))
@@ -207,7 +207,7 @@ end
 
 function getOSMData(filename::String; args...)
     osm = OSMdata()
-    
+
     callbacks = LibExpat.XPCallbacks()
     callbacks.start_element = parseElement
     callbacks.end_element = collectValues

src/routing.jl

@@ -143,55 +143,21 @@ end
 
 ### Get distance between two nodes ###
 # ENU Coordinates
-function distance(nodes::Dict{Int,ENU}, node0, node1)
+function Geodesy.distance{T<:Union(ENU,ECEF)}(nodes::Dict{Int,T}, node0, node1)
     loc0 = nodes[node0]
     loc1 = nodes[node1]
 
     return distance(loc0, loc1)
 end
 
-function distance(loc0::ENU, loc1::ENU)
-    x0 = loc0.east
-    y0 = loc0.north
-    z0 = loc0.up
-
-    x1 = loc1.east
-    y1 = loc1.north
-    z1 = loc1.up
-
-    return distance(x0, y0, z0, x1, y1, z1)
-end
-
-# ECEF Coordinates
-function distance(nodes::Dict{Int,ECEF}, node0, node1)
-    loc0 = nodes[node0]
-    loc1 = nodes[node1]
-
-    return distance(loc0, loc1)
-end
-
-function distance(loc0::ECEF, loc1::ECEF)
-    x0 = loc0.x
-    y0 = loc0.y
-    z0 = loc0.z
-
-    x1 = loc1.x
-    y1 = loc1.y
-    z1 = loc1.z
-
-    return distance(x0, y0, z0, x1, y1, z1)
-end
-
-# Cartesian coordinates
-function distance(x0, y0, z0, x1, y1, z1)
-    return sqrt((x1-x0)^2 + (y1-y0)^2 + (z1-z0)^2)
-end
-
 ### Compute the distance of a route ###
-function distance(nodes, route)
-    dist = 0
-    for n = 2:length(route)
-        dist += distance(nodes, route[n-1], route[n])
+function Geodesy.distance{T<:Union(ENU,ECEF)}(nodes::Dict{Int,T}, route::Vector{Int})
+    dist = 0.0
+    prev_point = nodes[route[1]]
+    for i = 2:length(route)
+        point = nodes[route[i]]
+        dist += distance(prev_point, point)
+        prev_point = point
     end
 
     return dist