From 1417e6cf41d27c594d2bdc92f87746d2383a1b1d Mon Sep 17 00:00:00 2001 From: Rutger Broekhoff Date: Tue, 8 Sep 2026 13:54:40 +0200 Subject: Some cleanups --- server/src/api.cpp | 3 +- server/src/api.cppm | 6 +- server/src/geo/multizonal.cppm | 307 --------------------------------- server/src/geo/utm.cppm | 140 +++++++-------- server/src/geo/utm_multizonal.cppm | 343 +++++++++++++++++++++++++++++++++++++ server/src/geo/utm_zone_local.cppm | 103 ++++++----- server/src/geo/wgs84.cppm | 51 ++---- 7 files changed, 494 insertions(+), 459 deletions(-) delete mode 100644 server/src/geo/multizonal.cppm create mode 100644 server/src/geo/utm_multizonal.cppm (limited to 'server/src') 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) { l_.debug("Checking part [{}/{}]", ++i, segments.size()); - auto part_zone_lss = geo::multizonal::zoned_linestring_seg_seq{part}; + auto part_zone_lss = + geo::utm::multizonal::split_linestring_across_zones(part); // TODO: consider buffering with min_distance_ auto part_box = geo::wgs84::box{}; 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; import std; import :datex2; import :geo.wgs84; -import :geo.multizonal; +import :geo.utm.multizonal; import :gpx; import :log; import :time; @@ -81,9 +81,9 @@ class handler constexpr static auto const min_distance_ = 5.0; using blse_index_value = std::shared_ptr; - using blse_index = geo::multizonal::linestring_rtree; + using blse_index = geo::utm::multizonal::linestring_rtree; using bpe_index_value = std::shared_ptr; - using bpe_index = geo::multizonal::point_rtree; + using bpe_index = geo::utm::multizonal::point_rtree; log::logger l_; 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 @@ -module; - -#include -#include -#include - -export module routemon:geo.multizonal; - -import std; -import :geo; -import :geo.utm; -import :geo.utm.zone_local; -import :geo.wgs84; - -namespace routemon::geo::multizonal { - -template -class multi_zone -{ - std::array zones_; - -public: - auto operator[](utm::zone z) -> T& - { - return zones_[z.as_index()]; - } - - auto operator[](utm::zone z) const -> T const& - { - return zones_[z.as_index()]; - } -}; - -namespace wgs84 -{ - -// Transformations from WGS 84 (EPSG:4326) -using transform_from_t = bgeo::srs::transformation>; - -class utm_transform : public transform_from_t -{ - utm::zone to_zone_; - -public: - explicit utm_transform(utm::zone to_zone) - : transform_from_t{{}, bgeo::srs::epsg{to_zone.wgs84_proj_epsg()}}, - to_zone_{to_zone} - { - } - - utm_transform() - : utm_transform{utm::zone::min()} - {} - - auto apply(utm::zonable_wgs84_point p) -> utm::zone_local::point - { - auto local_p = utm::zone_local::point{to_zone_, 0.0, 0.0}; - forward(p, local_p); - return local_p; - } -}; - -class utm_transforms : public multi_zone -{ - utm_transforms() - { - for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next()) - { - (*this)[z] = utm_transform{z}; - } - } - -public: - static utm_transforms const& instance() - { - static utm_transforms inst; - return inst; - } -}; - -auto neighbor_utm_zone(utm::zonable_wgs84_point p) -> std::pair -{ - auto separating_meridian_lon = std::round(p.lon() / 6.0) * 6.0; - auto closest_zone_middle = separating_meridian_lon < p.lon() - ? separating_meridian_lon - 3.0 - : separating_meridian_lon + 3.0; - 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()); - - auto separating_meridian_ls = geo::wgs84::linestring{ - geo::wgs84::from_lat_lon(90.0, 0.0), - geo::wgs84::from_lat_lon(0.0, separating_meridian_lon), - geo::wgs84::from_lat_lon(-90.0, 0.0), - }; - auto closest_zone_dist = bgeo::distance(p, separating_meridian_ls, geo::wgs84::vincenty_strategy{}); - - return std::make_pair(closest_zone, closest_zone_dist); -} - -} // namespace wgs84 - -struct multizone_linestring -{ - multi_zone> segments; - - // TODO: consider taking an input range instead - explicit multizone_linestring(utm::zonable_wgs84_linestring const& ls) - { - auto to_utm = wgs84::utm_transforms::instance(); - auto working = multi_zone{}; - auto mprev_p_zone = std::optional{}; - auto mprev_p_alt_zone = std::optional{}; - - auto push = [&](utm::zonable_wgs84_point p, utm::zone z) - { working[z].push_back(to_utm[z].apply(p)); }; - auto flush = [&](utm::zone z) - { - segments[z].push_back(std::move(working[z])); - working[z] = {}; - }; - - for (auto const& p : ls) - { - auto p_zone = utm::zone::for_wgs84_point(p); - auto mp_alt_zone = std::optional{}; - if (auto [neighbor_zone, neighbor_zone_dist] = wgs84::neighbor_utm_zone(p); - neighbor_zone_dist < 30.0 /* m */) - mp_alt_zone = neighbor_zone; - assert(!mp_alt_zone || *mp_alt_zone != p_zone); - - push(p, p_zone); - if (mp_alt_zone) - push(p, *mp_alt_zone); - if (mprev_p_zone && *mprev_p_zone != p_zone && mprev_p_zone != mp_alt_zone) - { - push(p, *mprev_p_zone); - flush(*mprev_p_zone); - } - if (mprev_p_alt_zone && *mprev_p_alt_zone != p_zone && mprev_p_alt_zone != mp_alt_zone) - { - push(p, *mprev_p_alt_zone); - flush(*mprev_p_alt_zone); - } - - mprev_p_zone = p_zone; - mprev_p_alt_zone = mp_alt_zone; - } - - if (mprev_p_zone && !working[*mprev_p_zone].empty()) - { - flush(*mprev_p_zone); - } - if (mprev_p_alt_zone && !working[*mprev_p_alt_zone].empty()) - { - flush(*mprev_p_alt_zone); - } - } -}; - -struct zoned_linestring_seg_seq -{ - std::vector segments; - - explicit zoned_linestring_seg_seq(utm::zonable_wgs84_linestring ls) - { - auto to_utm = wgs84::utm_transforms::instance(); - auto mworking_seg = std::optional{}; - - for (auto const& p : ls) - { - auto p_zone = utm::zone::for_wgs84_point(p); - if (mworking_seg && mworking_seg->zone != p_zone) - { - mworking_seg->push_back(to_utm[mworking_seg->zone].apply(p)); - segments.push_back(std::move(*mworking_seg)); - mworking_seg = std::nullopt; - } - if (!mworking_seg) - mworking_seg = utm::zone_local::linestring{p_zone}; - mworking_seg->push_back(to_utm[p_zone].apply(p)); - } - - if (mworking_seg && !mworking_seg->empty()) - { - segments.push_back(std::move(*mworking_seg)); - } - } -}; - -template -class linestring_rtree -{ -public: - using index_value = std::tuple; - -private: - multi_zone>> local_rtrees_; - -public: - template U> - auto insert(utm::zonable_wgs84_linestring const& ls, U&& arg) -> void - { - auto zone_segments = multizone_linestring{ls}.segments; - for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next()) - { - for (auto ls : zone_segments[z]) - { - auto blse = - bgeo::return_buffer(bgeo::return_envelope(ls), 30.0 /* m */); - local_rtrees_[z].insert(index_value{blse, std::move(ls), std::forward(arg)}); - } - } - } - - auto size() const -> std::size_t - { - auto total_size = 0uz; - for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next()) - total_size += local_rtrees_[z].size(); - return total_size; - } - - // Note: the same T may be generated more than once! - auto intersection(zoned_linestring_seg_seq const& lss) const -> std::generator - { - for (auto const& ls : lss.segments) - { - // auto ls_box = bgeo::return_envelope(ls); - for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls)); - it != local_rtrees_[ls.zone].qend(); it++) - { - using multi_linestring = bgeo::model::multi_linestring; - - auto relevant_ls_parts = multi_linestring{}; - bgeo::intersection(std::get<0>(*it), static_cast(ls), relevant_ls_parts); - if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */) - { - co_yield std::get<2>(*it); - } - } - } - } -}; - -template -class point_rtree -{ -public: - using index_value = std::tuple; - -private: - multi_zone>> local_rtrees_; - -public: - template U> - auto insert(utm::zonable_wgs84_point const& p, U&& arg) -> void - { - auto to_utm = wgs84::utm_transforms::instance(); - - auto p_zone = utm::zone::for_wgs84_point(p); - auto mp_alt_zone = std::optional{}; - if (auto [neighbor_zone, neighbor_zone_dist] = wgs84::neighbor_utm_zone(p); - neighbor_zone_dist < 30.0 /* m */) - mp_alt_zone = neighbor_zone; - - auto insert = [&](utm::zone z) -> void - { - auto p_utm = to_utm[z].apply(p); - auto bpe = bgeo::return_buffer(bgeo::return_envelope(p_utm), 30.0 /* m */); - local_rtrees_[z].insert(index_value{bpe, p_utm, arg}); - }; - - insert(p_zone); - if (mp_alt_zone) - insert(*mp_alt_zone); - } - - auto size() const -> std::size_t - { - auto total_size = 0uz; - for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next()) - total_size += local_rtrees_[z].size(); - return total_size; - } - - // Note: the same T may be generated more than once! - auto intersection(zoned_linestring_seg_seq const& lss) const -> std::generator - { - for (auto const& ls : lss.segments) - { - // auto ls_box = bgeo::return_envelope(ls); - for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls)); - it != local_rtrees_[ls.zone].qend(); it++) - { - using multi_linestring = bgeo::model::multi_linestring; - - auto relevant_ls_parts = multi_linestring{}; - bgeo::intersection(std::get<0>(*it), static_cast(ls), relevant_ls_parts); - if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */) - { - co_yield std::get<2>(*it); - } - } - } - } -}; - -} // 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 @@ module; -#include #include +#include export module routemon:geo.utm; @@ -9,46 +9,32 @@ import std; import :geo; import :geo.wgs84; -namespace routemon::geo::utm -{ +namespace routemon::geo::utm { class zonable_wgs84_point { wgs84::normalized_point p_; - explicit zonable_wgs84_point(wgs84::normalized_point p) - : p_{p} - {} + explicit inline zonable_wgs84_point(wgs84::normalized_point p) : p_{p} {} public: - zonable_wgs84_point() - : p_{} - { - } + inline zonable_wgs84_point() : p_{} {} - static auto from(wgs84::normalized_point p) -> std::optional + inline static auto from(wgs84::normalized_point p) + -> std::optional { if (p.lat() < -80 || p.lat() > 84) return std::nullopt; return zonable_wgs84_point{p}; } - auto lat() const -> double - { - return p_.lat(); - } + inline auto lat() const -> double { return p_.lat(); } - auto lon() const -> double - { - return p_.lon(); - } + inline auto lon() const -> double { return p_.lon(); } - auto lon(double lon) -> void - { - p_.lon(lon); - } + inline auto lon(double lon) -> void { p_.lon(lon); } - auto lat(double lat) -> void + inline auto lat(double lat) -> void { if (lat < -80 || lat > 84) throw std::range_error{"latitude out of range for UTM"}; @@ -56,48 +42,69 @@ public: } }; -} +} // namespace routemon::geo::utm namespace boost::geometry::traits { - namespace { +namespace { - using zonable_wgs84_point = routemon::geo::utm::zonable_wgs84_point; +using zonable_wgs84_point = routemon::geo::utm::zonable_wgs84_point; - } // namespace +} // namespace - template<> struct tag { using type = point_tag; }; - template<> struct dimension : boost::mpl::int_<2> {}; - template<> struct coordinate_type { using type = double; }; - template<> struct coordinate_system { using type = routemon::geo::wgs84::cs; }; +template <> +struct tag +{ + using type = point_tag; +}; - template - struct access { - static inline auto get(zonable_wgs84_point const& p) -> double - { - if constexpr (index == 0) - return p.lon(); - else if constexpr (index == 1) - return p.lat(); - else static_assert(false, "Out of range"); - } +template <> +struct dimension : boost::mpl::int_<2> +{ +}; - static inline auto set(zonable_wgs84_point& p, double v) -> void - { - if constexpr (index == 0) - p.lon(v); - else if constexpr (index == 1) - p.lat(v); - else static_assert(false, "Out of range"); - } - }; +template <> +struct coordinate_type +{ + using type = double; +}; -} // namespace boost::geometry::traits +template <> +struct coordinate_system +{ + using type = routemon::geo::wgs84::cs; +}; -namespace routemon::geo::utm +template +struct access { + static inline auto get(zonable_wgs84_point const& p) -> double + { + if constexpr (index == 0) + return p.lon(); + else if constexpr (index == 1) + return p.lat(); + else + static_assert(false, "Out of range"); + } + + static inline auto set(zonable_wgs84_point& p, double v) -> void + { + if constexpr (index == 0) + p.lon(v); + else if constexpr (index == 1) + p.lat(v); + else + static_assert(false, "Out of range"); + } +}; + +} // namespace boost::geometry::traits -using zonable_wgs84_linestring = bgeo::model::linestring; +namespace routemon::geo::utm { + +using zonable_wgs84_linestring = + bgeo::model::linestring; // In the sense of the common WGS84 subdivisions by simple northing // and easting (so no Norway/Svalbard exceptions). @@ -130,33 +137,30 @@ class zone public: explicit constexpr zone(std::uint8_t zone_no, hemisphere h) : zone_{from(zone_no, h)} - {} + { + } - auto hemisphere() const -> enum hemisphere + constexpr auto hemisphere() const -> enum hemisphere { return static_cast(zone_ % 2); } // Return value in range [1, 60] - auto zone_no() const -> std::uint8_t - { - return 1 + zone_ / 2; - } + constexpr auto zone_no() const -> std::uint8_t { return 1 + zone_ / 2; } // Return value in range [min(), max()] - constexpr auto as_index() const -> std::uint8_t - { - return zone_; - } + constexpr auto as_index() const -> std::uint8_t { return zone_; } static auto for_wgs84_point(zonable_wgs84_point p) noexcept -> zone { auto zone_no = static_cast(1 + (p.lon() + 180.0) / 6.0); auto northern = p.lat() >= 0.0; - return zone{zone_no, northern ? hemisphere::northern : hemisphere::southern}; + return zone{ + zone_no, northern ? hemisphere::northern : hemisphere::southern + }; } - auto wgs84_proj_epsg() const -> int + constexpr auto wgs84_proj_epsg() const -> int { switch (hemisphere()) { @@ -178,7 +182,7 @@ public: return zone{60, hemisphere::southern}; } - auto next() const -> zone + constexpr auto next() const -> zone { assert(zone_ <= max().as_index()); if (zone_ == max().as_index()) @@ -188,7 +192,7 @@ public: return copy; } - auto operator==(zone rhs) const -> bool + constexpr auto operator==(zone rhs) const -> bool { return zone_ == rhs.zone_; } 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 @@ +module; + +#include +#include +#include + +export module routemon:geo.utm.multizonal; + +import std; +import :geo; +import :geo.utm; +import :geo.utm.zone_local; +import :geo.wgs84; + +namespace routemon::geo::utm::multizonal { + +template +class multi_zone +{ + std::array zones_; + +public: + auto operator[](zone z) -> T& { return zones_[z.as_index()]; } + + auto operator[](zone z) const -> T const& { return zones_[z.as_index()]; } +}; + +// Transformations from WGS 84 (EPSG:4326) +using from_wgs84_transform_base = + bgeo::srs::transformation>; + +class from_wgs84_transform : public from_wgs84_transform_base +{ + zone to_zone_; + +public: + explicit from_wgs84_transform(zone to_zone) + : from_wgs84_transform_base{{}, bgeo::srs::epsg{to_zone.wgs84_proj_epsg()}}, + to_zone_{to_zone} + { + } + + from_wgs84_transform() : from_wgs84_transform{zone::min()} {} + + auto apply(zonable_wgs84_point p) -> zone_local::point + { + auto local_p = zone_local::point{to_zone_, 0.0, 0.0}; + forward(p, local_p); + return local_p; + } +}; + +class from_wgs84_transforms : public multi_zone +{ + from_wgs84_transforms() + { + for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next()) + { + (*this)[z] = from_wgs84_transform{z}; + } + } + +public: + static from_wgs84_transforms const& instance() + { + static from_wgs84_transforms inst; + return inst; + } +}; + +auto neighbor_utm_zone(utm::zonable_wgs84_point p) + -> std::pair +{ + auto separating_meridian_lon = std::round(p.lon() / 6.0) * 6.0; + auto closest_zone_middle = separating_meridian_lon < p.lon() + ? separating_meridian_lon - 3.0 + : separating_meridian_lon + 3.0; + 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()); + + auto separating_meridian_ls = geo::wgs84::linestring{ + geo::wgs84::from_lat_lon(90.0, 0.0), + geo::wgs84::from_lat_lon(0.0, separating_meridian_lon), + geo::wgs84::from_lat_lon(-90.0, 0.0), + }; + auto closest_zone_dist = bgeo::distance( + p, separating_meridian_ls, geo::wgs84::vincenty_strategy{}); + + return std::make_pair(closest_zone, closest_zone_dist); +} + +// Distribute a WGS 84 linestring over the UTM zones that it finds +// itself in or near to. +// +// Since it is not at all unthinkable that some road works will end up +// crossing different UTM zones, we need some machinery to work with +// such situations. Since we want to report all situations within a +// specified distance of the planned route, we must also consider the +// situation where this is indeed the case, but the situation is +// situated at the opposite side of a delineating UTM zone meridian +// w.r.t. the planned route. +// +// To solve this issue, we create segments of the input line string +// for multiple UTM zones when crossing zone boundaries or when very +// close to zone boundaries. +auto distribute_linestring_over_zones(utm::zonable_wgs84_linestring const& ls) + -> multi_zone> +{ + auto result = multi_zone>{}; + + auto to_utm = from_wgs84_transforms::instance(); + auto working = multi_zone{}; + auto mprev_p_zone = std::optional{}; + auto mprev_p_alt_zone = std::optional{}; + + auto push = [&](utm::zonable_wgs84_point p, utm::zone z) + { working[z].push_back(to_utm[z].apply(p)); }; + auto flush = [&](utm::zone z) + { + result[z].push_back(std::move(working[z])); + working[z] = {}; + }; + + for (auto const& p : ls) + { + auto p_zone = utm::zone::for_wgs84_point(p); + auto mp_alt_zone = std::optional{}; + if (auto [neighbor_zone, neighbor_zone_dist] = neighbor_utm_zone(p); + neighbor_zone_dist < 30.0 /* m */) + mp_alt_zone = neighbor_zone; + assert(!mp_alt_zone || *mp_alt_zone != p_zone); + + push(p, p_zone); + if (mp_alt_zone) + push(p, *mp_alt_zone); + if (mprev_p_zone && *mprev_p_zone != p_zone && mprev_p_zone != mp_alt_zone) + { + push(p, *mprev_p_zone); + flush(*mprev_p_zone); + } + if (mprev_p_alt_zone && *mprev_p_alt_zone != p_zone + && mprev_p_alt_zone != mp_alt_zone) + { + push(p, *mprev_p_alt_zone); + flush(*mprev_p_alt_zone); + } + + mprev_p_zone = p_zone; + mprev_p_alt_zone = mp_alt_zone; + } + + if (mprev_p_zone && !working[*mprev_p_zone].empty()) + { + flush(*mprev_p_zone); + } + if (mprev_p_alt_zone && !working[*mprev_p_alt_zone].empty()) + { + flush(*mprev_p_alt_zone); + } + + return result; +} + +// Split a WGS 84 line string into a sequence of UTM-zone-local line +// strings, starting a new UTM-zone-local line string when the input +// line string crosses a UTM zone boundary. The line segment that +// crosses the zone boundary can be found in the UTM-zone-local line +// strings for both zones which its points are in (we assume that line +// segments will be short enough to not cause significant distortion +// here, so we do not put in the effort to clip at the +// zone-delineating meridian here). +auto split_linestring_across_zones(utm::zonable_wgs84_linestring ls) + -> std::vector +{ + auto splits = std::vector{}; + + auto to_utm = from_wgs84_transforms::instance(); + auto mworking_seg = std::optional{}; + + for (auto const& p : ls) + { + auto p_zone = utm::zone::for_wgs84_point(p); + if (mworking_seg && mworking_seg->zone != p_zone) + { + mworking_seg->push_back(to_utm[mworking_seg->zone].apply(p)); + splits.push_back(std::move(*mworking_seg)); + mworking_seg = std::nullopt; + } + if (!mworking_seg) + mworking_seg = utm::zone_local::linestring{p_zone}; + mworking_seg->push_back(to_utm[p_zone].apply(p)); + } + + if (mworking_seg && !mworking_seg->empty()) + { + splits.push_back(std::move(*mworking_seg)); + } + + return splits; +} + +template +class linestring_rtree +{ +public: + using index_value = std:: + tuple; + +private: + multi_zone>> + local_rtrees_; + +public: + template U> + auto insert(utm::zonable_wgs84_linestring const& ls, U&& arg) -> void + { + auto zone_segments = distribute_linestring_over_zones(ls); + for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next()) + { + for (auto ls : zone_segments[z]) + { + auto blse = bgeo::return_buffer( + bgeo::return_envelope(ls), + 30.0 /* m */); + local_rtrees_[z].insert( + index_value{blse, std::move(ls), std::forward(arg)}); + } + } + } + + auto size() const -> std::size_t + { + auto total_size = 0uz; + for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next()) + total_size += local_rtrees_[z].size(); + return total_size; + } + + // Note: the same T may be generated more than once! + auto intersection(std::vector const& lss) const + -> std::generator + { + for (auto const& ls : lss) + { + for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls)); + it != local_rtrees_[ls.zone].qend(); it++) + { + using multi_linestring = + bgeo::model::multi_linestring; + + auto relevant_ls_parts = multi_linestring{}; + bgeo::intersection( + std::get<0>(*it), + static_cast(ls), + relevant_ls_parts); + // it is somewhat awkward and arbitrary that this measurement is + // performed here + if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */) + { + co_yield std::get<2>(*it); + } + } + } + } +}; + +template +class point_rtree +{ +public: + using index_value = + std::tuple; + +private: + multi_zone>> + local_rtrees_; + +public: + template U> + auto insert(utm::zonable_wgs84_point const& p, U&& arg) -> void + { + auto to_utm = from_wgs84_transforms::instance(); + + auto p_zone = utm::zone::for_wgs84_point(p); + auto mp_alt_zone = std::optional{}; + if (auto [neighbor_zone, neighbor_zone_dist] = neighbor_utm_zone(p); + neighbor_zone_dist < 30.0 /* m */) + mp_alt_zone = neighbor_zone; + + auto insert = [&](utm::zone z) -> void + { + auto p_utm = to_utm[z].apply(p); + auto bpe = bgeo::return_buffer( + bgeo::return_envelope(p_utm), + 30.0 /* m */); + local_rtrees_[z].insert(index_value{bpe, p_utm, arg}); + }; + + insert(p_zone); + if (mp_alt_zone) + insert(*mp_alt_zone); + } + + auto size() const -> std::size_t + { + auto total_size = 0uz; + for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next()) + total_size += local_rtrees_[z].size(); + return total_size; + } + + // Note: the same T may be generated more than once! + auto intersection(std::vector const& lss) const + -> std::generator + { + for (auto const& ls : lss) + { + for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls)); + it != local_rtrees_[ls.zone].qend(); it++) + { + using multi_linestring = + bgeo::model::multi_linestring; + + auto relevant_ls_parts = multi_linestring{}; + bgeo::intersection( + std::get<0>(*it), + static_cast(ls), + relevant_ls_parts); + // it is somewhat awkward and arbitrary that this measurement is + // performed here + if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */) + { + co_yield std::get<2>(*it); + } + } + } + } +}; + +} // 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; // Coordinates that are local to some known UTM zone. namespace routemon::geo::utm::zone_local { - using cs = bgeo::cs::cartesian; +using cs = bgeo::cs::cartesian; - namespace prim { +namespace prim { - using point = bgeo::model::d2::point_xy; - using linestring = bgeo::model::linestring; - using box = bgeo::model::box; +using point = bgeo::model::d2::point_xy; +using linestring = bgeo::model::linestring; +using box = bgeo::model::box; - } // namespace prim +} // namespace prim - struct point : public prim::point - { - zone zone; - - explicit point(class zone zone, double x, double y) - : prim::point{x, y}, zone{zone} - { - } - }; +struct point : public prim::point +{ + zone zone; - struct linestring : public prim::linestring + explicit inline point(class zone zone, double x, double y) + : prim::point{x, y}, zone{zone} { - zone zone; + } +}; + +struct linestring : public prim::linestring +{ + zone zone; - explicit linestring(class zone zone) - : zone{zone} - { - } - }; + explicit inline linestring(class zone zone) : zone{zone} {} +}; } // namespace routemon::geo::utm::zone_local @@ -47,28 +44,46 @@ BOOST_GEOMETRY_REGISTER_LINESTRING(routemon::geo::utm::zone_local::linestring) namespace boost::geometry::traits { - namespace { - - using zone_local_point = routemon::geo::utm::zone_local::point; - - } // namespace - - template<> struct tag { using type = point_tag; }; - template<> struct dimension : boost::mpl::int_<2> {}; - template<> struct coordinate_type { using type = double; }; - template<> struct coordinate_system { using type = routemon::geo::utm::zone_local::cs; }; - - template - struct access { - static inline auto get(zone_local_point const& p) -> double - { - return bgeo::get(static_cast(p)); - } +namespace { + +using zone_local_point = routemon::geo::utm::zone_local::point; + +} // namespace + +template <> +struct tag +{ + using type = point_tag; +}; +template <> +struct dimension : boost::mpl::int_<2> +{ +}; +template <> +struct coordinate_type +{ + using type = double; +}; +template <> +struct coordinate_system +{ + using type = routemon::geo::utm::zone_local::cs; +}; + +template +struct access +{ + static inline auto get(zone_local_point const& p) -> double + { + return bgeo::get( + static_cast(p)); + } - static inline auto set(zone_local_point& p, double v) -> void - { - return bgeo::set(static_cast(p), v); - } - }; + static inline auto set(zone_local_point& p, double v) -> void + { + return bgeo::set( + static_cast(p), v); + } +}; } // 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; // TODO: guarantee that linestring provides non-static member // auto reserve(std::size_t) -> void // Perhaps better yet: guarantee that the backing container is a std::vector. -using linestring = bgeo::model::linestring; +using linestring = + bgeo::model::linestring; using box = bgeo::model::box; using stype = bgeo::srs::spheroid; using vincenty_strategy = bgeo::strategy::distance::vincenty; -auto from_lat_lon(double lat, double lon) -> point +inline auto from_lat_lon(double lat, double lon) -> point { return point{lon, lat}; } -auto lat(point p) -> double -{ - return bgeo::get<1>(p); -} +inline auto lat(point p) -> double { return bgeo::get<1>(p); } -auto lon(point p) -> double -{ - return bgeo::get<0>(p); -} +inline auto lon(point p) -> double { return bgeo::get<0>(p); } // Point p with latitude in [-90, 90] and longitude in [-180, 180) auto normalize(point p) -> point @@ -44,7 +39,7 @@ auto normalize(point p) -> point auto p_lat = nonneg_mod360(lat(p)) - 90.0; if (90.0 < p_lat) { - assert(p_lat < 270.0); // by nonneg_mod360 + assert(p_lat < 270.0); // by nonneg_mod360 p_lon += 180.0; p_lat -= 180.0; } @@ -53,7 +48,7 @@ auto normalize(point p) -> point return from_lat_lon(p_lat, p_lon); } -auto is_normalized(point p) -> bool +inline auto is_normalized(point p) -> bool { return -90 <= lat(p) && lat(p) <= 90 && -180 <= lon(p) && lon(p) < 180; } @@ -63,32 +58,22 @@ class normalized_point point p_; public: - normalized_point(point p) - : p_{is_normalized(p) ? p : normalize(p)} - {} + inline normalized_point(point p) : p_{is_normalized(p) ? p : normalize(p)} {} - normalized_point() - : p_{} - {} + inline normalized_point() : p_{} {} - auto lat() const -> double - { - return bgeo::get<1>(p_); - } + inline auto lat() const -> double { return bgeo::get<1>(p_); } - auto lon() const -> double - { - return bgeo::get<0>(p_); - } + inline auto lon() const -> double { return bgeo::get<0>(p_); } - auto lat(double lat) -> void + inline auto lat(double lat) -> void { if (lat < -90 || lat > 90) throw std::range_error{"latitude out of range"}; bgeo::set<1>(p_, lat); } - auto lon(double lon) -> void + inline auto lon(double lon) -> void { if (lon < -180 || lon > 180) throw std::range_error{"longitude out of range"}; @@ -101,15 +86,9 @@ class normalized_linestring std::vector ls_; public: - auto empty() const -> bool - { - return ls_.empty(); - } + inline auto empty() const -> bool { return ls_.empty(); } - auto push_back(normalized_point p) -> void - { - ls_.push_back(p); - } + inline auto push_back(normalized_point p) -> void { ls_.push_back(p); } }; } // namespace routemon::geo::wgs84 -- cgit v1.3