summaryrefslogtreecommitdiffstats
path: root/server/src
diff options
context:
space:
mode:
Diffstat (limited to 'server/src')
-rw-r--r--server/src/api.cpp3
-rw-r--r--server/src/api.cppm6
-rw-r--r--server/src/geo/multizonal.cppm307
-rw-r--r--server/src/geo/utm.cppm140
-rw-r--r--server/src/geo/utm_multizonal.cppm343
-rw-r--r--server/src/geo/utm_zone_local.cppm95
-rw-r--r--server/src/geo/wgs84.cppm51
7 files changed, 490 insertions, 455 deletions
diff --git a/server/src/api.cpp b/server/src/api.cpp
index e35e1b3..41c49c8 100644
--- a/server/src/api.cpp
+++ b/server/src/api.cpp
@@ -172,7 +172,8 @@ auto handler::process_gpx(gpx::file&& gpx_file)
172 { 172 {
173 l_.debug("Checking part [{}/{}]", ++i, segments.size()); 173 l_.debug("Checking part [{}/{}]", ++i, segments.size());
174 174
175 auto part_zone_lss = geo::multizonal::zoned_linestring_seg_seq{part}; 175 auto part_zone_lss =
176 geo::utm::multizonal::split_linestring_across_zones(part);
176 // TODO: consider buffering with min_distance_ 177 // TODO: consider buffering with min_distance_
177 auto part_box = geo::wgs84::box{}; 178 auto part_box = geo::wgs84::box{};
178 bgeo::envelope(part, part_box); 179 bgeo::envelope(part, part_box);
diff --git a/server/src/api.cppm b/server/src/api.cppm
index 0bb0868..7166ad0 100644
--- a/server/src/api.cppm
+++ b/server/src/api.cppm
@@ -8,7 +8,7 @@ export module routemon:api;
8import std; 8import std;
9import :datex2; 9import :datex2;
10import :geo.wgs84; 10import :geo.wgs84;
11import :geo.multizonal; 11import :geo.utm.multizonal;
12import :gpx; 12import :gpx;
13import :log; 13import :log;
14import :time; 14import :time;
@@ -81,9 +81,9 @@ class handler
81 constexpr static auto const min_distance_ = 5.0; 81 constexpr static auto const min_distance_ = 5.0;
82 82
83 using blse_index_value = std::shared_ptr<datex2::road_closure>; 83 using blse_index_value = std::shared_ptr<datex2::road_closure>;
84 using blse_index = geo::multizonal::linestring_rtree<blse_index_value>; 84 using blse_index = geo::utm::multizonal::linestring_rtree<blse_index_value>;
85 using bpe_index_value = std::shared_ptr<datex2::road_closure>; 85 using bpe_index_value = std::shared_ptr<datex2::road_closure>;
86 using bpe_index = geo::multizonal::point_rtree<bpe_index_value>; 86 using bpe_index = geo::utm::multizonal::point_rtree<bpe_index_value>;
87 87
88 log::logger l_; 88 log::logger l_;
89 datex2::situation_publication pub_; 89 datex2::situation_publication pub_;
diff --git a/server/src/geo/multizonal.cppm b/server/src/geo/multizonal.cppm
deleted file mode 100644
index 7930367..0000000
--- a/server/src/geo/multizonal.cppm
+++ /dev/null
@@ -1,307 +0,0 @@
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 auto 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 using multi_linestring = bgeo::model::multi_linestring<utm::zone_local::prim::linestring>;
232
233 auto relevant_ls_parts = multi_linestring{};
234 bgeo::intersection(std::get<0>(*it), static_cast<utm::zone_local::prim::linestring const&>(ls), relevant_ls_parts);
235 if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */)
236 {
237 co_yield std::get<2>(*it);
238 }
239 }
240 }
241 }
242};
243
244template<class T>
245class point_rtree
246{
247public:
248 using index_value = std::tuple<utm::zone_local::prim::box, utm::zone_local::prim::point, T>;
249
250private:
251 multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>> local_rtrees_;
252
253public:
254 template<std::convertible_to<T> U>
255 auto insert(utm::zonable_wgs84_point const& p, U&& arg) -> void
256 {
257 auto to_utm = wgs84::utm_transforms::instance();
258
259 auto p_zone = utm::zone::for_wgs84_point(p);
260 auto mp_alt_zone = std::optional<utm::zone>{};
261 if (auto [neighbor_zone, neighbor_zone_dist] = wgs84::neighbor_utm_zone(p);
262 neighbor_zone_dist < 30.0 /* m */)
263 mp_alt_zone = neighbor_zone;
264
265 auto insert = [&](utm::zone z) -> void
266 {
267 auto p_utm = to_utm[z].apply(p);
268 auto bpe = bgeo::return_buffer<utm::zone_local::prim::box>(bgeo::return_envelope<utm::zone_local::prim::box>(p_utm), 30.0 /* m */);
269 local_rtrees_[z].insert(index_value{bpe, p_utm, arg});
270 };
271
272 insert(p_zone);
273 if (mp_alt_zone)
274 insert(*mp_alt_zone);
275 }
276
277 auto size() const -> std::size_t
278 {
279 auto total_size = 0uz;
280 for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
281 total_size += local_rtrees_[z].size();
282 return total_size;
283 }
284
285 // Note: the same T may be generated more than once!
286 auto intersection(zoned_linestring_seg_seq const& lss) const -> std::generator<T const&>
287 {
288 for (auto const& ls : lss.segments)
289 {
290 // auto ls_box = bgeo::return_envelope<utm::zone_local::prim::box>(ls);
291 for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls));
292 it != local_rtrees_[ls.zone].qend(); it++)
293 {
294 using multi_linestring = bgeo::model::multi_linestring<utm::zone_local::prim::linestring>;
295
296 auto relevant_ls_parts = multi_linestring{};
297 bgeo::intersection(std::get<0>(*it), static_cast<utm::zone_local::prim::linestring const&>(ls), relevant_ls_parts);
298 if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */)
299 {
300 co_yield std::get<2>(*it);
301 }
302 }
303 }
304 }
305};
306
307} // namespace routemon::geo::multizonal
diff --git a/server/src/geo/utm.cppm b/server/src/geo/utm.cppm
index 83da222..c58ffb4 100644
--- a/server/src/geo/utm.cppm
+++ b/server/src/geo/utm.cppm
@@ -1,7 +1,7 @@
1module; 1module;
2 2
3#include <cassert>
4#include <boost/geometry/geometries/geometries.hpp> 3#include <boost/geometry/geometries/geometries.hpp>
4#include <cassert>
5 5
6export module routemon:geo.utm; 6export module routemon:geo.utm;
7 7
@@ -9,46 +9,32 @@ import std;
9import :geo; 9import :geo;
10import :geo.wgs84; 10import :geo.wgs84;
11 11
12namespace routemon::geo::utm 12namespace routemon::geo::utm {
13{
14 13
15class zonable_wgs84_point 14class zonable_wgs84_point
16{ 15{
17 wgs84::normalized_point p_; 16 wgs84::normalized_point p_;
18 17
19 explicit zonable_wgs84_point(wgs84::normalized_point p) 18 explicit inline zonable_wgs84_point(wgs84::normalized_point p) : p_{p} {}
20 : p_{p}
21 {}
22 19
23public: 20public:
24 zonable_wgs84_point() 21 inline zonable_wgs84_point() : p_{} {}
25 : p_{}
26 {
27 }
28 22
29 static auto from(wgs84::normalized_point p) -> std::optional<zonable_wgs84_point> 23 inline static auto from(wgs84::normalized_point p)
24 -> std::optional<zonable_wgs84_point>
30 { 25 {
31 if (p.lat() < -80 || p.lat() > 84) 26 if (p.lat() < -80 || p.lat() > 84)
32 return std::nullopt; 27 return std::nullopt;
33 return zonable_wgs84_point{p}; 28 return zonable_wgs84_point{p};
34 } 29 }
35 30
36 auto lat() const -> double 31 inline auto lat() const -> double { return p_.lat(); }
37 {
38 return p_.lat();
39 }
40 32
41 auto lon() const -> double 33 inline auto lon() const -> double { return p_.lon(); }
42 {
43 return p_.lon();
44 }
45 34
46 auto lon(double lon) -> void 35 inline auto lon(double lon) -> void { p_.lon(lon); }
47 {
48 p_.lon(lon);
49 }
50 36
51 auto lat(double lat) -> void 37 inline auto lat(double lat) -> void
52 { 38 {
53 if (lat < -80 || lat > 84) 39 if (lat < -80 || lat > 84)
54 throw std::range_error{"latitude out of range for UTM"}; 40 throw std::range_error{"latitude out of range for UTM"};
@@ -56,48 +42,69 @@ public:
56 } 42 }
57}; 43};
58 44
59} 45} // namespace routemon::geo::utm
60 46
61namespace boost::geometry::traits { 47namespace boost::geometry::traits {
62 48
63 namespace { 49namespace {
64 50
65 using zonable_wgs84_point = routemon::geo::utm::zonable_wgs84_point; 51using zonable_wgs84_point = routemon::geo::utm::zonable_wgs84_point;
66 52
67 } // namespace <anonymous> 53} // namespace
68 54
69 template<> struct tag<zonable_wgs84_point> { using type = point_tag; }; 55template <>
70 template<> struct dimension<zonable_wgs84_point> : boost::mpl::int_<2> {}; 56struct tag<zonable_wgs84_point>
71 template<> struct coordinate_type<zonable_wgs84_point> { using type = double; }; 57{
72 template<> struct coordinate_system<zonable_wgs84_point> { using type = routemon::geo::wgs84::cs; }; 58 using type = point_tag;
59};
73 60
74 template<std::size_t index> 61template <>
75 struct access<zonable_wgs84_point, index> { 62struct dimension<zonable_wgs84_point> : boost::mpl::int_<2>
76 static inline auto get(zonable_wgs84_point const& p) -> double 63{
77 { 64};
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 65
85 static inline auto set(zonable_wgs84_point& p, double v) -> void 66template <>
86 { 67struct coordinate_type<zonable_wgs84_point>
87 if constexpr (index == 0) 68{
88 p.lon(v); 69 using type = double;
89 else if constexpr (index == 1) 70};
90 p.lat(v);
91 else static_assert(false, "Out of range");
92 }
93 };
94 71
95} // namespace boost::geometry::traits 72template <>
73struct coordinate_system<zonable_wgs84_point>
74{
75 using type = routemon::geo::wgs84::cs;
76};
96 77
97namespace routemon::geo::utm 78template <std::size_t index>
79struct access<zonable_wgs84_point, index>
98{ 80{
81 static inline auto get(zonable_wgs84_point const& p) -> double
82 {
83 if constexpr (index == 0)
84 return p.lon();
85 else if constexpr (index == 1)
86 return p.lat();
87 else
88 static_assert(false, "Out of range");
89 }
90
91 static inline auto set(zonable_wgs84_point& p, double v) -> void
92 {
93 if constexpr (index == 0)
94 p.lon(v);
95 else if constexpr (index == 1)
96 p.lat(v);
97 else
98 static_assert(false, "Out of range");
99 }
100};
101
102} // namespace boost::geometry::traits
99 103
100using zonable_wgs84_linestring = bgeo::model::linestring<zonable_wgs84_point, std::vector>; 104namespace routemon::geo::utm {
105
106using zonable_wgs84_linestring =
107 bgeo::model::linestring<zonable_wgs84_point, std::vector>;
101 108
102// In the sense of the common WGS84 subdivisions by simple northing 109// In the sense of the common WGS84 subdivisions by simple northing
103// and easting (so no Norway/Svalbard exceptions). 110// and easting (so no Norway/Svalbard exceptions).
@@ -130,33 +137,30 @@ class zone
130public: 137public:
131 explicit constexpr zone(std::uint8_t zone_no, hemisphere h) 138 explicit constexpr zone(std::uint8_t zone_no, hemisphere h)
132 : zone_{from(zone_no, h)} 139 : zone_{from(zone_no, h)}
133 {} 140 {
141 }
134 142
135 auto hemisphere() const -> enum hemisphere 143 constexpr auto hemisphere() const -> enum hemisphere
136 { 144 {
137 return static_cast<enum hemisphere>(zone_ % 2); 145 return static_cast<enum hemisphere>(zone_ % 2);
138 } 146 }
139 147
140 // Return value in range [1, 60] 148 // Return value in range [1, 60]
141 auto zone_no() const -> std::uint8_t 149 constexpr auto zone_no() const -> std::uint8_t { return 1 + zone_ / 2; }
142 {
143 return 1 + zone_ / 2;
144 }
145 150
146 // Return value in range [min(), max()] 151 // Return value in range [min(), max()]
147 constexpr auto as_index() const -> std::uint8_t 152 constexpr auto as_index() const -> std::uint8_t { return zone_; }
148 {
149 return zone_;
150 }
151 153
152 static auto for_wgs84_point(zonable_wgs84_point p) noexcept -> zone 154 static auto for_wgs84_point(zonable_wgs84_point p) noexcept -> zone
153 { 155 {
154 auto zone_no = static_cast<std::uint8_t>(1 + (p.lon() + 180.0) / 6.0); 156 auto zone_no = static_cast<std::uint8_t>(1 + (p.lon() + 180.0) / 6.0);
155 auto northern = p.lat() >= 0.0; 157 auto northern = p.lat() >= 0.0;
156 return zone{zone_no, northern ? hemisphere::northern : hemisphere::southern}; 158 return zone{
159 zone_no, northern ? hemisphere::northern : hemisphere::southern
160 };
157 } 161 }
158 162
159 auto wgs84_proj_epsg() const -> int 163 constexpr auto wgs84_proj_epsg() const -> int
160 { 164 {
161 switch (hemisphere()) 165 switch (hemisphere())
162 { 166 {
@@ -178,7 +182,7 @@ public:
178 return zone{60, hemisphere::southern}; 182 return zone{60, hemisphere::southern};
179 } 183 }
180 184
181 auto next() const -> zone 185 constexpr auto next() const -> zone
182 { 186 {
183 assert(zone_ <= max().as_index()); 187 assert(zone_ <= max().as_index());
184 if (zone_ == max().as_index()) 188 if (zone_ == max().as_index())
@@ -188,7 +192,7 @@ public:
188 return copy; 192 return copy;
189 } 193 }
190 194
191 auto operator==(zone rhs) const -> bool 195 constexpr auto operator==(zone rhs) const -> bool
192 { 196 {
193 return zone_ == rhs.zone_; 197 return zone_ == rhs.zone_;
194 } 198 }
diff --git a/server/src/geo/utm_multizonal.cppm b/server/src/geo/utm_multizonal.cppm
new file mode 100644
index 0000000..540e309
--- /dev/null
+++ b/server/src/geo/utm_multizonal.cppm
@@ -0,0 +1,343 @@
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.utm.multizonal;
8
9import std;
10import :geo;
11import :geo.utm;
12import :geo.utm.zone_local;
13import :geo.wgs84;
14
15namespace routemon::geo::utm::multizonal {
16
17template <std::default_initializable T>
18class multi_zone
19{
20 std::array<T, zone::max().as_index() + 1> zones_;
21
22public:
23 auto operator[](zone z) -> T& { return zones_[z.as_index()]; }
24
25 auto operator[](zone z) const -> T const& { return zones_[z.as_index()]; }
26};
27
28// Transformations from WGS 84 (EPSG:4326)
29using from_wgs84_transform_base =
30 bgeo::srs::transformation<bgeo::srs::static_epsg<4326>>;
31
32class from_wgs84_transform : public from_wgs84_transform_base
33{
34 zone to_zone_;
35
36public:
37 explicit from_wgs84_transform(zone to_zone)
38 : from_wgs84_transform_base{{}, bgeo::srs::epsg{to_zone.wgs84_proj_epsg()}},
39 to_zone_{to_zone}
40 {
41 }
42
43 from_wgs84_transform() : from_wgs84_transform{zone::min()} {}
44
45 auto apply(zonable_wgs84_point p) -> zone_local::point
46 {
47 auto local_p = zone_local::point{to_zone_, 0.0, 0.0};
48 forward(p, local_p);
49 return local_p;
50 }
51};
52
53class from_wgs84_transforms : public multi_zone<from_wgs84_transform>
54{
55 from_wgs84_transforms()
56 {
57 for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
58 {
59 (*this)[z] = from_wgs84_transform{z};
60 }
61 }
62
63public:
64 static from_wgs84_transforms const& instance()
65 {
66 static from_wgs84_transforms inst;
67 return inst;
68 }
69};
70
71auto neighbor_utm_zone(utm::zonable_wgs84_point p)
72 -> std::pair<utm::zone, double>
73{
74 auto separating_meridian_lon = std::round(p.lon() / 6.0) * 6.0;
75 auto closest_zone_middle = separating_meridian_lon < p.lon()
76 ? separating_meridian_lon - 3.0
77 : separating_meridian_lon + 3.0;
78 auto closest_zone = utm::zone::for_wgs84_point(
79 utm::zonable_wgs84_point::from(
80 geo::wgs84::normalized_point{geo::wgs84::from_lat_lon(
81 p.lat(), closest_zone_middle)})
82 .value());
83
84 auto separating_meridian_ls = geo::wgs84::linestring{
85 geo::wgs84::from_lat_lon(90.0, 0.0),
86 geo::wgs84::from_lat_lon(0.0, separating_meridian_lon),
87 geo::wgs84::from_lat_lon(-90.0, 0.0),
88 };
89 auto closest_zone_dist = bgeo::distance(
90 p, separating_meridian_ls, geo::wgs84::vincenty_strategy{});
91
92 return std::make_pair(closest_zone, closest_zone_dist);
93}
94
95// Distribute a WGS 84 linestring over the UTM zones that it finds
96// itself in or near to.
97//
98// Since it is not at all unthinkable that some road works will end up
99// crossing different UTM zones, we need some machinery to work with
100// such situations. Since we want to report all situations within a
101// specified distance of the planned route, we must also consider the
102// situation where this is indeed the case, but the situation is
103// situated at the opposite side of a delineating UTM zone meridian
104// w.r.t. the planned route.
105//
106// To solve this issue, we create segments of the input line string
107// for multiple UTM zones when crossing zone boundaries or when very
108// close to zone boundaries.
109auto distribute_linestring_over_zones(utm::zonable_wgs84_linestring const& ls)
110 -> multi_zone<std::vector<utm::zone_local::prim::linestring>>
111{
112 auto result = multi_zone<std::vector<utm::zone_local::prim::linestring>>{};
113
114 auto to_utm = from_wgs84_transforms::instance();
115 auto working = multi_zone<utm::zone_local::prim::linestring>{};
116 auto mprev_p_zone = std::optional<utm::zone>{};
117 auto mprev_p_alt_zone = std::optional<utm::zone>{};
118
119 auto push = [&](utm::zonable_wgs84_point p, utm::zone z)
120 { working[z].push_back(to_utm[z].apply(p)); };
121 auto flush = [&](utm::zone z)
122 {
123 result[z].push_back(std::move(working[z]));
124 working[z] = {};
125 };
126
127 for (auto const& p : ls)
128 {
129 auto p_zone = utm::zone::for_wgs84_point(p);
130 auto mp_alt_zone = std::optional<utm::zone>{};
131 if (auto [neighbor_zone, neighbor_zone_dist] = neighbor_utm_zone(p);
132 neighbor_zone_dist < 30.0 /* m */)
133 mp_alt_zone = neighbor_zone;
134 assert(!mp_alt_zone || *mp_alt_zone != p_zone);
135
136 push(p, p_zone);
137 if (mp_alt_zone)
138 push(p, *mp_alt_zone);
139 if (mprev_p_zone && *mprev_p_zone != p_zone && mprev_p_zone != mp_alt_zone)
140 {
141 push(p, *mprev_p_zone);
142 flush(*mprev_p_zone);
143 }
144 if (mprev_p_alt_zone && *mprev_p_alt_zone != p_zone
145 && mprev_p_alt_zone != mp_alt_zone)
146 {
147 push(p, *mprev_p_alt_zone);
148 flush(*mprev_p_alt_zone);
149 }
150
151 mprev_p_zone = p_zone;
152 mprev_p_alt_zone = mp_alt_zone;
153 }
154
155 if (mprev_p_zone && !working[*mprev_p_zone].empty())
156 {
157 flush(*mprev_p_zone);
158 }
159 if (mprev_p_alt_zone && !working[*mprev_p_alt_zone].empty())
160 {
161 flush(*mprev_p_alt_zone);
162 }
163
164 return result;
165}
166
167// Split a WGS 84 line string into a sequence of UTM-zone-local line
168// strings, starting a new UTM-zone-local line string when the input
169// line string crosses a UTM zone boundary. The line segment that
170// crosses the zone boundary can be found in the UTM-zone-local line
171// strings for both zones which its points are in (we assume that line
172// segments will be short enough to not cause significant distortion
173// here, so we do not put in the effort to clip at the
174// zone-delineating meridian here).
175auto split_linestring_across_zones(utm::zonable_wgs84_linestring ls)
176 -> std::vector<utm::zone_local::linestring>
177{
178 auto splits = std::vector<utm::zone_local::linestring>{};
179
180 auto to_utm = from_wgs84_transforms::instance();
181 auto mworking_seg = std::optional<utm::zone_local::linestring>{};
182
183 for (auto const& p : ls)
184 {
185 auto p_zone = utm::zone::for_wgs84_point(p);
186 if (mworking_seg && mworking_seg->zone != p_zone)
187 {
188 mworking_seg->push_back(to_utm[mworking_seg->zone].apply(p));
189 splits.push_back(std::move(*mworking_seg));
190 mworking_seg = std::nullopt;
191 }
192 if (!mworking_seg)
193 mworking_seg = utm::zone_local::linestring{p_zone};
194 mworking_seg->push_back(to_utm[p_zone].apply(p));
195 }
196
197 if (mworking_seg && !mworking_seg->empty())
198 {
199 splits.push_back(std::move(*mworking_seg));
200 }
201
202 return splits;
203}
204
205template <class T>
206class linestring_rtree
207{
208public:
209 using index_value = std::
210 tuple<utm::zone_local::prim::box, utm::zone_local::prim::linestring, T>;
211
212private:
213 multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>>
214 local_rtrees_;
215
216public:
217 template <std::convertible_to<T> U>
218 auto insert(utm::zonable_wgs84_linestring const& ls, U&& arg) -> void
219 {
220 auto zone_segments = distribute_linestring_over_zones(ls);
221 for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
222 {
223 for (auto ls : zone_segments[z])
224 {
225 auto blse = bgeo::return_buffer<utm::zone_local::prim::box>(
226 bgeo::return_envelope<utm::zone_local::prim::box>(ls),
227 30.0 /* m */);
228 local_rtrees_[z].insert(
229 index_value{blse, std::move(ls), std::forward<U>(arg)});
230 }
231 }
232 }
233
234 auto size() const -> std::size_t
235 {
236 auto total_size = 0uz;
237 for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
238 total_size += local_rtrees_[z].size();
239 return total_size;
240 }
241
242 // Note: the same T may be generated more than once!
243 auto intersection(std::vector<utm::zone_local::linestring> const& lss) const
244 -> std::generator<T const&>
245 {
246 for (auto const& ls : lss)
247 {
248 for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls));
249 it != local_rtrees_[ls.zone].qend(); it++)
250 {
251 using multi_linestring =
252 bgeo::model::multi_linestring<utm::zone_local::prim::linestring>;
253
254 auto relevant_ls_parts = multi_linestring{};
255 bgeo::intersection(
256 std::get<0>(*it),
257 static_cast<utm::zone_local::prim::linestring const&>(ls),
258 relevant_ls_parts);
259 // it is somewhat awkward and arbitrary that this measurement is
260 // performed here
261 if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */)
262 {
263 co_yield std::get<2>(*it);
264 }
265 }
266 }
267 }
268};
269
270template <class T>
271class point_rtree
272{
273public:
274 using index_value =
275 std::tuple<utm::zone_local::prim::box, utm::zone_local::prim::point, T>;
276
277private:
278 multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>>
279 local_rtrees_;
280
281public:
282 template <std::convertible_to<T> U>
283 auto insert(utm::zonable_wgs84_point const& p, U&& arg) -> void
284 {
285 auto to_utm = from_wgs84_transforms::instance();
286
287 auto p_zone = utm::zone::for_wgs84_point(p);
288 auto mp_alt_zone = std::optional<utm::zone>{};
289 if (auto [neighbor_zone, neighbor_zone_dist] = neighbor_utm_zone(p);
290 neighbor_zone_dist < 30.0 /* m */)
291 mp_alt_zone = neighbor_zone;
292
293 auto insert = [&](utm::zone z) -> void
294 {
295 auto p_utm = to_utm[z].apply(p);
296 auto bpe = bgeo::return_buffer<utm::zone_local::prim::box>(
297 bgeo::return_envelope<utm::zone_local::prim::box>(p_utm),
298 30.0 /* m */);
299 local_rtrees_[z].insert(index_value{bpe, p_utm, arg});
300 };
301
302 insert(p_zone);
303 if (mp_alt_zone)
304 insert(*mp_alt_zone);
305 }
306
307 auto size() const -> std::size_t
308 {
309 auto total_size = 0uz;
310 for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
311 total_size += local_rtrees_[z].size();
312 return total_size;
313 }
314
315 // Note: the same T may be generated more than once!
316 auto intersection(std::vector<utm::zone_local::linestring> const& lss) const
317 -> std::generator<T const&>
318 {
319 for (auto const& ls : lss)
320 {
321 for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls));
322 it != local_rtrees_[ls.zone].qend(); it++)
323 {
324 using multi_linestring =
325 bgeo::model::multi_linestring<utm::zone_local::prim::linestring>;
326
327 auto relevant_ls_parts = multi_linestring{};
328 bgeo::intersection(
329 std::get<0>(*it),
330 static_cast<utm::zone_local::prim::linestring const&>(ls),
331 relevant_ls_parts);
332 // it is somewhat awkward and arbitrary that this measurement is
333 // performed here
334 if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */)
335 {
336 co_yield std::get<2>(*it);
337 }
338 }
339 }
340 }
341};
342
343} // namespace routemon::geo::utm::multizonal
diff --git a/server/src/geo/utm_zone_local.cppm b/server/src/geo/utm_zone_local.cppm
index 98a1344..b2147ad 100644
--- a/server/src/geo/utm_zone_local.cppm
+++ b/server/src/geo/utm_zone_local.cppm
@@ -11,35 +11,32 @@ import :geo.utm;
11// Coordinates that are local to some known UTM zone. 11// Coordinates that are local to some known UTM zone.
12namespace routemon::geo::utm::zone_local { 12namespace routemon::geo::utm::zone_local {
13 13
14 using cs = bgeo::cs::cartesian; 14using cs = bgeo::cs::cartesian;
15 15
16 namespace prim { 16namespace prim {
17 17
18 using point = bgeo::model::d2::point_xy<double>; 18using point = bgeo::model::d2::point_xy<double>;
19 using linestring = bgeo::model::linestring<point>; 19using linestring = bgeo::model::linestring<point>;
20 using box = bgeo::model::box<point>; 20using box = bgeo::model::box<point>;
21 21
22 } // namespace prim 22} // namespace prim
23 23
24 struct point : public prim::point 24struct point : public prim::point
25 { 25{
26 zone zone; 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 27
34 struct linestring : public prim::linestring 28 explicit inline point(class zone zone, double x, double y)
29 : prim::point{x, y}, zone{zone}
35 { 30 {
36 zone zone; 31 }
32};
33
34struct linestring : public prim::linestring
35{
36 zone zone;
37 37
38 explicit linestring(class zone zone) 38 explicit inline linestring(class zone zone) : zone{zone} {}
39 : zone{zone} 39};
40 {
41 }
42 };
43 40
44} // namespace routemon::geo::utm::zone_local 41} // namespace routemon::geo::utm::zone_local
45 42
@@ -47,28 +44,46 @@ BOOST_GEOMETRY_REGISTER_LINESTRING(routemon::geo::utm::zone_local::linestring)
47 44
48namespace boost::geometry::traits { 45namespace boost::geometry::traits {
49 46
50 namespace { 47namespace {
51 48
52 using zone_local_point = routemon::geo::utm::zone_local::point; 49using zone_local_point = routemon::geo::utm::zone_local::point;
53 50
54 } // namespace <anonymous> 51} // namespace
55 52
56 template<> struct tag<zone_local_point> { using type = point_tag; }; 53template <>
57 template<> struct dimension<zone_local_point> : boost::mpl::int_<2> {}; 54struct tag<zone_local_point>
58 template<> struct coordinate_type<zone_local_point> { using type = double; }; 55{
59 template<> struct coordinate_system<zone_local_point> { using type = routemon::geo::utm::zone_local::cs; }; 56 using type = point_tag;
57};
58template <>
59struct dimension<zone_local_point> : boost::mpl::int_<2>
60{
61};
62template <>
63struct coordinate_type<zone_local_point>
64{
65 using type = double;
66};
67template <>
68struct coordinate_system<zone_local_point>
69{
70 using type = routemon::geo::utm::zone_local::cs;
71};
60 72
61 template<std::size_t index> 73template <std::size_t index>
62 struct access<zone_local_point, index> { 74struct access<zone_local_point, index>
63 static inline auto get(zone_local_point const& p) -> double 75{
64 { 76 static inline auto get(zone_local_point const& p) -> double
65 return bgeo::get<index>(static_cast<routemon::geo::utm::zone_local::prim::point const&>(p)); 77 {
66 } 78 return bgeo::get<index>(
79 static_cast<routemon::geo::utm::zone_local::prim::point const&>(p));
80 }
67 81
68 static inline auto set(zone_local_point& p, double v) -> void 82 static inline auto set(zone_local_point& p, double v) -> void
69 { 83 {
70 return bgeo::set<index>(static_cast<routemon::geo::utm::zone_local::prim::point&>(p), v); 84 return bgeo::set<index>(
71 } 85 static_cast<routemon::geo::utm::zone_local::prim::point&>(p), v);
72 }; 86 }
87};
73 88
74} // namespace boost::geometry::traits 89} // namespace boost::geometry::traits
diff --git a/server/src/geo/wgs84.cppm b/server/src/geo/wgs84.cppm
index d574d39..b53a963 100644
--- a/server/src/geo/wgs84.cppm
+++ b/server/src/geo/wgs84.cppm
@@ -13,25 +13,20 @@ using point = bgeo::model::point<double, 2, cs>;
13// TODO: guarantee that linestring provides non-static member 13// TODO: guarantee that linestring provides non-static member
14// auto reserve(std::size_t) -> void 14// auto reserve(std::size_t) -> void
15// Perhaps better yet: guarantee that the backing container is a std::vector. 15// Perhaps better yet: guarantee that the backing container is a std::vector.
16using linestring = bgeo::model::linestring<point, std::vector /* the default */>; 16using linestring =
17 bgeo::model::linestring<point, std::vector /* the default */>;
17using box = bgeo::model::box<point>; 18using box = bgeo::model::box<point>;
18using stype = bgeo::srs::spheroid<double>; 19using stype = bgeo::srs::spheroid<double>;
19using vincenty_strategy = bgeo::strategy::distance::vincenty<stype>; 20using vincenty_strategy = bgeo::strategy::distance::vincenty<stype>;
20 21
21auto from_lat_lon(double lat, double lon) -> point 22inline auto from_lat_lon(double lat, double lon) -> point
22{ 23{
23 return point{lon, lat}; 24 return point{lon, lat};
24} 25}
25 26
26auto lat(point p) -> double 27inline auto lat(point p) -> double { return bgeo::get<1>(p); }
27{
28 return bgeo::get<1>(p);
29}
30 28
31auto lon(point p) -> double 29inline auto lon(point p) -> double { return bgeo::get<0>(p); }
32{
33 return bgeo::get<0>(p);
34}
35 30
36// Point p with latitude in [-90, 90] and longitude in [-180, 180) 31// Point p with latitude in [-90, 90] and longitude in [-180, 180)
37auto normalize(point p) -> point 32auto normalize(point p) -> point
@@ -44,7 +39,7 @@ auto normalize(point p) -> point
44 auto p_lat = nonneg_mod360(lat(p)) - 90.0; 39 auto p_lat = nonneg_mod360(lat(p)) - 90.0;
45 if (90.0 < p_lat) 40 if (90.0 < p_lat)
46 { 41 {
47 assert(p_lat < 270.0); // by nonneg_mod360 42 assert(p_lat < 270.0); // by nonneg_mod360
48 p_lon += 180.0; 43 p_lon += 180.0;
49 p_lat -= 180.0; 44 p_lat -= 180.0;
50 } 45 }
@@ -53,7 +48,7 @@ auto normalize(point p) -> point
53 return from_lat_lon(p_lat, p_lon); 48 return from_lat_lon(p_lat, p_lon);
54} 49}
55 50
56auto is_normalized(point p) -> bool 51inline auto is_normalized(point p) -> bool
57{ 52{
58 return -90 <= lat(p) && lat(p) <= 90 && -180 <= lon(p) && lon(p) < 180; 53 return -90 <= lat(p) && lat(p) <= 90 && -180 <= lon(p) && lon(p) < 180;
59} 54}
@@ -63,32 +58,22 @@ class normalized_point
63 point p_; 58 point p_;
64 59
65public: 60public:
66 normalized_point(point p) 61 inline normalized_point(point p) : p_{is_normalized(p) ? p : normalize(p)} {}
67 : p_{is_normalized(p) ? p : normalize(p)}
68 {}
69 62
70 normalized_point() 63 inline normalized_point() : p_{} {}
71 : p_{}
72 {}
73 64
74 auto lat() const -> double 65 inline auto lat() const -> double { return bgeo::get<1>(p_); }
75 {
76 return bgeo::get<1>(p_);
77 }
78 66
79 auto lon() const -> double 67 inline auto lon() const -> double { return bgeo::get<0>(p_); }
80 {
81 return bgeo::get<0>(p_);
82 }
83 68
84 auto lat(double lat) -> void 69 inline auto lat(double lat) -> void
85 { 70 {
86 if (lat < -90 || lat > 90) 71 if (lat < -90 || lat > 90)
87 throw std::range_error{"latitude out of range"}; 72 throw std::range_error{"latitude out of range"};
88 bgeo::set<1>(p_, lat); 73 bgeo::set<1>(p_, lat);
89 } 74 }
90 75
91 auto lon(double lon) -> void 76 inline auto lon(double lon) -> void
92 { 77 {
93 if (lon < -180 || lon > 180) 78 if (lon < -180 || lon > 180)
94 throw std::range_error{"longitude out of range"}; 79 throw std::range_error{"longitude out of range"};
@@ -101,15 +86,9 @@ class normalized_linestring
101 std::vector<normalized_point> ls_; 86 std::vector<normalized_point> ls_;
102 87
103public: 88public:
104 auto empty() const -> bool 89 inline auto empty() const -> bool { return ls_.empty(); }
105 {
106 return ls_.empty();
107 }
108 90
109 auto push_back(normalized_point p) -> void 91 inline auto push_back(normalized_point p) -> void { ls_.push_back(p); }
110 {
111 ls_.push_back(p);
112 }
113}; 92};
114 93
115} // namespace routemon::geo::wgs84 94} // namespace routemon::geo::wgs84
generated by cgit v1.3 (git 2.53.0) at 2026-09-10 14:04:40 +0000