summaryrefslogtreecommitdiffstats
path: root/server/src/geo/utm.cppm
blob: c58ffb465d6ffbb5f9de8b944dbc9b73a682562f (about) (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
module;

#include <boost/geometry/geometries/geometries.hpp>
#include <cassert>

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 inline zonable_wgs84_point(wgs84::normalized_point p) : p_{p} {}

public:
  inline zonable_wgs84_point() : p_{} {}

  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};
  }

  inline auto lat() const -> double { return p_.lat(); }

  inline auto lon() const -> double { return p_.lon(); }

  inline auto lon(double lon) -> void { p_.lon(lon); }

  inline auto lat(double lat) -> void
  {
    if (lat < -80 || lat > 84)
      throw std::range_error{"latitude out of range for UTM"};
    p_.lat(lat);
  }
};

} // namespace routemon::geo::utm

namespace boost::geometry::traits {

namespace {

using zonable_wgs84_point = routemon::geo::utm::zonable_wgs84_point;

} // 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 <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)}
  {
  }

  constexpr auto hemisphere() const -> enum hemisphere
  {
    return static_cast<enum hemisphere>(zone_ % 2);
  }

  // Return value in range [1, 60]
  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_; }

  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
    };
  }

  constexpr 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};
  }

  constexpr 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;
  }

  constexpr 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