Some cleanups

9 files changed, 494 insertions, 457 deletions

diff --git a/server/CMakeLists.txt b/server/CMakeLists.txt
index 53cda76..fad973f 100644
--- a/server/CMakeLists.txt
+++ b/server/CMakeLists.txt
@@ -64,8 +64,8 @@ target_sources(routemon_lib
     src/xml.cppm
     src/xml.cpp
     src/geo/geo.cppm
-    src/geo/multizonal.cppm
     src/geo/utm.cppm
+    src/geo/utm_multizonal.cppm
     src/geo/utm_zone_local.cppm
     src/geo/wgs84.cppm
 )
diff --git a/server/cf.sh b/server/cf.sh
index b8bd392..5c5ec3a 100755
--- a/server/cf.sh
+++ b/server/cf.sh
@@ -1,3 +1,5 @@
 #!/bin/bash
 
-clang-format -i src/*.cpp src/*.cppm
+shopt -s globstar
+
+clang-format -i src/**/*.cpp src/**/*.cppm
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<datex2::road_closure>;
-  using blse_index = geo::multizonal::linestring_rtree<blse_index_value>;
+  using blse_index = geo::utm::multizonal::linestring_rtree<blse_index_value>;
   using bpe_index_value = std::shared_ptr<datex2::road_closure>;
-  using bpe_index = geo::multizonal::point_rtree<bpe_index_value>;
+  using bpe_index = geo::utm::multizonal::point_rtree<bpe_index_value>;
 
   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 <boost/geometry.hpp>
