1 files changed, 199 insertions, 0 deletions

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 <cassert>
+#include <boost/geometry/geometries/geometries.hpp>
+
+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<zonable_wgs84_point>
+  {
+    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 <anonymous>
+
+  template<> struct tag<zonable_wgs84_point> { using type = point_tag; };
+  template<> struct dimension<zonable_wgs84_point> : boost::mpl::int_<2> {};
+  template<> struct coordinate_type<zonable_wgs84_point> { using type = double; };
+  template<> struct coordinate_system<zonable_wgs84_point> { using type = routemon::geo::wgs84::cs; };
+
+  template<std::size_t index>
+  struct access<zonable_wgs84_point, index> {
+    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<zonable_wgs84_point, std::vector>;
+
+// 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<std::uint8_t>(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<enum hemisphere>(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<std::uint8_t>(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