1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2
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.
6
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)
10 //
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)
14
15 #include <iostream>
16 #include <fstream>
17 #include <iomanip>
18 #include <limits>
19
20 #include <boost/tuple/tuple.hpp>
21 #include <boost/foreach.hpp>
22
23 #include <boost/graph/adjacency_list.hpp>
24 #include <boost/graph/dijkstra_shortest_paths.hpp>
25
26 #include <boost/geometry/geometry.hpp>
27 #include <boost/geometry/geometries/linestring.hpp>
28 #include <boost/geometry/io/wkt/read.hpp>
29
30
31 // For output:
32 #include <boost/geometry/io/svg/svg_mapper.hpp>
33
34 // For distance-calculations over the Earth:
35 //#include <boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp>
36
37 // Read an ASCII file containing WKT's, fill a vector of tuples
38 // The tuples consist of at least <0> a geometry and <1> an identifying string
39 template <typename Geometry, typename Tuple, typename Box>
read_wkt(std::string const & filename,std::vector<Tuple> & tuples,Box & box)40 void read_wkt(std::string const& filename, std::vector<Tuple>& tuples, Box& box)
41 {
42 std::ifstream cpp_file(filename.c_str());
43 if (cpp_file.is_open())
44 {
45 while (! cpp_file.eof() )
46 {
47 std::string line;
48 std::getline(cpp_file, line);
49 Geometry geometry;
50 boost::trim(line);
51 if (! line.empty() && ! boost::starts_with(line, "#"))
52 {
53 std::string name;
54
55 // Split at ';', if any
56 std::string::size_type pos = line.find(";");
57 if (pos != std::string::npos)
58 {
59 name = line.substr(pos + 1);
60 line.erase(pos);
61
62 boost::trim(line);
63 boost::trim(name);
64 }
65
66 Geometry geometry;
67 boost::geometry::read_wkt(line, geometry);
68
69 Tuple tuple(geometry, name);
70
71 tuples.push_back(tuple);
72 boost::geometry::expand(box, boost::geometry::return_envelope<Box>(geometry));
73 }
74 }
75 }
76 }
77
78
79
80 // Code to define properties for Boost Graph's
81 enum vertex_bg_property_t { vertex_bg_property };
82 enum edge_bg_property_t { edge_bg_property };
83 namespace boost
84 {
85 BOOST_INSTALL_PROPERTY(vertex, bg_property);
86 BOOST_INSTALL_PROPERTY(edge, bg_property);
87 }
88
89 // To calculate distance, declare and construct a strategy with average earth radius
90 boost::geometry::strategy::distance::haversine<double> const haversine(6372795.0);
91
92 // Define properties for vertex
93 template <typename Point>
94 struct bg_vertex_property
95 {
bg_vertex_propertybg_vertex_property96 bg_vertex_property()
97 {
98 boost::geometry::assign_zero(location);
99 }
bg_vertex_propertybg_vertex_property100 bg_vertex_property(Point const& loc)
101 : location(loc)
102 {
103 }
104
105 Point location;
106 };
107
108 // Define properties for edge
109 template <typename Linestring>
110 struct bg_edge_property
111 {
bg_edge_propertybg_edge_property112 bg_edge_property(Linestring const& line)
113 : m_length(boost::geometry::length(line, haversine))
114 , m_line(line)
115 {
116 }
117
operator doublebg_edge_property118 inline operator double() const
119 {
120 return m_length;
121 }
122
linebg_edge_property123 inline Linestring const& line() const
124 {
125 return m_line;
126 }
127
128 private :
129 double m_length;
130 Linestring m_line;
131 };
132
133 // Utility function to add a vertex to a graph. It might exist already. Then do not insert,
134 // but return vertex descriptor back. It might not exist. Then add it (and return).
135 // To efficiently handle this, a std::map is used.
136 template <typename M, typename K, typename G>
find_or_insert(M & map,K const & key,G & graph)137 inline typename boost::graph_traits<G>::vertex_descriptor find_or_insert(M& map, K const& key, G& graph)
138 {
139 typename M::const_iterator it = map.find(key);
140 if (it == map.end())
141 {
142 // Add a vertex to the graph
143 typename boost::graph_traits<G>::vertex_descriptor new_vertex
144 = boost::add_vertex(graph);
145
146 // Set the property (= location)
147 boost::put(boost::get(vertex_bg_property, graph), new_vertex,
148 bg_vertex_property<typename M::key_type>(key));
149
150 // Add to the map, using POINT as key
151 map[key] = new_vertex;
152 return new_vertex;
153 }
154 return it->second;
155 }
156
157 template
158 <
159 typename Graph,
160 typename RoadTupleVector,
161 typename CityTupleVector
162 >
add_roads_and_connect_cities(Graph & graph,RoadTupleVector const & roads,CityTupleVector & cities)163 void add_roads_and_connect_cities(Graph& graph,
164 RoadTupleVector const& roads,
165 CityTupleVector& cities)
166 {
167 typedef typename boost::range_value<RoadTupleVector>::type road_type;
168 typedef typename boost::tuples::element<0, road_type>::type line_type;
169 typedef typename boost::geometry::point_type<line_type>::type point_type;
170
171 typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_type;
172
173 // Define a map to be used during graph filling
174 // Maps from point to vertex-id's
175 typedef std::map<point_type, vertex_type, boost::geometry::less<point_type> > map_type;
176 map_type map;
177
178
179 // Fill the graph
180 BOOST_FOREACH(road_type const& road, roads)
181 {
182 line_type const& line = road.template get<0>();
183 // Find or add begin/end point of these line
184 vertex_type from = find_or_insert(map, line.front(), graph);
185 vertex_type to = find_or_insert(map, line.back(), graph);
186 boost::add_edge(from, to, bg_edge_property<line_type>(line), graph);
187 }
188
189 // Find nearest graph vertex for each city, using the map
190 typedef typename boost::range_value<CityTupleVector>::type city_type;
191 BOOST_FOREACH(city_type& city, cities)
192 {
193 double min_distance = 1e300;
194 for(typename map_type::const_iterator it = map.begin(); it != map.end(); ++it)
195 {
196 double dist = boost::geometry::distance(it->first, city.template get<0>());
197 if (dist < min_distance)
198 {
199 min_distance = dist;
200 // Set the vertex
201 city.template get<2>() = it->second;
202 }
203 }
204 }
205 }
206
207 template <typename Graph, typename Route>
add_edge_to_route(Graph const & graph,typename boost::graph_traits<Graph>::vertex_descriptor vertex1,typename boost::graph_traits<Graph>::vertex_descriptor vertex2,Route & route)208 inline void add_edge_to_route(Graph const& graph,
209 typename boost::graph_traits<Graph>::vertex_descriptor vertex1,
210 typename boost::graph_traits<Graph>::vertex_descriptor vertex2,
211 Route& route)
212 {
213 std::pair
214 <
215 typename boost::graph_traits<Graph>::edge_descriptor,
216 bool
217 > opt_edge = boost::edge(vertex1, vertex2, graph);
218 if (opt_edge.second)
219 {
220 // Get properties of edge and of vertex
221 bg_edge_property<Route> const& edge_prop =
222 boost::get(boost::get(edge_bg_property, graph), opt_edge.first);
223
224 bg_vertex_property<typename boost::geometry::point_type<Route>::type> const& vertex_prop =
225 boost::get(boost::get(vertex_bg_property, graph), vertex2);
226
227 // Depending on how edge connects to vertex, copy it forward or backward
228 if (boost::geometry::equals(edge_prop.line().front(), vertex_prop.location))
229 {
230 std::copy(edge_prop.line().begin(), edge_prop.line().end(),
231 std::back_inserter(route));
232 }
233 else
234 {
235 std::reverse_copy(edge_prop.line().begin(), edge_prop.line().end(),
236 std::back_inserter(route));
237 }
238 }
239 }
240
241
242 template <typename Graph, typename Route>
build_route(Graph const & graph,std::vector<typename boost::graph_traits<Graph>::vertex_descriptor> const & predecessors,typename boost::graph_traits<Graph>::vertex_descriptor vertex1,typename boost::graph_traits<Graph>::vertex_descriptor vertex2,Route & route)243 inline void build_route(Graph const& graph,
244 std::vector<typename boost::graph_traits<Graph>::vertex_descriptor> const& predecessors,
245 typename boost::graph_traits<Graph>::vertex_descriptor vertex1,
246 typename boost::graph_traits<Graph>::vertex_descriptor vertex2,
247 Route& route)
248 {
249 typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_type;
250 vertex_type pred = predecessors[vertex2];
251
252 add_edge_to_route(graph, vertex2, pred, route);
253 while (pred != vertex1)
254 {
255 add_edge_to_route(graph, predecessors[pred], pred, route);
256 pred = predecessors[pred];
257 }
258 }
259
260
main()261 int main()
262 {
263 // Define a point in the Geographic coordinate system (currently Spherical)
264 // (geographic calculations are in an extension; for sample it makes no difference)
265 typedef boost::geometry::model::point
266 <
267 double, 2, boost::geometry::cs::spherical_equatorial<boost::geometry::degree>
268 > point_type;
269
270 typedef boost::geometry::model::linestring<point_type> line_type;
271
272 // Define the graph, lateron containing the road network
273 typedef boost::adjacency_list
274 <
275 boost::vecS, boost::vecS, boost::undirectedS
276 , boost::property<vertex_bg_property_t, bg_vertex_property<point_type> >
277 , boost::property<edge_bg_property_t, bg_edge_property<line_type> >
278 > graph_type;
279
280 typedef boost::graph_traits<graph_type>::vertex_descriptor vertex_type;
281
282
283 // Init a bounding box, lateron used to define SVG map
284 boost::geometry::model::box<point_type> box;
285 boost::geometry::assign_inverse(box);
286
287 // Read the cities
288 typedef boost::tuple<point_type, std::string, vertex_type> city_type;
289 std::vector<city_type> cities;
290 read_wkt<point_type>("data/cities.wkt", cities, box);
291
292 // Read the road network
293 typedef boost::tuple<line_type, std::string> road_type;
294 std::vector<road_type> roads;
295 read_wkt<line_type>("data/roads.wkt", roads, box);
296
297
298 graph_type graph;
299
300 // Add roads and connect cities
301 add_roads_and_connect_cities(graph, roads, cities);
302
303 double const km = 1000.0;
304 std::cout << "distances, all in KM" << std::endl
305 << std::fixed << std::setprecision(0);
306
307 // Main functionality: calculate shortest routes from/to all cities
308
309
310 // For the first one, the complete route is stored as a linestring
311 bool first = true;
312 line_type route;
313
314 int const n = boost::num_vertices(graph);
315 BOOST_FOREACH(city_type const& city1, cities)
316 {
317 std::vector<vertex_type> predecessors(n);
318 std::vector<double> costs(n);
319
320 // Call Dijkstra (without named-parameter to be compatible with all VC)
321 boost::dijkstra_shortest_paths(graph, city1.get<2>(),
322 &predecessors[0], &costs[0],
323 boost::get(edge_bg_property, graph),
324 boost::get(boost::vertex_index, graph),
325 std::less<double>(), std::plus<double>(),
326 (std::numeric_limits<double>::max)(), double(),
327 boost::dijkstra_visitor<boost::null_visitor>());
328
329 BOOST_FOREACH(city_type const& city2, cities)
330 {
331 if (! boost::equals(city1.get<1>(), city2.get<1>()))
332 {
333 double distance = costs[city2.get<2>()] / km;
334 double acof = boost::geometry::distance(city1.get<0>(), city2.get<0>(), haversine) / km;
335
336 std::cout
337 << std::setiosflags (std::ios_base::left) << std::setw(15)
338 << city1.get<1>() << " - "
339 << std::setiosflags (std::ios_base::left) << std::setw(15)
340 << city2.get<1>()
341 << " -> through the air: " << std::setw(4) << acof
342 << " , over the road: " << std::setw(4) << distance
343 << std::endl;
344
345 if (first)
346 {
347 build_route(graph, predecessors,
348 city1.get<2>(), city2.get<2>(),
349 route);
350 first = false;
351 }
352 }
353 }
354 }
355
356 #if defined(HAVE_SVG)
357 // Create the SVG
358 std::ofstream stream("routes.svg");
359 boost::geometry::svg_mapper<point_type> mapper(stream, 600, 600);
360
361 // Map roads
362 BOOST_FOREACH(road_type const& road, roads)
363 {
364 mapper.add(road.get<0>());
365 }
366
367 BOOST_FOREACH(road_type const& road, roads)
368 {
369 mapper.map(road.get<0>(),
370 "stroke:rgb(128,128,128);stroke-width:1");
371 }
372
373 mapper.map(route,
374 "stroke:rgb(0, 255, 0);stroke-width:6;opacity:0.5");
375
376 // Map cities
377 BOOST_FOREACH(city_type const& city, cities)
378 {
379 mapper.map(city.get<0>(),
380 "fill:rgb(255,255,0);stroke:rgb(0,0,0);stroke-width:1");
381 mapper.text(city.get<0>(), city.get<1>(),
382 "fill:rgb(0,0,0);font-family:Arial;font-size:10px", 5, 5);
383 }
384 #endif
385
386 return 0;
387 }
388