-#include <boost/geometry/srs/epsg.hpp>
-#include <boost/geometry/srs/transformation.hpp>
-
-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<std::default_initializable T>
-class multi_zone
-{
-  std::array<T, utm::zone::max().as_index() + 1> 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<bgeo::srs::static_epsg<4326>>;
-
-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_transform>
-{
-  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<utm::zone, double>
-{
-  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<std::vector<utm::zone_local::prim::linestring>> 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<utm::zone_local::prim::linestring>{};
-    auto mprev_p_zone = std::optional<utm::zone>{};
-    auto mprev_p_alt_zone = std::optional<utm::zone>{};
-
-    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<utm::zone>{};
-      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<utm::zone_local::linestring> segments;
-
-  explicit zoned_linestring_seg_seq(utm::zonable_wgs84_linestring ls)
-  {
-    auto to_utm = wgs84::utm_transforms::instance();
-    auto mworking_seg = std::optional<utm::zone_local::linestring>{};
-
-    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 T>
-class linestring_rtree
-{
-public:
-  using index_value = std::tuple<utm::zone_local::prim::box, utm::zone_local::prim::linestring, T>;
-
-private:
-  multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>> local_rtrees_;
-
-public:
-  template<std::convertible_to<T> 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<utm::zone_local::prim::box>(bgeo::return_envelope<utm::zone_local::prim::box>(ls), 30.0 /* m */);
-        local_rtrees_[z].insert(index_value{blse, std::move(ls), std::forward<U>(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<T const&>
-  {
-    for (auto const& ls : lss.segments)
-    {
-      // auto ls_box = bgeo::return_envelope<utm::zone_local::prim::box>(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<utm::zone_local::prim::linestring>;
-
-        auto relevant_ls_parts = multi_linestring{};
-        bgeo::intersection(std::get<0>(*it), static_cast<utm::zone_local::prim::linestring const&>(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 T>
-class point_rtree
-{
-public:
-  using index_value = std::tuple<utm::zone_local::prim::box, utm::zone_local::prim::point, T>;
-
-private:
-  multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>> local_rtrees_;
-
-public:
-  template<std::convertible_to<T> 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<utm::zone>{};
-    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<utm::zone_local::prim::box>(bgeo::return_envelope<utm::zone_local::prim::box>(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<T const&>
-  {
-    for (auto const& ls : lss.segments)
-    {
-      // auto ls_box = bgeo::return_envelope<utm::zone_local::prim::box>(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<utm::zone_local::prim::linestring>;
-
-        auto relevant_ls_parts = multi_linestring{};
-        bgeo::intersection(std::get<0>(*it), static_cast<utm::zone_local::prim::linestring const&>(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 <cassert>
 #include <boost/geometry/geometries/geometries.hpp>
+#include <cassert>
 
 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<zonable_wgs84_point>
+  inline 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();
-  }
+  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 <anonymous>
+} // namespace
 
-  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 <>
+struct tag<zonable_wgs84_point>
+{
+  using type = point_tag;
+};
 
-  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");
-    }
+template <>
+struct dimension<zonable_wgs84_point> : 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<zonable_wgs84_point>
+{
+  using type = double;
+};
 
-} // namespace boost::geometry::traits
+template <>
+struct coordinate_system<zonable_wgs84_point>
+{
+  using type = routemon::geo::wgs84::cs;
+};
 
-namespace routemon::geo::utm
+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
 
-using zonable_wgs84_linestring = bgeo::model::linestring<zonable_wgs84_point, std::vector>;
+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).
@@ -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<enum hemisphere>(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<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};
+    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 <boost/geometry.hpp>
+#include <boost/geometry/srs/epsg.hpp>
+#include <boost/geometry/srs/transformation.hpp>
+
+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 <std::default_initializable T>
+class multi_zone
+{
+  std::array<T, zone::max().as_index() + 1> 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<bgeo::srs::static_epsg<4326>>;
+
+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_transform>
+{
+  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<utm::zone, double>
+{
+  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<std::vector<utm::zone_local::prim::linestring>>
+{
+  auto result = multi_zone<std::vector<utm::zone_local::prim::linestring>>{};
+
+  auto to_utm = from_wgs84_transforms::instance();
+  auto working = multi_zone<utm::zone_local::prim::linestring>{};
+  auto mprev_p_zone = std::optional<utm::zone>{};
+  auto mprev_p_alt_zone = std::optional<utm::zone>{};
+
+  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<utm::zone>{};
+    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<utm::zone_local::linestring>
+{
+  auto splits = std::vector<utm::zone_local::linestring>{};
+
+  auto to_utm = from_wgs84_transforms::instance();
+  auto mworking_seg = std::optional<utm::zone_local::linestring>{};
+
+  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 T>
+class linestring_rtree
+{
+public:
+  using index_value = std::
+      tuple<utm::zone_local::prim::box, utm::zone_local::prim::linestring, T>;
+
+private:
+  multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>>
+      local_rtrees_;
+
+public:
+  template <std::convertible_to<T> 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<utm::zone_local::prim::box>(
+            bgeo::return_envelope<utm::zone_local::prim::box>(ls),
+            30.0 /* m */);
+        local_rtrees_[z].insert(
+            index_value{blse, std::move(ls), std::forward<U>(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<utm::zone_local::linestring> const& lss) const
+      -> std::generator<T const&>
+  {
+    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<utm::zone_local::prim::linestring>;
+
+        auto relevant_ls_parts = multi_linestring{};
+        bgeo::intersection(
+            std::get<0>(*it),
+            static_cast<utm::zone_local::prim::linestring const&>(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 T>
+class point_rtree
+{
+public:
+  using index_value =
+      std::tuple<utm::zone_local::prim::box, utm::zone_local::prim::point, T>;
+
+private:
+  multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>>
+      local_rtrees_;
+
+public:
+  template <std::convertible_to<T> 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<utm::zone>{};
+    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<utm::zone_local::prim::box>(
+          bgeo::return_envelope<utm::zone_local::prim::box>(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<utm::zone_local::linestring> const& lss) const
+      -> std::generator<T const&>
+  {
+    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<utm::zone_local::prim::linestring>;
+
+        auto relevant_ls_parts = multi_linestring{};
+        bgeo::intersection(
+            std::get<0>(*it),
+            static_cast<utm::zone_local::prim::linestring const&>(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<double>;
-    using linestring = bgeo::model::linestring<point>;
-    using box = bgeo::model::box<point>;
+using point = bgeo::model::d2::point_xy<double>;
+using linestring = bgeo::model::linestring<point>;
+using box = bgeo::model::box<point>;
 
-  } // 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 {
+namespace {
 
-    using zone_local_point = routemon::geo::utm::zone_local::point;
+using zone_local_point = routemon::geo::utm::zone_local::point;
 
-  } // namespace <anonymous>
+} // namespace
 
-  template<> struct tag<zone_local_point> { using type = point_tag; };
-  template<> struct dimension<zone_local_point> : boost::mpl::int_<2> {};
-  template<> struct coordinate_type<zone_local_point> { using type = double; };
-  template<> struct coordinate_system<zone_local_point> { using type = routemon::geo::utm::zone_local::cs; };
+template <>
+struct tag<zone_local_point>
+{
+  using type = point_tag;
+};
+template <>
+struct dimension<zone_local_point> : boost::mpl::int_<2>
+{
+};
+template <>
+struct coordinate_type<zone_local_point>
+{
+  using type = double;
+};
+template <>
+struct coordinate_system<zone_local_point>
+{
+  using type = routemon::geo::utm::zone_local::cs;
+};
 
-  template<std::size_t index>
-  struct access<zone_local_point, index> {
-    static inline auto get(zone_local_point const& p) -> double
-    {
-      return bgeo::get<index>(static_cast<routemon::geo::utm::zone_local::prim::point const&>(p));
-    }
+template <std::size_t index>
+struct access<zone_local_point, index>
+{
+  static inline auto get(zone_local_point const& p) -> double
+  {
+    return bgeo::get<index>(
+        static_cast<routemon::geo::utm::zone_local::prim::point const&>(p));
+  }
 
-    static inline auto set(zone_local_point& p, double v) -> void
-    {
-      return bgeo::set<index>(static_cast<routemon::geo::utm::zone_local::prim::point&>(p), v);
-    }
-  };
+  static inline auto set(zone_local_point& p, double v) -> void
+  {
+    return bgeo::set<index>(
+        static_cast<routemon::geo::utm::zone_local::prim::point&>(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<double, 2, cs>;
 // 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<point, std::vector /* the default */>;
+using linestring =
+    bgeo::model::linestring<point, std::vector /* the default */>;
 using box = bgeo::model::box<point>;
 using stype = bgeo::srs::spheroid<double>;
 using vincenty_strategy = bgeo::strategy::distance::vincenty<stype>;
 
-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<normalized_point> 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