1 // Boost.Geometry (aka GGL, Generic Geometry Library)
3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
4 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
5 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
7 // Use, modification and distribution is subject to the Boost Software License,
8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
9 // http://www.boost.org/LICENSE_1_0.txt)
11 // Example showing Boost.Geometry combined with Boost.Graph, calculating shortest routes
12 // input: two WKT's, provided in subfolder data
13 // output: text, + an SVG, displayable in e.g. Firefox)
16 // Second variant, with bundled properties
24 #include <boost/tuple/tuple.hpp>
25 #include <boost/foreach.hpp>
27 #include <boost/graph/adjacency_list.hpp>
28 #include <boost/graph/dijkstra_shortest_paths.hpp>
30 #include <boost/geometry/geometry.hpp>
31 #include <boost/geometry/geometries/linestring.hpp>
32 #include <boost/geometry/io/wkt/read.hpp>
35 #include <boost/geometry/io/svg/svg_mapper.hpp>
37 // For distance-calculations over the Earth:
38 //#include <boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp>
42 // Read an ASCII file containing WKT's, fill a vector of tuples
43 // The tuples consist of at least <0> a geometry and <1> an identifying string
44 template <typename Geometry
, typename Tuple
, typename Box
>
45 void read_wkt(std::string
const& filename
, std::vector
<Tuple
>& tuples
, Box
& box
)
47 std::ifstream
cpp_file(filename
.c_str());
48 if (cpp_file
.is_open())
50 while (! cpp_file
.eof() )
53 std::getline(cpp_file
, line
);
56 if (! line
.empty() && ! boost::starts_with(line
, "#"))
60 // Split at ';', if any
61 std::string::size_type pos
= line
.find(";");
62 if (pos
!= std::string::npos
)
64 name
= line
.substr(pos
+ 1);
72 boost::geometry::read_wkt(line
, geometry
);
74 Tuple
tuple(geometry
, name
);
76 tuples
.push_back(tuple
);
77 boost::geometry::expand(box
, boost::geometry::return_envelope
<Box
>(geometry
));
83 // To calculate distance, declare and construct a strategy with average earth radius
84 boost::geometry::strategy::distance::haversine
<double> const haversine(6372795.0);
86 // Define properties for vertex
87 template <typename Point
>
88 struct bg_vertex_property
92 boost::geometry::assign_zero(location
);
94 bg_vertex_property(Point
const& loc
)
102 // Define properties for edge
103 template <typename Linestring
>
104 struct bg_edge_property
106 bg_edge_property(Linestring
const& line
)
107 : length(boost::geometry::length(line
, haversine
))
112 inline Linestring
const& line() const
122 // Utility function to add a vertex to a graph. It might exist already. Then do not insert,
123 // but return vertex descriptor back. It might not exist. Then add it (and return).
124 // To efficiently handle this, a std::map is used.
125 template <typename M
, typename K
, typename G
>
126 inline typename
boost::graph_traits
<G
>::vertex_descriptor
find_or_insert(M
& map
, K
const& key
, G
& graph
)
128 typename
M::const_iterator it
= map
.find(key
);
131 // Add a vertex to the graph
132 typename
boost::graph_traits
<G
>::vertex_descriptor new_vertex
133 = boost::add_vertex(bg_vertex_property
<typename
M::key_type
>(key
), graph
);
135 // Add to the map, using POINT as key
136 map
[key
] = new_vertex
;
145 typename RoadTupleVector
,
146 typename CityTupleVector
148 void add_roads_and_connect_cities(Graph
& graph
,
149 RoadTupleVector
const& roads
,
150 CityTupleVector
& cities
)
152 typedef typename
boost::range_value
<RoadTupleVector
>::type road_type
;
153 typedef typename
boost::tuples::element
<0, road_type
>::type line_type
;
154 typedef typename
boost::geometry::point_type
<line_type
>::type point_type
;
156 typedef typename
boost::graph_traits
<Graph
>::vertex_descriptor vertex_type
;
158 // Define a map to be used during graph filling
159 // Maps from point to vertex-id's
160 typedef std::map
<point_type
, vertex_type
, boost::geometry::less
<point_type
> > map_type
;
165 BOOST_FOREACH(road_type
const& road
, roads
)
167 line_type
const& line
= road
.template get
<0>();
168 // Find or add begin/end point of these line
169 vertex_type from
= find_or_insert(map
, line
.front(), graph
);
170 vertex_type to
= find_or_insert(map
, line
.back(), graph
);
171 boost::add_edge(from
, to
, bg_edge_property
<line_type
>(line
), graph
);
174 // Find nearest graph vertex for each city, using the map
175 typedef typename
boost::range_value
<CityTupleVector
>::type city_type
;
176 BOOST_FOREACH(city_type
& city
, cities
)
178 double min_distance
= 1e300
;
179 for(typename
map_type::const_iterator it
= map
.begin(); it
!= map
.end(); ++it
)
181 double dist
= boost::geometry::distance(it
->first
, city
.template get
<0>());
182 if (dist
< min_distance
)
186 city
.template get
<2>() = it
->second
;
192 template <typename Graph
, typename Route
>
193 inline void add_edge_to_route(Graph
const& graph
,
194 typename
boost::graph_traits
<Graph
>::vertex_descriptor vertex1
,
195 typename
boost::graph_traits
<Graph
>::vertex_descriptor vertex2
,
200 typename
boost::graph_traits
<Graph
>::edge_descriptor
,
202 > opt_edge
= boost::edge(vertex1
, vertex2
, graph
);
205 // Get properties of edge and of vertex
206 bg_edge_property
<Route
> const& edge_prop
= graph
[opt_edge
.first
];
209 typename
boost::geometry::point_type
<Route
>::type
210 > const& vertex_prop
= graph
[vertex2
];
212 // Depending on how edge connects to vertex, copy it forward or backward
213 if (boost::geometry::equals(edge_prop
.line().front(), vertex_prop
.location
))
215 std::copy(edge_prop
.line().begin(), edge_prop
.line().end(),
216 std::back_inserter(route
));
220 std::reverse_copy(edge_prop
.line().begin(), edge_prop
.line().end(),
221 std::back_inserter(route
));
227 template <typename Graph
, typename Route
>
228 inline void build_route(Graph
const& graph
,
229 std::vector
<typename
boost::graph_traits
<Graph
>::vertex_descriptor
> const& predecessors
,
230 typename
boost::graph_traits
<Graph
>::vertex_descriptor vertex1
,
231 typename
boost::graph_traits
<Graph
>::vertex_descriptor vertex2
,
234 typedef typename
boost::graph_traits
<Graph
>::vertex_descriptor vertex_type
;
235 vertex_type pred
= predecessors
[vertex2
];
237 add_edge_to_route(graph
, vertex2
, pred
, route
);
238 while (pred
!= vertex1
)
240 add_edge_to_route(graph
, predecessors
[pred
], pred
, route
);
241 pred
= predecessors
[pred
];
248 // Define a point in the Geographic coordinate system (currently spherical-equatorial)
249 // (geographic calculations are in an extension; for sample it makes no difference)
250 typedef boost::geometry::model::point
252 double, 2, boost::geometry::cs::spherical_equatorial
<boost::geometry::degree
>
255 typedef boost::geometry::model::linestring
<point_type
> line_type
;
257 // Define the graph, lateron containing the road network
258 // With this, specify bundled properties for vertex and edge,
259 // as in http://www.boost.org/doc/libs/1_43_0/libs/graph/doc/bundles.html
260 typedef boost::adjacency_list
262 boost::vecS
, boost::vecS
, boost::undirectedS
263 , bg_vertex_property
<point_type
> // bundled
264 , bg_edge_property
<line_type
>
267 typedef boost::graph_traits
<graph_type
>::vertex_descriptor vertex_type
;
270 // Init a bounding box, lateron used to define SVG map
271 boost::geometry::model::box
<point_type
> box
;
272 boost::geometry::assign_inverse(box
);
277 typedef boost::tuple
<point_type
, std::string
, vertex_type
> city_type
;
278 std::vector
<city_type
> cities
;
279 read_wkt
<point_type
>("data/cities.wkt", cities
, box
);
281 // Read the road network
282 typedef boost::tuple
<line_type
, std::string
> road_type
;
283 std::vector
<road_type
> roads
;
284 read_wkt
<line_type
>("data/roads.wkt", roads
, box
);
287 // Add roads and connect cities
288 add_roads_and_connect_cities(graph
, roads
, cities
);
290 double const km
= 1000.0;
291 std::cout
<< "distances, all in KM" << std::endl
292 << std::fixed
<< std::setprecision(0);
294 // Main functionality: calculate shortest routes from/to all cities
296 // For the first one, the complete route is stored as a linestring
300 int const n
= boost::num_vertices(graph
);
301 BOOST_FOREACH(city_type
const& city1
, cities
)
303 std::vector
<vertex_type
> predecessors(n
);
304 std::vector
<double> costs(n
);
306 // Call Dijkstra (without named-parameter to be compatible with all VC)
307 boost::dijkstra_shortest_paths(graph
, city1
.get
<2>(),
308 &predecessors
[0], &costs
[0],
309 boost::get(&bg_edge_property
<line_type
>::length
, graph
),
310 boost::get(boost::vertex_index
, graph
),
311 std::less
<double>(), std::plus
<double>(),
312 (std::numeric_limits
<double>::max
)(), double(),
313 boost::dijkstra_visitor
<boost::null_visitor
>());
315 BOOST_FOREACH(city_type
const& city2
, cities
)
317 if (! boost::equals(city1
.get
<1>(), city2
.get
<1>()))
319 double distance
= costs
[city2
.get
<2>()] / km
;
320 double acof
= boost::geometry::distance(city1
.get
<0>(), city2
.get
<0>(), haversine
) / km
;
323 << std::setiosflags (std::ios_base::left
) << std::setw(15)
324 << city1
.get
<1>() << " - "
325 << std::setiosflags (std::ios_base::left
) << std::setw(15)
327 << " -> through the air: " << std::setw(4) << acof
328 << " , over the road: " << std::setw(4) << distance
333 build_route(graph
, predecessors
,
334 city1
.get
<2>(), city2
.get
<2>(),
342 #if defined(HAVE_SVG)
344 std::ofstream
stream("routes.svg");
345 boost::geometry::svg_mapper
<point_type
> mapper(stream
, 600, 600);
348 BOOST_FOREACH(road_type
const& road
, roads
)
350 mapper
.add(road
.get
<0>());
353 BOOST_FOREACH(road_type
const& road
, roads
)
355 mapper
.map(road
.get
<0>(),
356 "stroke:rgb(128,128,128);stroke-width:1");
360 "stroke:rgb(0, 255, 0);stroke-width:6;opacity:0.5");
363 BOOST_FOREACH(city_type
const& city
, cities
)
365 mapper
.map(city
.get
<0>(),
366 "fill:rgb(255,255,0);stroke:rgb(0,0,0);stroke-width:1");
367 mapper
.text(city
.get
<0>(), city
.get
<1>(),
368 "fill:rgb(0,0,0);font-family:Arial;font-size:10px", 5, 5);