From 3e8ce840c70ee00925210b96738696d1dc445f8f Mon Sep 17 00:00:00 2001 From: Rutger Broekhoff Date: Tue, 8 Sep 2026 01:36:38 +0200 Subject: UTM projection --- server/src/api.cpp | 86 ++++--------- server/src/api.cppm | 19 ++- server/src/datex2.cpp | 62 +++++++--- server/src/datex2.cppm | 12 +- server/src/geo.cpp | 42 ------- server/src/geo.cppm | 23 ---- server/src/geo/geo.cppm | 7 ++ server/src/geo/multizonal.cppm | 240 +++++++++++++++++++++++++++++++++++++ server/src/geo/utm.cppm | 199 ++++++++++++++++++++++++++++++ server/src/geo/utm_zone_local.cppm | 74 ++++++++++++ server/src/geo/wgs84.cppm | 115 ++++++++++++++++++ server/src/gpx.cpp | 16 ++- server/src/gpx.cppm | 4 +- server/src/log.cpp | 2 +- 14 files changed, 738 insertions(+), 163 deletions(-) delete mode 100644 server/src/geo.cpp delete mode 100644 server/src/geo.cppm create mode 100644 server/src/geo/geo.cppm create mode 100644 server/src/geo/multizonal.cppm create mode 100644 server/src/geo/utm.cppm create mode 100644 server/src/geo/utm_zone_local.cppm create mode 100644 server/src/geo/wgs84.cppm (limited to 'server/src') 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; namespace routemon::api { -auto json_value_from_point(geo::point const& p) -> json::value +auto json_value_from_point(geo::utm::zonable_wgs84_point const& p) + -> json::value { - return json::array{bgeo::get<1>(p), bgeo::get<0>(p)}; + return json::array{p.lat(), p.lon()}; } -auto json_value_from_linestring(geo::linestring const& ls) -> json::value +auto json_value_from_linestring(geo::utm::zonable_wgs84_linestring const& ls) + -> json::value { json::array a; for (auto const& p : ls) @@ -30,8 +32,8 @@ auto json_value_from_linestring(geo::linestring const& ls) -> json::value return a; } -auto json_value_from_linestrings(std::vector const& lss) - -> json::value +auto json_value_from_linestrings( + std::vector const& lss) -> json::value { json::array a; for (auto const& ls : lss) @@ -108,24 +110,7 @@ handler::handler(log::logger const& l, datex2::situation_publication pub) static_cast(i) / static_cast(nsit) * 100.0f); }; - l_.info("Building indices for {} situations", nsit); - - auto const buffer_distance = min_distance_; - auto const points_per_circle = 8; - // Note: thomas strategy does not work for geographic_join_round; - // need to use andoyer for that. - using formula = bgeo::strategy::andoyer; - bgeo::strategy::buffer::distance_symmetric distance_strategy{ - buffer_distance - }; - bgeo::strategy::buffer::geographic_join_round join_strategy{ - points_per_circle - }; - bgeo::strategy::buffer::geographic_end_round end_strategy{4}; - bgeo::strategy::buffer::geographic_point_circle circle_strategy{ - points_per_circle - }; - bgeo::strategy::buffer::geographic_side_straight side_strategy; + l_.info("Building indexes for {} situations", nsit); l_.info("Progress: 0%"); auto const before_build = chrono::steady_clock::now(); @@ -135,17 +120,11 @@ handler::handler(log::logger const& l, datex2::situation_publication pub) { for (auto const& ls : rc->relevant_line_strings) { - auto buffered_ls = geo::multi_polygon{}; - bgeo::buffer( - *ls, buffered_ls, distance_strategy, side_strategy, join_strategy, - end_strategy, circle_strategy); - auto blse = geo::box{}; - bgeo::envelope(buffered_ls, blse); - blse_index_.insert(std::make_tuple(blse, buffered_ls, rc)); + blse_index_.insert(*ls, rc); } for (auto p : rc->relevant_points) { - p_index_.insert(std::make_pair(p, rc)); + p_index_.insert(p_index_value{p, rc}); } } @@ -160,7 +139,7 @@ handler::handler(log::logger const& l, datex2::situation_publication pub) auto const after_build = chrono::steady_clock::now(); auto const dur_build_s = chrono::duration{after_build - before_build}; - l_.info("Indices built in {}", dur_build_s); + l_.info("Indexes built in {}", dur_build_s); l_.info("BLSE index size: {}", blse_index_.size()); l_.info("Point index size: {}", p_index_.size()); } @@ -175,16 +154,13 @@ auto handler::process_gpx(gpx::file&& gpx_file) auto const check_periods = time::period_seq{relevant.begin(), relevant.end()}; // TODO: eliminate use of overlap segments - auto splits_with_overlap_segments = std::vector{}; + auto segments = std::vector{}; for (auto const& track : gpx_file.tracks) for (auto const& seg : track.segments) - geo::split_linestring_with_overlap_segments( - seg.waypoints, - 5000 /* meters max total dist until a new split is forced */, - splits_with_overlap_segments); + segments.push_back(seg.waypoints); auto const before_query = chrono::steady_clock::now(); - auto vincenty_strategy = geo::vincenty_strategy{}; + auto vincenty_strategy = geo::wgs84::vincenty_strategy{}; l_.debug("Querying for relevant situations"); auto relevant_road_closures = @@ -192,28 +168,17 @@ auto handler::process_gpx(gpx::file&& gpx_file) auto ls_checked = 0uz; auto p_checked = 0uz; auto i = 0; - for (geo::linestring const& part : splits_with_overlap_segments) + for (auto const& part : segments) { - l_.debug("Checking part [{}/{}]", ++i, splits_with_overlap_segments.size()); + l_.debug("Checking part [{}/{}]", ++i, segments.size()); + auto part_zone_lss = geo::multizonal::zoned_linestring_seg_seq{part}; // TODO: consider buffering with min_distance_ - auto part_box = geo::box{}; + auto part_box = geo::wgs84::box{}; bgeo::envelope(part, part_box); - for (auto it = blse_index_.qbegin(bgeo::index::intersects(part_box)); - it != blse_index_.qend(); it++) + for (auto const& rc : blse_index_.intersection(part_zone_lss)) { - // Cannot use structured bindings here, as boost::geometry::get - // interferes with ADL. It is a candidate as the namespace - // boost::geometry is part of the associated namespace set, which - // happens because geo::linestring ≡ - // boost::geometry::model::linestring is part of the - // whole tuple type (blse_index_value) that is the value_type of the - // iterator. - - geo::box const& blse = std::get<0>(*it); - geo::multi_polygon const& bls = std::get<1>(*it); - std::shared_ptr const& rc = std::get<2>(*it); if (rc->validity && rc->validity->intersect(check_periods).periods().empty()) continue; @@ -221,8 +186,7 @@ auto handler::process_gpx(gpx::file&& gpx_file) // auto small_parts = bgeo::model::multi_linestring{}; // bgeo::intersection(part, lse, small_parts); - if (bgeo::intersects(blse, part) && bgeo::intersects(bls, part)) - relevant_road_closures.emplace(rc); + relevant_road_closures.emplace(rc); // if (!bgeo::is_empty(small_parts)) // if (bgeo::distance(*ls, small_parts, vincenty_strategy) < @@ -234,7 +198,7 @@ auto handler::process_gpx(gpx::file&& gpx_file) it != p_index_.qend(); it++) { // Cannot use structured bindings here for the same reason as above. - geo::point const& p = std::get<0>(*it); + geo::utm::zonable_wgs84_point const& p = std::get<0>(*it); std::shared_ptr const& rc = std::get<1>(*it); if (rc->validity && rc->validity->intersect(check_periods).periods().empty()) @@ -301,10 +265,12 @@ auto handler::process_gpx(gpx::file&& gpx_file) .relevant_lss = rc->relevant_line_strings | views::transform( - [](auto const& lsp) -> geo::linestring + [](auto const& lsp) + -> geo::utm::zonable_wgs84_linestring { return *lsp; }) - | std::ranges:: - to>(), + | std::ranges::to>(), }; }) | std::ranges::to>(), 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; import std; import :datex2; -import :geo; +import :geo.wgs84; +import :geo.multizonal; import :gpx; import :log; import :time; @@ -25,20 +26,20 @@ export namespace routemon::api { struct relevant_road_closure { - std::vector relevant_lss; + std::vector relevant_lss; }; struct relevant_situation { std::string id; - std::optional location; + std::optional location; std::vector comments; std::vector relevant_road_closures; }; struct track_segment { - geo::linestring points; + geo::utm::zonable_wgs84_linestring points; }; struct track @@ -79,13 +80,11 @@ class handler { constexpr static auto const min_distance_ = 5.0; - using blse_index_value = std::tuple< - geo::box, geo::multi_polygon, std::shared_ptr + using blse_index_value = std::shared_ptr; + using blse_index = geo::multizonal::linestring_rtree; + using p_index_value = std::pair< + geo::utm::zonable_wgs84_point, std::shared_ptr >; - using p_index_value = - std::pair>; - using blse_index = - bgeo::index::rtree>; using p_index = bgeo::index::rtree>; 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( } auto const pos_list_str = std::string_view{loc_gml_xml.child_value("loc:posList")}; - // lat1 long1 lat2 long2 ... lat(n-1) long(n-1) latn longn + // lat1 lon1 lat2 lon2 ... lat(n-1) lon(n-1) latn lonn - auto ls = std::make_shared(); + auto ls = std::make_shared(); auto lat_set = false; auto lat = 0.0; - for (auto const lat_or_long_str : std::views::split(pos_list_str, " "sv)) + for (auto const lat_or_lon_str : std::views::split(pos_list_str, " "sv)) { - auto mlat_or_long = util::parse_double(std::string_view{lat_or_long_str}); - if (!mlat_or_long) + auto mlat_or_lon = util::parse_double(std::string_view{lat_or_lon_str}); + if (!mlat_or_lon) { warnings_.insert( std::format( "failed to parse coordinate {:?}", - std::string_view{lat_or_long_str})); + std::string_view{lat_or_lon_str})); return; } if (!lat_set) { - lat = *mlat_or_long; + lat = *mlat_or_lon; lat_set = true; } else { - bgeo::append(*ls, geo::point{*mlat_or_long, lat}); + auto mp = geo::utm::zonable_wgs84_point::from( + geo::wgs84::normalized_point{geo::wgs84::from_lat_lon( + lat, *mlat_or_lon)}); + if (!mp) + warnings_.insert( + std::format( + "point ({}, {}) outside of UTM range; ignoring it", lat, + *mlat_or_lon)); + else + bgeo::append(*ls, *mp); lat = 0; lat_set = false; } @@ -110,11 +119,24 @@ auto loader::add_location_from_xml( // maar heeft het UML-model van DATEX II v3 het over ETRS 89: // https://docs.datex2.eu/_static/data/v3.7/umlmodel/html/EARoot/EA3/EA3/EA5/EA676.htm - auto const coords_etrs89 = geo::point{*mlon, *mlat}; - auto coords_wgs84 = geo::point{}; + auto const coords_etrs89 = geo::wgs84::point{*mlon, *mlat}; + auto coords_wgs84 = geo::wgs84::point{}; etrs89_to_wgs84_.forward(coords_etrs89, coords_wgs84); - rc.relevant_points.push_back(coords_wgs84); + auto mp = geo::utm::zonable_wgs84_point::from( + geo::wgs84::normalized_point{coords_wgs84}); + if (!mp) + { + warnings_.insert( + std::format( + "(ETRS89) point ({}, {}) outside of UTM range after " + "transformation to WGS 84; ignoring it", + *mlat, *mlon)); + } + else + { + rc.relevant_points.push_back(*mp); + } } else { @@ -347,13 +369,23 @@ auto loader::load_situation_publication(std::string const& filename) // het over ETRS 89: // https://docs.datex2.eu/_static/data/v3.7/umlmodel/html/EARoot/EA3/EA3/EA5/EA676.htm - auto const coords_etrs89 = geo::point{*mlon, *mlat}; - auto coords_wgs84 = geo::point{}; + auto const coords_etrs89 = geo::wgs84::point{*mlon, *mlat}; + auto coords_wgs84 = geo::wgs84::point{}; etrs89_to_wgs84_.forward(coords_etrs89, coords_wgs84); - if (!sit->location) + auto mp = geo::utm::zonable_wgs84_point::from( + geo::wgs84::normalized_point{coords_wgs84}); + if (!mp) + { + warnings_.insert( + std::format( + "(ETRS89) point ({}, {}) outside of UTM range after " + "transformation to WGS 84; ignoring it", + *mlat, *mlon)); + } + else if (!sit->location) { - sit->location = coords_wgs84; + sit->location = *mp; } } } 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 @@ module; -#include #include #include @@ -9,7 +8,7 @@ module; export module routemon:datex2; import std; -import :geo; +import :geo.utm; import :time; import :util; @@ -23,14 +22,15 @@ export struct road_closure { std::weak_ptr parent; std::optional validity; - std::vector relevant_points = {}; - std::vector> relevant_line_strings = {}; + std::vector relevant_points = {}; + std::vector> + relevant_line_strings = {}; }; export struct situation { std::string id; - std::optional location = + std::optional location = std::nullopt; // as shown on the map, not used for querying std::vector comments = {}; std::vector> road_closures = {}; @@ -47,7 +47,7 @@ export class loader // ETRS 89 (EPSG:4258) -> WGS 84 (EPSG:4326) bgeo::srs:: transformation, bgeo::srs::static_epsg<4326>> - etrs89_to_wgs84_{}; + etrs89_to_wgs84_ = {}; std::multiset warnings_; 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 @@ -module; - -#include - -module routemon:geo$impl; - -import :geo; - -namespace routemon::geo { - -auto split_linestring_with_overlap_segments( - linestring const& ls, double max_split_distance_m, - std::vector& append_to) -> void -{ - if (bgeo::is_empty(ls)) - return; - - auto current_ls = linestring{}; - auto current_ls_length = 0.0; - auto previous = std::optional{}; - bgeo::for_each_point( - ls, - [&](point p) -> void - { - bgeo::append(current_ls, p); - if (previous) - { - auto d = bgeo::distance(*previous, p, vincenty_strategy()); - current_ls_length += d; - if (current_ls_length > max_split_distance_m) - { - append_to.push_back(std::move(current_ls)); - current_ls = linestring{*previous, p}; - current_ls_length = d; - } - } - previous = p; - }); - append_to.emplace_back(std::move(current_ls)); -} - -} // 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 @@ -module; - -#include - -export module routemon:geo; - -export namespace bgeo = boost::geometry; - -export namespace routemon::geo { - -using point = bgeo::model::point>; -using linestring = bgeo::model::linestring; -using box = bgeo::model::box; -using polygon = bgeo::model::polygon; -using multi_polygon = bgeo::model::multi_polygon; -using stype = bgeo::srs::spheroid; -using vincenty_strategy = bgeo::strategy::distance::vincenty; - -auto split_linestring_with_overlap_segments( - linestring const& ls, double max_split_distance_m, - std::vector& append_to) -> void; - -} // 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 @@ +module; + +#include + +export module routemon:geo; + +namespace 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 @@ +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]) + { + utm::zone_local::prim::box 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++) + { + if (bgeo::intersects(ls, std::get<1>(*it))) + { + co_yield std::get<2>(*it); + } + } + } + } +}; + +} // 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 @@ +module; + +#include +#include + +export module routemon:geo.utm; + +import std; +import :geo; +import :geo.wgs84; + +namespace routemon::geo::utm +{ + +class zonable_wgs84_point +{ + wgs84::normalized_point p_; + + explicit zonable_wgs84_point(wgs84::normalized_point p) + : p_{p} + {} + +public: + zonable_wgs84_point() + : p_{} + { + } + + 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(); + } + + auto lon() const -> double + { + return p_.lon(); + } + + auto lon(double lon) -> void + { + p_.lon(lon); + } + + auto lat(double lat) -> void + { + if (lat < -80 || lat > 84) + throw std::range_error{"latitude out of range for UTM"}; + p_.lat(lat); + } +}; + +} + +namespace boost::geometry::traits { + + namespace { + + using zonable_wgs84_point = routemon::geo::utm::zonable_wgs84_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::wgs84::cs; }; + + 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 + +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). +class zone +{ + // Note: calculations heavily depend on the values of the variants. + enum class hemisphere : std::uint8_t + { + northern = 0, + southern = 1, + }; + + // Negative if in the southern hemisphere + // Equal to (zone_no - 1) * 2 + hemisphere + // In range [0, 119] + std::uint8_t zone_; + + static constexpr auto zone_no_valid(std::uint8_t zone_no) -> bool + { + return 0 < zone_no && zone_no <= 60; + } + + static constexpr auto from(std::uint8_t zone_no, hemisphere h) -> std::uint8_t + { + if (!zone_no_valid(zone_no)) + throw std::invalid_argument{"invalid UTM zone number"}; + return (zone_no - 1) * 2 + static_cast(h); + } + +public: + explicit constexpr zone(std::uint8_t zone_no, hemisphere h) + : zone_{from(zone_no, h)} + {} + + 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; + } + + // Return value in range [min(), max()] + 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}; + } + + auto wgs84_proj_epsg() const -> int + { + switch (hemisphere()) + { + case hemisphere::northern: + return 32600 + zone_no(); + case hemisphere::southern: + return 32700 + zone_no(); + } + } + + // Guarantee: min().as_index() == 0. + static constexpr auto min() noexcept -> zone + { + return zone{1, hemisphere::northern}; + } + + static constexpr auto max() noexcept -> zone + { + return zone{60, hemisphere::southern}; + } + + auto next() const -> zone + { + assert(zone_ <= max().as_index()); + if (zone_ == max().as_index()) + throw std::range_error{"cannot take next of greatest UTM zone"}; + auto copy = zone{*this}; + copy.zone_++; + return copy; + } + + auto operator==(zone rhs) const -> bool + { + return zone_ == rhs.zone_; + } +}; +static_assert(zone::min().as_index() == 0); +static_assert(zone::max().as_index() == 119); + +} // 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 @@ +module; + +#include +#include + +export module routemon:geo.utm.zone_local; + +import :geo; +import :geo.utm; + +// Coordinates that are local to some known UTM zone. +namespace routemon::geo::utm::zone_local { + + using cs = bgeo::cs::cartesian; + + namespace prim { + + using point = bgeo::model::d2::point_xy; + using linestring = bgeo::model::linestring; + using box = bgeo::model::box; + + } // 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 linestring : public prim::linestring + { + zone zone; + + explicit linestring(class zone zone) + : zone{zone} + { + } + }; + +} // namespace routemon::geo::utm::zone_local + +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)); + } + + 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 new file mode 100644 index 0000000..d574d39 --- /dev/null +++ b/server/src/geo/wgs84.cppm @@ -0,0 +1,115 @@ +module; + +#include + +export module routemon:geo.wgs84; + +import :geo; + +namespace routemon::geo::wgs84 { + +using cs = bgeo::cs::geographic; +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 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 +{ + return point{lon, lat}; +} + +auto lat(point p) -> double +{ + return bgeo::get<1>(p); +} + +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 +{ + // Equivalent degrees in the range [0, 360) + auto nonneg_mod360 = [](double t) -> double + { return std::remainder(std::remainder(t, 360.0) + 360.0, 360.0); }; + + auto p_lon = lon(p); + auto p_lat = nonneg_mod360(lat(p)) - 90.0; + if (90.0 < p_lat) + { + assert(p_lat < 270.0); // by nonneg_mod360 + p_lon += 180.0; + p_lat -= 180.0; + } + p_lon = nonneg_mod360(p_lon + 180.0) - 180.0; + + return from_lat_lon(p_lat, p_lon); +} + +auto is_normalized(point p) -> bool +{ + return -90 <= lat(p) && lat(p) <= 90 && -180 <= lon(p) && lon(p) < 180; +} + +class normalized_point +{ + point p_; + +public: + normalized_point(point p) + : p_{is_normalized(p) ? p : normalize(p)} + {} + + normalized_point() + : p_{} + {} + + auto lat() const -> double + { + return bgeo::get<1>(p_); + } + + auto lon() const -> double + { + return bgeo::get<0>(p_); + } + + 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 + { + if (lon < -180 || lon > 180) + throw std::range_error{"longitude out of range"}; + bgeo::set<0>(p_, lon); + } +}; + +class normalized_linestring +{ + std::vector ls_; + +public: + auto empty() const -> bool + { + return ls_.empty(); + } + + auto push_back(normalized_point p) -> void + { + ls_.push_back(p); + } +}; + +} // 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 } auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs) - -> xml::parser + -> xml::parser { auto parse_xml_double = [](std::string_view sv) -> std::optional { return util::parse_double(sv, std::chars_format::fixed); }; @@ -43,8 +43,12 @@ auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs) throw std::runtime_error{ "expected valid latitude and longitude for waypoint" }; + auto mp = geo::utm::zonable_wgs84_point::from( + geo::wgs84::normalized_point{geo::wgs84::from_lat_lon(*mlat, *mlon)}); + if (!mp) + throw std::runtime_error{"expected waypoint to be within UTM range"}; co_await xml::ignore_contents(e); - co_return geo::point{*mlon, *mlat}; + co_return *mp; } auto parse_trkseg(xml::executor_ref e, xml::attribute_view) @@ -121,7 +125,7 @@ auto parse_metadata(xml::executor_ref e, xml::attribute_view) } auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs) - -> xml::parser + -> xml::parser { auto parse_xml_double = [](std::string_view sv) -> std::optional { return util::parse_double(sv, std::chars_format::fixed); }; @@ -142,8 +146,12 @@ auto parse_wpt(xml::executor_ref e, xml::attribute_view attrs) throw std::runtime_error{ "expected valid latitude and longitude for waypoint" }; + auto mp = geo::utm::zonable_wgs84_point::from( + geo::wgs84::normalized_point{geo::wgs84::from_lat_lon(*mlat, *mlon)}); + if (!mp) + throw std::runtime_error{"expected waypoint to be within UTM range"}; co_await xml::ignore_contents(e); - co_return geo::point{*mlon, *mlat}; + co_return *mp; } auto 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 @@ export module routemon:gpx; import std; -import :geo; +import :geo.utm; import :util; import :xml; @@ -15,7 +15,7 @@ struct metadata struct track_segment { - geo::linestring waypoints; + geo::utm::zonable_wgs84_linestring waypoints; }; struct 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 sos << "] " << msg.txt; for (auto const& [k, v] : msg.attrs) sos << " " << k << "=" << std::quoted(v); - sos << '\n'; + sos << std::endl; // we want the flush! } auto make_sink(level lvl) -> std::shared_ptr -- cgit v1.3