summaryrefslogtreecommitdiffstats
path: root/server/src
diff options
context:
space:
mode:
Diffstat (limited to 'server/src')
-rw-r--r--server/src/api.cpp86
-rw-r--r--server/src/api.cppm19
-rw-r--r--server/src/datex2.cpp62
-rw-r--r--server/src/datex2.cppm12
-rw-r--r--server/src/geo.cpp42
-rw-r--r--server/src/geo.cppm23
-rw-r--r--server/src/geo/geo.cppm7
-rw-r--r--server/src/geo/multizonal.cppm240
-rw-r--r--server/src/geo/utm.cppm199
-rw-r--r--server/src/geo/utm_zone_local.cppm74
-rw-r--r--server/src/geo/wgs84.cppm115
-rw-r--r--server/src/gpx.cpp16
-rw-r--r--server/src/gpx.cppm4
-rw-r--r--server/src/log.cpp2
14 files changed, 738 insertions, 163 deletions
diff --git a/server/src/api.cpp b/server/src/api.cpp
index 5dd3fc5..137666b 100644
--- a/server/src/api.cpp
+++ b/server/src/api.cpp
@@ -17,12 +17,14 @@ namespace views = std::views;
17 17
18namespace routemon::api { 18namespace routemon::api {
19 19
20auto json_value_from_point(geo::point const& p) -> json::value 20auto json_value_from_point(geo::utm::zonable_wgs84_point const& p)
21 -> json::value
21{ 22{
22 return json::array{bgeo::get<1>(p), bgeo::get<0>(p)}; 23 return json::array{p.lat(), p.lon()};
23} 24}
24 25
25auto json_value_from_linestring(geo::linestring const& ls) -> json::value 26auto json_value_from_linestring(geo::utm::zonable_wgs84_linestring const& ls)
27 -> json::value
26{ 28{
27 json::array a; 29 json::array a;
28 for (auto const& p : ls) 30 for (auto const& p : ls)
@@ -30,8 +32,8 @@ auto json_value_from_linestring(geo::linestring const& ls) -> json::value
30 return a; 32 return a;
31} 33}
32 34
33auto json_value_from_linestrings(std::vector<geo::linestring> const& lss) 35auto json_value_from_linestrings(
34 -> json::value 36 std::vector<geo::utm::zonable_wgs84_linestring> const& lss) -> json::value
35{ 37{
36 json::array a; 38 json::array a;
37 for (auto const& ls : lss) 39 for (auto const& ls : lss)
@@ -108,24 +110,7 @@ handler::handler(log::logger const& l, datex2::situation_publication pub)
108 static_cast<float>(i) / static_cast<float>(nsit) * 100.0f); 110 static_cast<float>(i) / static_cast<float>(nsit) * 100.0f);
109 }; 111 };
110 112
111 l_.info("Building indices for {} situations", nsit); 113 l_.info("Building indexes for {} situations", nsit);
112
113 auto const buffer_distance = min_distance_;
114 auto const points_per_circle = 8;
115 // Note: thomas strategy does not work for geographic_join_round;
116 // need to use andoyer for that.
117 using formula = bgeo::strategy::andoyer;
118 bgeo::strategy::buffer::distance_symmetric<double> distance_strategy{
119 buffer_distance
120 };
121 bgeo::strategy::buffer::geographic_join_round<formula> join_strategy{
122 points_per_circle
123 };
124 bgeo::strategy::buffer::geographic_end_round<formula> end_strategy{4};
125 bgeo::strategy::buffer::geographic_point_circle<formula> circle_strategy{
126 points_per_circle
127 };
128 bgeo::strategy::buffer::geographic_side_straight<formula> side_strategy;
129 114
130 l_.info("Progress: 0%"); 115 l_.info("Progress: 0%");
131 auto const before_build = chrono::steady_clock::now(); 116 auto const before_build = chrono::steady_clock::now();
@@ -135,17 +120,11 @@ handler::handler(log::logger const& l, datex2::situation_publication pub)
135 { 120 {
136 for (auto const& ls : rc->relevant_line_strings) 121 for (auto const& ls : rc->relevant_line_strings)
137 { 122 {
138 auto buffered_ls = geo::multi_polygon{}; 123 blse_index_.insert(*ls, rc);
139 bgeo::buffer(
140 *ls, buffered_ls, distance_strategy, side_strategy, join_strategy,
141 end_strategy, circle_strategy);
142 auto blse = geo::box{};
143 bgeo::envelope(buffered_ls, blse);
144 blse_index_.insert(std::make_tuple(blse, buffered_ls, rc));
145 } 124 }
146 for (auto p : rc->relevant_points) 125 for (auto p : rc->relevant_points)
147 { 126 {
148 p_index_.insert(std::make_pair(p, rc)); 127 p_index_.insert(p_index_value{p, rc});
149 } 128 }
150 } 129 }
151 130
@@ -160,7 +139,7 @@ handler::handler(log::logger const& l, datex2::situation_publication pub)
160 139
161 auto const after_build = chrono::steady_clock::now(); 140 auto const after_build = chrono::steady_clock::now();
162 auto const dur_build_s = chrono::duration<double>{after_build - before_build}; 141 auto const dur_build_s = chrono::duration<double>{after_build - before_build};
163 l_.info("Indices built in {}", dur_build_s); 142 l_.info("Indexes built in {}", dur_build_s);
164 l_.info("BLSE index size: {}", blse_index_.size()); 143 l_.info("BLSE index size: {}", blse_index_.size());
165 l_.info("Point index size: {}", p_index_.size()); 144 l_.info("Point index size: {}", p_index_.size());
166} 145}
@@ -175,16 +154,13 @@ auto handler::process_gpx(gpx::file&& gpx_file)
175 auto const check_periods = time::period_seq{relevant.begin(), relevant.end()}; 154 auto const check_periods = time::period_seq{relevant.begin(), relevant.end()};
176 155
177 // TODO: eliminate use of overlap segments 156 // TODO: eliminate use of overlap segments
178 auto splits_with_overlap_segments = std::vector<geo::linestring>{}; 157 auto segments = std::vector<geo::utm::zonable_wgs84_linestring>{};
179 for (auto const& track : gpx_file.tracks) 158 for (auto const& track : gpx_file.tracks)
180 for (auto const& seg : track.segments) 159 for (auto const& seg : track.segments)
181 geo::split_linestring_with_overlap_segments( 160 segments.push_back(seg.waypoints);
182 seg.waypoints,
183 5000 /* meters max total dist until a new split is forced */,
184 splits_with_overlap_segments);
185 auto const before_query = chrono::steady_clock::now(); 161 auto const before_query = chrono::steady_clock::now();
186 162
187 auto vincenty_strategy = geo::vincenty_strategy{}; 163 auto vincenty_strategy = geo::wgs84::vincenty_strategy{};
188 164
189 l_.debug("Querying for relevant situations"); 165 l_.debug("Querying for relevant situations");
190 auto relevant_road_closures = 166 auto relevant_road_closures =
@@ -192,28 +168,17 @@ auto handler::process_gpx(gpx::file&& gpx_file)
192 auto ls_checked = 0uz; 168 auto ls_checked = 0uz;
193 auto p_checked = 0uz; 169 auto p_checked = 0uz;
194 auto i = 0; 170 auto i = 0;
195 for (geo::linestring const& part : splits_with_overlap_segments) 171 for (auto const& part : segments)
196 { 172 {
197 l_.debug("Checking part [{}/{}]", ++i, splits_with_overlap_segments.size()); 173 l_.debug("Checking part [{}/{}]", ++i, segments.size());
198 174
175 auto part_zone_lss = geo::multizonal::zoned_linestring_seg_seq{part};
199 // TODO: consider buffering with min_distance_ 176 // TODO: consider buffering with min_distance_
200 auto part_box = geo::box{}; 177 auto part_box = geo::wgs84::box{};
201 bgeo::envelope(part, part_box); 178 bgeo::envelope(part, part_box);
202 179
203 for (auto it = blse_index_.qbegin(bgeo::index::intersects(part_box)); 180 for (auto const& rc : blse_index_.intersection(part_zone_lss))
204 it != blse_index_.qend(); it++)
205 { 181 {
206 // Cannot use structured bindings here, as boost::geometry::get
207 // interferes with ADL. It is a candidate as the namespace
208 // boost::geometry is part of the associated namespace set, which
209 // happens because geo::linestring ≡
210 // boost::geometry::model::linestring<geo::point> is part of the
211 // whole tuple type (blse_index_value) that is the value_type of the
212 // iterator.
213
214 geo::box const& blse = std::get<0>(*it);
215 geo::multi_polygon const& bls = std::get<1>(*it);
216 std::shared_ptr<datex2::road_closure> const& rc = std::get<2>(*it);
217 if (rc->validity 182 if (rc->validity
218 && rc->validity->intersect(check_periods).periods().empty()) 183 && rc->validity->intersect(check_periods).periods().empty())
219 continue; 184 continue;
@@ -221,8 +186,7 @@ auto handler::process_gpx(gpx::file&& gpx_file)
221 // auto small_parts = bgeo::model::multi_linestring<geo::linestring>{}; 186 // auto small_parts = bgeo::model::multi_linestring<geo::linestring>{};
222 // bgeo::intersection(part, lse, small_parts); 187 // bgeo::intersection(part, lse, small_parts);
223 188
224 if (bgeo::intersects(blse, part) && bgeo::intersects(bls, part)) 189 relevant_road_closures.emplace(rc);
225 relevant_road_closures.emplace(rc);
226 190
227 // if (!bgeo::is_empty(small_parts)) 191 // if (!bgeo::is_empty(small_parts))
228 // if (bgeo::distance(*ls, small_parts, vincenty_strategy) < 192 // if (bgeo::distance(*ls, small_parts, vincenty_strategy) <
@@ -234,7 +198,7 @@ auto handler::process_gpx(gpx::file&& gpx_file)
234 it != p_index_.qend(); it++) 198 it != p_index_.qend(); it++)
235 { 199 {
236 // Cannot use structured bindings here for the same reason as above. 200 // Cannot use structured bindings here for the same reason as above.
237 geo::point const& p = std::get<0>(*it); 201 geo::utm::zonable_wgs84_point const& p = std::get<0>(*it);
238 std::shared_ptr<datex2::road_closure> const& rc = std::get<1>(*it); 202 std::shared_ptr<datex2::road_closure> const& rc = std::get<1>(*it);
239 if (rc->validity 203 if (rc->validity
240 && rc->validity->intersect(check_periods).periods().empty()) 204 && rc->validity->intersect(check_periods).periods().empty())
@@ -301,10 +265,12 @@ auto handler::process_gpx(gpx::file&& gpx_file)
301 .relevant_lss = 265 .relevant_lss =
302 rc->relevant_line_strings 266 rc->relevant_line_strings
303 | views::transform( 267 | views::transform(
304 [](auto const& lsp) -> geo::linestring 268 [](auto const& lsp)
269 -> geo::utm::zonable_wgs84_linestring
305 { return *lsp; }) 270 { return *lsp; })
306 | std::ranges:: 271 | std::ranges::to<std::vector<
307 to<std::vector<geo::linestring>>(), 272 geo::utm::zonable_wgs84_linestring
273 >>(),
308 }; 274 };
309 }) 275 })
310 | std::ranges::to<std::vector<relevant_road_closure>>(), 276 | std::ranges::to<std::vector<relevant_road_closure>>(),
diff --git a/server/src/api.cppm b/server/src/api.cppm
index b0f3944..77e0f85 100644
--- a/server/src/api.cppm
+++ b/server/src/api.cppm
@@ -7,7 +7,8 @@ export module routemon:api;
7 7
8import std; 8import std;
9import :datex2; 9import :datex2;
10import :geo; 10import :geo.wgs84;
11import :geo.multizonal;
11import :gpx; 12import :gpx;
12import :log; 13import :log;
13import :time; 14import :time;
@@ -25,20 +26,20 @@ export namespace routemon::api {
25 26
26struct relevant_road_closure 27struct relevant_road_closure
27{ 28{
28 std::vector<geo::linestring> relevant_lss; 29 std::vector<geo::utm::zonable_wgs84_linestring> relevant_lss;
29}; 30};
30 31
31struct relevant_situation 32struct relevant_situation
32{ 33{
33 std::string id; 34 std::string id;
34 std::optional<geo::point> location; 35 std::optional<geo::utm::zonable_wgs84_point> location;
35 std::vector<std::string> comments; 36 std::vector<std::string> comments;
36 std::vector<relevant_road_closure> relevant_road_closures; 37 std::vector<relevant_road_closure> relevant_road_closures;
37}; 38};
38 39
39struct track_segment 40struct track_segment
40{ 41{
41 geo::linestring points; 42 geo::utm::zonable_wgs84_linestring points;
42}; 43};
43 44
44struct track 45struct track
@@ -79,13 +80,11 @@ class handler
79{ 80{
80 constexpr static auto const min_distance_ = 5.0; 81 constexpr static auto const min_distance_ = 5.0;
81 82
82 using blse_index_value = std::tuple< 83 using blse_index_value = std::shared_ptr<datex2::road_closure>;
83 geo::box, geo::multi_polygon, std::shared_ptr<datex2::road_closure> 84 using blse_index = geo::multizonal::linestring_rtree<blse_index_value>;
85 using p_index_value = std::pair<
86 geo::utm::zonable_wgs84_point, std::shared_ptr<datex2::road_closure>
84 >; 87 >;
85 using p_index_value =
86 std::pair<geo::point, std::shared_ptr<datex2::road_closure>>;
87 using blse_index =
88 bgeo::index::rtree<blse_index_value, bgeo::index::quadratic<16>>;
89 using p_index = bgeo::index::rtree<p_index_value, bgeo::index::quadratic<16>>; 88 using p_index = bgeo::index::rtree<p_index_value, bgeo::index::quadratic<16>>;
90 89
91 log::logger l_; 90 log::logger l_;
diff --git a/server/src/datex2.cpp b/server/src/datex2.cpp
index 4c8000b..70c5f94 100644
--- a/server/src/datex2.cpp
+++ b/server/src/datex2.cpp
@@ -51,32 +51,41 @@ auto loader::add_location_from_xml(
51 } 51 }
52 auto const pos_list_str = 52 auto const pos_list_str =
53 std::string_view{loc_gml_xml.child_value("loc:posList")}; 53 std::string_view{loc_gml_xml.child_value("loc:posList")};
54 // lat1 long1 lat2 long2 ... lat(n-1) long(n-1) latn longn 54 // lat1 lon1 lat2 lon2 ... lat(n-1) lon(n-1) latn lonn
55 55
56 auto ls = std::make_shared<geo::linestring>(); 56 auto ls = std::make_shared<geo::utm::zonable_wgs84_linestring>();
57 57
58 auto lat_set = false; 58 auto lat_set = false;
59 auto lat = 0.0; 59 auto lat = 0.0;
60 for (auto const lat_or_long_str : std::views::split(pos_list_str, " "sv)) 60 for (auto const lat_or_lon_str : std::views::split(pos_list_str, " "sv))
61 { 61 {
62 auto mlat_or_long = util::parse_double(std::string_view{lat_or_long_str}); 62 auto mlat_or_lon = util::parse_double(std::string_view{lat_or_lon_str});
63 if (!mlat_or_long) 63 if (!mlat_or_lon)
64 { 64 {
65 warnings_.insert( 65 warnings_.insert(
66 std::format( 66 std::format(
67 "failed to parse coordinate {:?}", 67 "failed to parse coordinate {:?}",
68 std::string_view{lat_or_long_str})); 68 std::string_view{lat_or_lon_str}));
69 return; 69 return;
70 } 70 }
71 71
72 if (!lat_set) 72 if (!lat_set)
73 { 73 {
74 lat = *mlat_or_long; 74 lat = *mlat_or_lon;
75 lat_set = true; 75 lat_set = true;
76 } 76 }
77 else 77 else
78 { 78 {
79 bgeo::append(*ls, geo::point{*mlat_or_long, lat}); 79 auto mp = geo::utm::zonable_wgs84_point::from(
80 geo::wgs84::normalized_point{geo::wgs84::from_lat_lon(
81 lat, *mlat_or_lon)});
82 if (!mp)
83 warnings_.insert(
84 std::format(
85 "point ({}, {}) outside of UTM range; ignoring it", lat,
86 *mlat_or_lon));
87 else
88 bgeo::append(*ls, *mp);
80 lat = 0; 89 lat = 0;
81 lat_set = false; 90 lat_set = false;
82 } 91 }
@@ -110,11 +119,24 @@ auto loader::add_location_from_xml(
110 // maar heeft het UML-model van DATEX II v3 het over ETRS 89: 119 // maar heeft het UML-model van DATEX II v3 het over ETRS 89:
111 // https://docs.datex2.eu/_static/data/v3.7/umlmodel/html/EARoot/EA3/EA3/EA5/EA676.htm 120 // https://docs.datex2.eu/_static/data/v3.7/umlmodel/html/EARoot/EA3/EA3/EA5/EA676.htm
112 121
113 auto const coords_etrs89 = geo::point{*mlon, *mlat}; 122 auto const coords_etrs89 = geo::wgs84::point{*mlon, *mlat};
114 auto coords_wgs84 = geo::point{}; 123 auto coords_wgs84 = geo::wgs84::point{};
115 etrs89_to_wgs84_.forward(coords_etrs89, coords_wgs84); 124 etrs89_to_wgs84_.forward(coords_etrs89, coords_wgs84);
116 125
117 rc.relevant_points.push_back(coords_wgs84); 126 auto mp = geo::utm::zonable_wgs84_point::from(
127 geo::wgs84::normalized_point{coords_wgs84});
128 if (!mp)
129 {
130 warnings_.insert(
131 std::format(
132 "(ETRS89) point ({}, {}) outside of UTM range after "
133 "transformation to WGS 84; ignoring it",
134 *mlat, *mlon));
135 }
136 else
137 {
138 rc.relevant_points.push_back(*mp);
139 }
118 } 140 }
119 else 141 else
120 { 142 {
@@ -347,13 +369,23 @@ auto loader::load_situation_publication(std::string const& filename)
347 // het over ETRS 89: 369 // het over ETRS 89:
348 // https://docs.datex2.eu/_static/data/v3.7/umlmodel/html/EARoot/EA3/EA3/EA5/EA676.htm 370 // https://docs.datex2.eu/_static/data/v3.7/umlmodel/html/EARoot/EA3/EA3/EA5/EA676.htm
349 371
350 auto const coords_etrs89 = geo::point{*mlon, *mlat}; 372 auto const coords_etrs89 = geo::wgs84::point{*mlon, *mlat};
351 auto coords_wgs84 = geo::point{}; 373 auto coords_wgs84 = geo::wgs84::point{};
352 etrs89_to_wgs84_.forward(coords_etrs89, coords_wgs84); 374 etrs89_to_wgs84_.forward(coords_etrs89, coords_wgs84);
353 375
354 if (!sit->location) 376 auto mp = geo::utm::zonable_wgs84_point::from(
377 geo::wgs84::normalized_point{coords_wgs84});
378 if (!mp)
379 {
380 warnings_.insert(
381 std::format(
382 "(ETRS89) point ({}, {}) outside of UTM range after "
383 "transformation to WGS 84; ignoring it",
384 *mlat, *mlon));
385 }
386 else if (!sit->location)
355 { 387 {
356 sit->location = coords_wgs84; 388 sit->location = *mp;
357 } 389 }
358 } 390 }
359 } 391 }
diff --git a/server/src/datex2.cppm b/server/src/datex2.cppm
index 9c75013..ae66631 100644
--- a/server/src/datex2.cppm
+++ b/server/src/datex2.cppm
@@ -1,6 +1,5 @@
1module; 1module;
2 2
3#include <boost/geometry/algorithms/is_empty.hpp>
4#include <boost/geometry/srs/epsg.hpp> 3#include <boost/geometry/srs/epsg.hpp>
5#include <boost/geometry/srs/transformation.hpp> 4#include <boost/geometry/srs/transformation.hpp>
6 5
@@ -9,7 +8,7 @@ module;
9export module routemon:datex2; 8export module routemon:datex2;
10 9
11import std; 10import std;
12import :geo; 11import :geo.utm;
13import :time; 12import :time;
14import :util; 13import :util;
15 14
@@ -23,14 +22,15 @@ export struct road_closure
23{ 22{
24 std::weak_ptr<situation> parent; 23 std::weak_ptr<situation> parent;
25 std::optional<time::period_seq> validity; 24 std::optional<time::period_seq> validity;
26 std::vector<geo::point> relevant_points = {}; 25 std::vector<geo::utm::zonable_wgs84_point> relevant_points = {};
27 std::vector<std::shared_ptr<geo::linestring>> relevant_line_strings = {}; 26 std::vector<std::shared_ptr<geo::utm::zonable_wgs84_linestring>>
27 relevant_line_strings = {};
28}; 28};
29 29
30export struct situation 30export struct situation
31{ 31{
32 std::string id; 32 std::string id;
33 std::optional<geo::point> location = 33 std::optional<geo::utm::zonable_wgs84_point> location =
34 std::nullopt; // as shown on the map, not used for querying 34 std::nullopt; // as shown on the map, not used for querying
35 std::vector<std::string> comments = {}; 35 std::vector<std::string> comments = {};
36 std::vector<std::shared_ptr<road_closure>> road_closures = {}; 36 std::vector<std::shared_ptr<road_closure>> road_closures = {};
@@ -47,7 +47,7 @@ export class loader
47 // ETRS 89 (EPSG:4258) -> WGS 84 (EPSG:4326) 47 // ETRS 89 (EPSG:4258) -> WGS 84 (EPSG:4326)
48 bgeo::srs:: 48 bgeo::srs::
49 transformation<bgeo::srs::static_epsg<4258>, bgeo::srs::static_epsg<4326>> 49 transformation<bgeo::srs::static_epsg<4258>, bgeo::srs::static_epsg<4326>>
50 etrs89_to_wgs84_{}; 50 etrs89_to_wgs84_ = {};
51 51
52 std::multiset<std::string> warnings_; 52 std::multiset<std::string> warnings_;
53 53
diff --git a/server/src/geo.cpp b/server/src/geo.cpp
deleted file mode 100644
index 1776f3c..0000000
--- a/server/src/geo.cpp
+++ /dev/null
@@ -1,42 +0,0 @@
1module;
2
3#include <boost/geometry.hpp>
4
5module routemon:geo$impl;
6
7import :geo;
8
9namespace routemon::geo {
10
11auto split_linestring_with_overlap_segments(
12 linestring const& ls, double max_split_distance_m,
13 std::vector<linestring>& append_to) -> void
14{
15 if (bgeo::is_empty(ls))
16 return;
17
18 auto current_ls = linestring{};
19 auto current_ls_length = 0.0;
20 auto previous = std::optional<point>{};
21 bgeo::for_each_point(
22 ls,
23 [&](point p) -> void
24 {
25 bgeo::append(current_ls, p);
26 if (previous)
27 {
28 auto d = bgeo::distance(*previous, p, vincenty_strategy());
29 current_ls_length += d;
30 if (current_ls_length > max_split_distance_m)
31 {
32 append_to.push_back(std::move(current_ls));
33 current_ls = linestring{*previous, p};
34 current_ls_length = d;
35 }
36 }
37 previous = p;
38 });
39 append_to.emplace_back(std::move(current_ls));
40}
41
42} // namespace routemon::geo
diff --git a/server/src/geo.cppm b/server/src/geo.cppm
deleted file mode 100644
index b3513c8..0000000
--- a/server/src/geo.cppm
+++ /dev/null
@@ -1,23 +0,0 @@
1module;
2
3#include <boost/geometry.hpp>
4
5export module routemon:geo;
6
7export namespace bgeo = boost::geometry;
8
9export namespace routemon::geo {
10
11using point = bgeo::model::point<double, 2, bgeo::cs::geographic<bgeo::degree>>;
12using linestring = bgeo::model::linestring<point>;
13using box = bgeo::model::box<point>;
14using polygon = bgeo::model::polygon<geo::point>;
15using multi_polygon = bgeo::model::multi_polygon<polygon>;
16using stype = bgeo::srs::spheroid<double>;
17using vincenty_strategy = bgeo::strategy::distance::vincenty<stype>;
18
19auto split_linestring_with_overlap_segments(
20 linestring const& ls, double max_split_distance_m,
21 std::vector<linestring>& append_to) -> void;
22
23} // namespace routemon::geo
diff --git a/server/src/geo/geo.cppm b/server/src/geo/geo.cppm
new file mode 100644
index 0000000..446f1eb
--- /dev/null
+++ b/server/src/geo/geo.cppm
@@ -0,0 +1,7 @@
1module;
2
3#include <boost/geometry.hpp>
4
5export module routemon:geo;
6
7namespace bgeo = boost::geometry;
diff --git a/server/src/geo/multizonal.cppm b/server/src/geo/multizonal.cppm
new file mode 100644
index 0000000..f732994
--- /dev/null
+++ b/server/src/geo/multizonal.cppm
@@ -0,0 +1,240 @@
1module;
2
3#include <boost/geometry.hpp>
4#include <boost/geometry/srs/epsg.hpp>
5#include <boost/geometry/srs/transformation.hpp>
6
7export module routemon:geo.multizonal;
8
9import std;
10import :geo;
11import :geo.utm;
12import :geo.utm.zone_local;
13import :geo.wgs84;
14
15namespace routemon::geo::multizonal {
16
17template<std::default_initializable T>
18class multi_zone
19{
20 std::array<T, utm::zone::max().as_index() + 1> zones_;
21
22public:
23 auto operator[](utm::zone z) -> T&
24 {
25 return zones_[z.as_index()];
26 }
27
28 auto operator[](utm::zone z) const -> T const&
29 {
30 return zones_[z.as_index()];
31 }
32};
33
34namespace wgs84
35{
36
37// Transformations from WGS 84 (EPSG:4326)
38using transform_from_t = bgeo::srs::transformation<bgeo::srs::static_epsg<4326>>;
39
40class utm_transform : public transform_from_t
41{
42 utm::zone to_zone_;
43
44public:
45 explicit utm_transform(utm::zone to_zone)
46 : transform_from_t{{}, bgeo::srs::epsg{to_zone.wgs84_proj_epsg()}},
47 to_zone_{to_zone}
48 {
49 }
50
51 utm_transform()
52 : utm_transform{utm::zone::min()}
53 {}
54
55 auto apply(utm::zonable_wgs84_point p) -> utm::zone_local::point
56 {
57 auto local_p = utm::zone_local::point{to_zone_, 0.0, 0.0};
58 forward(p, local_p);
59 return local_p;
60 }
61};
62
63class utm_transforms : public multi_zone<utm_transform>
64{
65 utm_transforms()
66 {
67 for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
68 {
69 (*this)[z] = utm_transform{z};
70 }
71 }
72
73public:
74 static utm_transforms const& instance()
75 {
76 static utm_transforms inst;
77 return inst;
78 }
79};
80
81auto neighbor_utm_zone(utm::zonable_wgs84_point p) -> std::pair<utm::zone, double>
82{
83 auto separating_meridian_lon = std::round(p.lon() / 6.0) * 6.0;
84 auto closest_zone_middle = separating_meridian_lon < p.lon()
85 ? separating_meridian_lon - 3.0
86 : separating_meridian_lon + 3.0;
87 auto closest_zone = utm::zone::for_wgs84_point(utm::zonable_wgs84_point::from(geo::wgs84::normalized_point{geo::wgs84::from_lat_lon(p.lat(), closest_zone_middle)}).value());
88
89 auto separating_meridian_ls = geo::wgs84::linestring{
90 geo::wgs84::from_lat_lon(90.0, 0.0),
91 geo::wgs84::from_lat_lon(0.0, separating_meridian_lon),
92 geo::wgs84::from_lat_lon(-90.0, 0.0),
93 };
94 auto closest_zone_dist = bgeo::distance(p, separating_meridian_ls, geo::wgs84::vincenty_strategy{});
95
96 return std::make_pair(closest_zone, closest_zone_dist);
97}
98
99} // namespace wgs84
100
101struct multizone_linestring
102{
103 multi_zone<std::vector<utm::zone_local::prim::linestring>> segments;
104
105 // TODO: consider taking an input range instead
106 explicit multizone_linestring(utm::zonable_wgs84_linestring const& ls)
107 {
108 auto to_utm = wgs84::utm_transforms::instance();
109 auto working = multi_zone<utm::zone_local::prim::linestring>{};
110 auto mprev_p_zone = std::optional<utm::zone>{};
111 auto mprev_p_alt_zone = std::optional<utm::zone>{};
112
113 auto push = [&](utm::zonable_wgs84_point p, utm::zone z)
114 { working[z].push_back(to_utm[z].apply(p)); };
115 auto flush = [&](utm::zone z)
116 {
117 segments[z].push_back(std::move(working[z]));
118 working[z] = {};
119 };
120
121 for (auto const& p : ls)
122 {
123 auto p_zone = utm::zone::for_wgs84_point(p);
124 auto mp_alt_zone = std::optional<utm::zone>{};
125 if (auto [neighbor_zone, neighbor_zone_dist] = wgs84::neighbor_utm_zone(p);
126 neighbor_zone_dist < 30.0 /* m */)
127 mp_alt_zone = neighbor_zone;
128 assert(!mp_alt_zone || *mp_alt_zone != p_zone);
129
130 push(p, p_zone);
131 if (mp_alt_zone)
132 push(p, *mp_alt_zone);
133 if (mprev_p_zone && *mprev_p_zone != p_zone && mprev_p_zone != mp_alt_zone)
134 {
135 push(p, *mprev_p_zone);
136 flush(*mprev_p_zone);
137 }
138 if (mprev_p_alt_zone && *mprev_p_alt_zone != p_zone && mprev_p_alt_zone != mp_alt_zone)
139 {
140 push(p, *mprev_p_alt_zone);
141 flush(*mprev_p_alt_zone);
142 }
143
144 mprev_p_zone = p_zone;
145 mprev_p_alt_zone = mp_alt_zone;
146 }
147
148 if (mprev_p_zone && !working[*mprev_p_zone].empty())
149 {
150 flush(*mprev_p_zone);
151 }
152 if (mprev_p_alt_zone && !working[*mprev_p_alt_zone].empty())
153 {
154 flush(*mprev_p_alt_zone);
155 }
156 }
157};
158
159struct zoned_linestring_seg_seq
160{
161 std::vector<utm::zone_local::linestring> segments;
162
163 explicit zoned_linestring_seg_seq(utm::zonable_wgs84_linestring ls)
164 {
165 auto to_utm = wgs84::utm_transforms::instance();
166 auto mworking_seg = std::optional<utm::zone_local::linestring>{};
167
168 for (auto const& p : ls)
169 {
170 auto p_zone = utm::zone::for_wgs84_point(p);
171 if (mworking_seg && mworking_seg->zone != p_zone)
172 {
173 mworking_seg->push_back(to_utm[mworking_seg->zone].apply(p));
174 segments.push_back(std::move(*mworking_seg));
175 mworking_seg = std::nullopt;
176 }
177 if (!mworking_seg)
178 mworking_seg = utm::zone_local::linestring{p_zone};
179 mworking_seg->push_back(to_utm[p_zone].apply(p));
180 }
181
182 if (mworking_seg && !mworking_seg->empty())
183 {
184 segments.push_back(std::move(*mworking_seg));
185 }
186 }
187};
188
189template<class T>
190class linestring_rtree
191{
192public:
193 using index_value = std::tuple<utm::zone_local::prim::box, utm::zone_local::prim::linestring, T>;
194
195private:
196 multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>> local_rtrees_;
197
198public:
199 template<std::convertible_to<T> U>
200 auto insert(utm::zonable_wgs84_linestring const& ls, U&& arg) -> void
201 {
202 auto zone_segments = multizone_linestring{ls}.segments;
203 for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
204 {
205 for (auto ls : zone_segments[z])
206 {
207 utm::zone_local::prim::box blse =
208 bgeo::return_buffer<utm::zone_local::prim::box>(bgeo::return_envelope<utm::zone_local::prim::box>(ls), 30.0 /* m */);
209 local_rtrees_[z].insert(index_value{blse, std::move(ls), std::forward<U>(arg)});
210 }
211 }
212 }
213
214 auto size() const -> std::size_t
215 {
216 auto total_size = 0uz;
217 for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
218 total_size += local_rtrees_[z].size();
219 return total_size;
220 }
221
222 // Note: the same T may be generated more than once!
223 auto intersection(zoned_linestring_seg_seq const& lss) const -> std::generator<T const&>
224 {
225 for (auto const& ls : lss.segments)
226 {
227 // auto ls_box = bgeo::return_envelope<utm::zone_local::prim::box>(ls);
228 for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls));
229 it != local_rtrees_[ls.zone].qend(); it++)
230 {
231 if (bgeo::intersects(ls, std::get<1>(*it)))
232 {
233 co_yield std::get<2>(*it);
234 }
235 }
236 }
237 }
238};
239
240} // namespace routemon::geo::multizonal
diff --git a/server/src/geo/utm.cppm b/server/src/geo/utm.cppm
new file mode 100644
index 0000000..83da222
--- /dev/null
+++ b/server/src/geo/utm.cppm
@@ -0,0 +1,199 @@
1module;
2
3#include <cassert>
4#include <boost/geometry/geometries/geometries.hpp>
5
6export module routemon:geo.utm;
7
8import std;
9import :geo;
10import :geo.wgs84;
11
12namespace routemon::geo::utm
13{
14
15class zonable_wgs84_point
16{
17 wgs84::normalized_point p_;
18
19 explicit zonable_wgs84_point(wgs84::normalized_point p)
20 : p_{p}
21 {}
22
23public:
24 zonable_wgs84_point()
25 : p_{}
26 {
27 }
28
29 static auto from(wgs84::normalized_point p) -> std::optional<zonable_wgs84_point>
30 {
31 if (p.lat() < -80 || p.lat() > 84)
32 return std::nullopt;
33 return zonable_wgs84_point{p};
34 }
35
36 auto lat() const -> double
37 {
38 return p_.lat();
39 }
40
41 auto lon() const -> double
42 {
43 return p_.lon();
44 }
45
46 auto lon(double lon) -> void
47 {
48 p_.lon(lon);
49 }
50
51 auto lat(double lat) -> void
52 {
53 if (lat < -80 || lat > 84)
54 throw std::range_error{"latitude out of range for UTM"};
55 p_.lat(lat);
56 }
57};
58
59}
60
61namespace boost::geometry::traits {
62
63 namespace {
64
65 using zonable_wgs84_point = routemon::geo::utm::zonable_wgs84_point;
66
67 } // namespace <anonymous>
68
69 template<> struct tag<zonable_wgs84_point> { using type = point_tag; };
70 template<> struct dimension<zonable_wgs84_point> : boost::mpl::int_<2> {};
71 template<> struct coordinate_type<zonable_wgs84_point> { using type = double; };
72 template<> struct coordinate_system<zonable_wgs84_point> { using type = routemon::geo::wgs84::cs; };
73
74 template<std::size_t index>
75 struct access<zonable_wgs84_point, index> {
76 static inline auto get(zonable_wgs84_point const& p) -> double
77 {
78 if constexpr (index == 0)
79 return p.lon();
80 else if constexpr (index == 1)
81 return p.lat();
82 else static_assert(false, "Out of range");
83 }
84
85 static inline auto set(zonable_wgs84_point& p, double v) -> void
86 {
87 if constexpr (index == 0)
88 p.lon(v);
89 else if constexpr (index == 1)
90 p.lat(v);
91 else static_assert(false, "Out of range");
92 }
93 };
94
95} // namespace boost::geometry::traits
96
97namespace routemon::geo::utm
98{
99
100using zonable_wgs84_linestring = bgeo::model::linestring<zonable_wgs84_point, std::vector>;
101
102// In the sense of the common WGS84 subdivisions by simple northing
103// and easting (so no Norway/Svalbard exceptions).
104class zone
105{
106 // Note: calculations heavily depend on the values of the variants.
107 enum class hemisphere : std::uint8_t
108 {
109 northern = 0,
110 southern = 1,
111 };
112
113 // Negative if in the southern hemisphere
114 // Equal to (zone_no - 1) * 2 + hemisphere
115 // In range [0, 119]
116 std::uint8_t zone_;
117
118 static constexpr auto zone_no_valid(std::uint8_t zone_no) -> bool
119 {
120 return 0 < zone_no && zone_no <= 60;
121 }
122
123 static constexpr auto from(std::uint8_t zone_no, hemisphere h) -> std::uint8_t
124 {
125 if (!zone_no_valid(zone_no))
126 throw std::invalid_argument{"invalid UTM zone number"};
127 return (zone_no - 1) * 2 + static_cast<std::uint8_t>(h);
128 }
129
130public:
131 explicit constexpr zone(std::uint8_t zone_no, hemisphere h)
132 : zone_{from(zone_no, h)}
133 {}
134
135 auto hemisphere() const -> enum hemisphere
136 {
137 return static_cast<enum hemisphere>(zone_ % 2);
138 }
139
140 // Return value in range [1, 60]
141 auto zone_no() const -> std::uint8_t
142 {
143 return 1 + zone_ / 2;
144 }
145
146 // Return value in range [min(), max()]
147 constexpr auto as_index() const -> std::uint8_t
148 {
149 return zone_;
150 }
151
152 static auto for_wgs84_point(zonable_wgs84_point p) noexcept -> zone
153 {
154 auto zone_no = static_cast<std::uint8_t>(1 + (p.lon() + 180.0) / 6.0);
155 auto northern = p.lat() >= 0.0;
156 return zone{zone_no, northern ? hemisphere::northern : hemisphere::southern};
157 }
158
159 auto wgs84_proj_epsg() const -> int
160 {
161 switch (hemisphere())
162 {
163 case hemisphere::northern:
164 return 32600 + zone_no();
165 case hemisphere::southern:
166 return 32700 + zone_no();
167 }
168 }
169
170 // Guarantee: min().as_index() == 0.
171 static constexpr auto min() noexcept -> zone
172 {
173 return zone{1, hemisphere::northern};
174 }
175
176 static constexpr auto max() noexcept -> zone
177 {
178 return zone{60, hemisphere::southern};
179 }
180
181 auto next() const -> zone
182 {
183 assert(zone_ <= max().as_index());
184 if (zone_ == max().as_index())
185 throw std::range_error{"cannot take next of greatest UTM zone"};
186 auto copy = zone{*this};
187 copy.zone_++;
188 return copy;
189 }
190
191 auto operator==(zone rhs) const -> bool
192 {
193 return zone_ == rhs.zone_;
194 }
195};
196static_assert(zone::min().as_index() == 0);
197static_assert(zone::max().as_index() == 119);
198
199} // namespace routemon::geo::utm
diff --git a/server/src/geo/utm_zone_local.cppm b/server/src/geo/utm_zone_local.cppm
new file mode 100644
index 0000000..98a1344
--- /dev/null
+++ b/server/src/geo/utm_zone_local.cppm
@@ -0,0 +1,74 @@
1module;
2
3#include <boost/geometry.hpp>
4#include <boost/geometry/geometries/register/linestring.hpp>
5
6export module routemon:geo.utm.zone_local;
7
8import :geo;
9import :geo.utm;
10
11// Coordinates that are local to some known UTM zone.
12namespace routemon::geo::utm::zone_local {
13
14 using cs = bgeo::cs::cartesian;
15
16 namespace prim {
17
18 using point = bgeo::model::d2::point_xy<double>;
19 using linestring = bgeo::model::linestring<point>;
20 using box = bgeo::model::box<point>;
21
22 } // namespace prim
23
24 struct point : public prim::point
25 {
26 zone zone;
27
28 explicit point(class zone zone, double x, double y)
29 : prim::point{x, y}, zone{zone}
30 {
31 }
32 };
33
34 struct linestring : public prim::linestring
35 {
36 zone zone;
37
38 explicit linestring(class zone zone)
39 : zone{zone}
40 {
41 }
42 };
43
44} // namespace routemon::geo::utm::zone_local
45
46BOOST_GEOMETRY_REGISTER_LINESTRING(routemon::geo::utm::zone_local::linestring)
47
48namespace boost::geometry::traits {
49
50 namespace {
51
52 using zone_local_point = routemon::geo::utm::zone_local::point;
53
54 } // namespace <anonymous>
55
56 template<> struct tag<zone_local_point> { using type = point_tag; };
57 template<> struct dimension<zone_local_point> : boost::mpl::int_<2> {};
58 template<> struct coordinate_type<zone_local_point> { using type = double; };
59 template<> struct coordinate_system<zone_local_point> { using type = routemon::geo::utm::zone_local::cs; };
60
61 template<std::size_t index>
62 struct access<zone_local_point, index> {
63 static inline auto get(zone_local_point const& p) -> double
64 {
65 return bgeo::get<index>(static_cast<routemon::geo::utm::zone_local::prim::point const&>(p));
66 }
67
68 static inline auto set(zone_local_point& p, double v) -> void
69 {
70 return bgeo::set<index>(static_cast<routemon::geo::utm::zone_local::prim::point&>(p), v);
71 }
72 };
73
74} // namespace boost::geometry::traits
diff --git a/server/src/geo/wgs84.cppm b/server/src/geo/wgs84.cppm
new file mode 100644
index 0000000..d574d39
--- /dev/null
+++ b/server/src/geo/wgs84.cppm
@@ -0,0 +1,115 @@
1module;
2
3#include <boost/geometry.hpp>
4
5export module routemon:geo.wgs84;
6
7import :geo;
8
9namespace routemon::geo::wgs84 {
10
11using cs = bgeo::cs::geographic<bgeo::degree>;
12using point = bgeo::model::point<double, 2, cs>;
13// TODO: guarantee that linestring provides non-static member
14// auto reserve(std::size_t) -> void
15// Perhaps better yet: guarantee that the backing container is a std::vector.
16using linestring = bgeo::model::linestring<point, std::vector /* the default */>;
17using box = bgeo::model::box<point>;
18using stype = bgeo::srs::spheroid<double>;
19using vincenty_strategy = bgeo::strategy::distance::vincenty<stype>;
20
21auto from_lat_lon(double lat, double lon) -> point
22{
23 return point{lon, lat};
24}
25
26auto lat(point p) -> double
27{
28 return bgeo::get<1>(p);
29}
30
31auto lon(point p) -> double
32{
33 return bgeo::get<0>(p);
34}
35
36// Point p with latitude in [-90, 90] and longitude in [-180, 180)
37auto normalize(point p) -> point
38{
39 // Equivalent degrees in the range [0, 360)
40 auto nonneg_mod360 = [](double t) -> double
41 { return std::remainder(std::remainder(t, 360.0) + 360.0, 360.0); };
42
43 auto p_lon = lon(p);
44 auto p_lat = nonneg_mod360(lat(p)) - 90.0;
45 if (90.0 < p_lat)
46 {
47 assert(p_lat < 270.0); // by nonneg_mod360
48 p_lon += 180.0;
49 p_lat -= 180.0;
50 }
51 p_lon = nonneg_mod360(p_lon + 180.0) - 180.0;
52
53 return from_lat_lon(p_lat, p_lon);
54}
55
56auto is_normalized(point p) -> bool
57{
58 return -90 <= lat(p) && lat(p) <= 90 && -180 <= lon(p) && lon(p) < 180;
59}
60
61class normalized_point
62{
63 point p_;
64
65public:
66 normalized_point(point p)
67 : p_{is_normalized(p) ? p : normalize(p)}
68 {}
69
70 normalized_point()
71 : p_{}
72 {}
73
74 auto lat() const -> double
75 {
76 return bgeo::get<1>(p_);
77 }
78
79 auto lon() const -> double
80 {
81 return bgeo::get<0>(p_);
82 }
83
84 auto lat(double lat) -> void
85 {
86 if (lat < -90 || lat > 90)
87 throw std::range_error{"latitude out of range"};
88 bgeo::set<1>(p_, lat);
89 }
90
91 auto lon(double lon) -> void
92 {
93 if (lon < -180 || lon > 180)
94 throw std::range_error{"longitude out of range"};
95 bgeo::set<0>(p_, lon);
96 }
97};
98
99class normalized_linestring
100{
101 std::vector<normalized_point> ls_;
102
103public:
104 auto empty() const -> bool
105 {
106 return ls_.empty();
107 }
108
109 auto push_back(normalized_point p) -> void
110 {
111 ls_.push_back(p);
112 }
113};
114
115} // namespace routemon::geo::wgs84
diff --git a/server/src/gpx.cpp b/server/src/gpx.cpp
index 078248f..2193e6a 100644
--- a/server/src/gpx.cpp
+++ b/server/src/gpx.cpp
@@ -22,7 +22,7 @@ auto qname(std::string_view local) -> xml::qname_view
22} 22}
23 23
24auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs) 24auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs)
25 -> xml::parser<geo::point> 25 -> xml::parser<geo::utm::zonable_wgs84_point>
26{ 26{
27 auto parse_xml_double = [](std::string_view sv) -> std::optional<double> 27 auto parse_xml_double = [](std::string_view sv) -> std::optional<double>
28 { return util::parse_double(sv, std::chars_format::fixed); }; 28 { return util::parse_double(sv, std::chars_format::fixed); };
@@ -43,8 +43,12 @@ auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs)
43 throw std::runtime_error{ 43 throw std::runtime_error{
44 "expected valid latitude and longitude for waypoint" 44 "expected valid latitude and longitude for waypoint"
45 }; 45 };
46 auto mp = geo::utm::zonable_wgs84_point::from(
47 geo::wgs84::normalized_point{geo::wgs84::from_lat_lon(*mlat, *mlon)});
48 if (!mp)
49 throw std::runtime_error{"expected waypoint to be within UTM range"};
46 co_await xml::ignore_contents(e); 50 co_await xml::ignore_contents(e);
47 co_return geo::point{*mlon, *mlat}; 51 co_return *mp;
48} 52}
49 53
50auto parse_trkseg(xml::executor_ref e, xml::attribute_view) 54auto parse_trkseg(xml::executor_ref e, xml::attribute_view)
@@ -121,7 +125,7 @@ auto parse_metadata(xml::executor_ref e, xml::attribute_view)
121} 125}
122 126
123auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs) 127auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs)
124 -> xml::parser<geo::point> 128 -> xml::parser<geo::utm::zonable_wgs84_point>
125{ 129{
126 auto parse_xml_double = [](std::string_view sv) -> std::optional<double> 130 auto parse_xml_double = [](std::string_view sv) -> std::optional<double>
127 { return util::parse_double(sv, std::chars_format::fixed); }; 131 { return util::parse_double(sv, std::chars_format::fixed); };
@@ -142,8 +146,12 @@ auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs)
142 throw std::runtime_error{ 146 throw std::runtime_error{
143 "expected valid latitude and longitude for waypoint" 147 "expected valid latitude and longitude for waypoint"
144 }; 148 };
149 auto mp = geo::utm::zonable_wgs84_point::from(
150 geo::wgs84::normalized_point{geo::wgs84::from_lat_lon(*mlat, *mlon)});
151 if (!mp)
152 throw std::runtime_error{"expected waypoint to be within UTM range"};
145 co_await xml::ignore_contents(e); 153 co_await xml::ignore_contents(e);
146 co_return geo::point{*mlon, *mlat}; 154 co_return *mp;
147} 155}
148 156
149auto parse_trkseg(xml::executor_ref e, xml::attribute_view) 157auto parse_trkseg(xml::executor_ref e, xml::attribute_view)
diff --git a/server/src/gpx.cppm b/server/src/gpx.cppm
index 5872805..6267443 100644
--- a/server/src/gpx.cppm
+++ b/server/src/gpx.cppm
@@ -1,7 +1,7 @@
1export module routemon:gpx; 1export module routemon:gpx;
2 2
3import std; 3import std;
4import :geo; 4import :geo.utm;
5import :util; 5import :util;
6import :xml; 6import :xml;
7 7
@@ -15,7 +15,7 @@ struct metadata
15 15
16struct track_segment 16struct track_segment
17{ 17{
18 geo::linestring waypoints; 18 geo::utm::zonable_wgs84_linestring waypoints;
19}; 19};
20 20
21struct track 21struct track
diff --git a/server/src/log.cpp b/server/src/log.cpp
index 7fe725f..bb9fde2 100644
--- a/server/src/log.cpp
+++ b/server/src/log.cpp
@@ -39,7 +39,7 @@ auto sink::write(tmp_message msg) -> void
39 sos << "] " << msg.txt; 39 sos << "] " << msg.txt;
40 for (auto const& [k, v] : msg.attrs) 40 for (auto const& [k, v] : msg.attrs)
41 sos << " " << k << "=" << std::quoted(v); 41 sos << " " << k << "=" << std::quoted(v);
42 sos << '\n'; 42 sos << std::endl; // we want the flush!
43} 43}
44 44
45auto make_sink(level lvl) -> std::shared_ptr<sink> 45auto make_sink(level lvl) -> std::shared_ptr<sink>