Branch data Line data Source code
1 : : #pragma once
2 : :
3 : : /*
4 : : MIT License
5 : :
6 : : Copyright (c) 2014-2024 Stephane Cuillerdier (aka aiekick)
7 : :
8 : : Permission is hereby granted, free of charge, to any person obtaining a copy
9 : : of this software and associated documentation files (the "Software"), to deal
10 : : in the Software without restriction, including without limitation the rights
11 : : to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 : : copies of the Software, and to permit persons to whom the Software is
13 : : furnished to do so, subject to the following conditions:
14 : :
15 : : The above copyright notice and this permission notice shall be included in all
16 : : copies or substantial portions of the Software.
17 : :
18 : : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 : : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 : : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 : : AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 : : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 : : OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24 : : SOFTWARE.
25 : : */
26 : :
27 : : // ezGeo is part of the ezLibs project : https://github.com/aiekick/ezLibs.git
28 : :
29 : : #include <cmath>
30 : : #include <string>
31 : : #include <cctype>
32 : : #include <sstream>
33 : : #include <iomanip>
34 : : #include <utility> // std::pair
35 : :
36 : : #include "../ezMath.hpp"
37 : :
38 : : namespace ez {
39 : : namespace geo {
40 : :
41 : : // WGS84 spherical earth radius in meters
42 : : constexpr double EARTH_RADIUS = 6378137.0;
43 : :
44 : : // Convert Mercator coordinates (x,y in meters) -> (lon,lat) in degrees
45 : 3 : inline ez::dvec2 fromMercatorMetersToDegrees(const ez::dvec2& vPointMeters) {
46 : 3 : ez::dvec2 ret;
47 : 3 : double lonRad = vPointMeters.x / EARTH_RADIUS;
48 : 3 : double latRad = 2.0 * std::atan(std::exp(vPointMeters.y / EARTH_RADIUS)) - M_PI / 2.0;
49 : 3 : ret.x = lonRad * 180.0 / M_PI;
50 : 3 : ret.y = latRad * 180.0 / M_PI;
51 : 3 : return ret;
52 : 3 : }
53 : :
54 : : // Convert (lon,lat) in degrees -> Mercator coordinates (x,y) in meters
55 : 10 : inline ez::dvec2 fromDegressToMercatorMeters(const ez::dvec2& vPointDegree) {
56 : 10 : ez::dvec2 ret{vPointDegree};
57 : : // Clamp latitude to avoid infinities at the poles
58 : 10 : const double maxLat = 89.9999;
59 [ - + ]: 10 : if (ret.y > maxLat)
60 : 0 : ret.y = maxLat;
61 [ - + ]: 10 : if (ret.y < -maxLat)
62 : 0 : ret.y = -maxLat;
63 : :
64 : 10 : double lonRad = ret.x * M_PI / 180.0;
65 : 10 : double latRad = ret.y * M_PI / 180.0;
66 : 10 : ret.x = EARTH_RADIUS * lonRad;
67 : 10 : ret.y = EARTH_RADIUS * std::log(std::tan(M_PI / 4.0 + latRad / 2.0));
68 : 10 : return ret;
69 : 10 : }
70 : :
71 : : // Convert Mercator Y (meters) -> unrolled Y (linearized meters)
72 : 7 : inline double mercatorYToUnrolledY(double yMerc) {
73 : : // 1) Convert Mercator Y -> latitude φ in radians
74 : 7 : double phi = 2.0 * std::atan(std::exp(yMerc / EARTH_RADIUS)) - M_PI / 2.0;
75 : : // 2) Unroll: linearized Y = R * φ
76 : 7 : return EARTH_RADIUS * phi;
77 : 7 : }
78 : :
79 : : // Rescale latitude axis for WGS84 approximated ellipse
80 : 2 : inline double unscaleWGS84Yaxis(double vLatDeg, double vRefDeg = 45.0) {
81 : 2 : const double phi0 = vRefDeg * M_PI / 180.0;
82 : 2 : const double k = 111132.0 / (111320.0 * std::cos(phi0)); // ≈0.9983 / cos(phi0)
83 : 2 : return vLatDeg * k;
84 : 2 : }
85 : :
86 : : /*
87 : : * Parse DEM coordinate component like "S01", "N10", "E001", "W112".
88 : : * Returns the signed numeric value.
89 : : * Example:
90 : : * "N10" -> +10
91 : : * "S10" -> -10
92 : : * "E001" -> +1
93 : : * "W112" -> -112
94 : : */
95 : 4 : inline int16_t parseDemCoordinate(const std::string& vStr) {
96 : 4 : const int16_t value = std::stoi(vStr.substr(1));
97 : 4 : const char c = vStr.at(0);
98 [ + + ][ - + ]: 4 : if (c == 'S' || c == 's' || c == 'W' || c == 'w') {
[ - + ][ + + ]
99 : 2 : return -value;
100 : 2 : }
101 : 2 : return value;
102 : 4 : }
103 : :
104 : : /*
105 : : * Check if a DEM filename is valid.
106 : : * Expected format: N00E000 (7 chars total)
107 : : * - N|S followed by 2 digits
108 : : * - E|W followed by 3 digits
109 : : * Returns true if valid, also fills voLat and voLon with longitude/latitude codes.
110 : : */
111 : 5 : inline bool checkDemFileName(const std::string& vFileName, std::string& voLat, std::string& voLon) {
112 [ + + ]: 5 : if (vFileName.size() != 7U)
113 : 2 : return false;
114 : 3 : try {
115 : 3 : bool ret = true;
116 : 3 : const auto& sLat = vFileName.substr(0, 3);
117 : 3 : const char clat = sLat.at(0);
118 [ - + ][ + + ]: 3 : if (clat == 'S' || clat == 's' || clat == 'N' || clat == 'n') {
[ + + ][ - + ]
119 : 2 : const auto& nLat = sLat.substr(1);
120 [ - + ]: 2 : if (!ez::isInteger(nLat)) {
121 : 0 : ret = false;
122 : 0 : }
123 : 2 : } else {
124 : 1 : ret = false;
125 : 1 : }
126 : 3 : const auto& sLon = vFileName.substr(3);
127 : 3 : const char cLon = sLon.at(0);
128 [ + + ][ - + ]: 3 : if (cLon == 'E' || cLon == 'e' || cLon == 'W' || cLon == 'w') {
[ - + ][ + - ]
129 : 3 : const auto& cx_num = sLon.substr(1);
130 [ - + ]: 3 : if (!ez::isInteger(cx_num)) {
131 : 0 : ret = false;
132 : 0 : }
133 : 3 : } else {
134 : 0 : ret = false;
135 : 0 : }
136 [ + + ]: 3 : if (ret) {
137 : 2 : voLat = sLat;
138 : 2 : voLon = sLon;
139 : 2 : }
140 : 3 : return ret;
141 : 3 : } catch (std::exception& ex) {
142 : 0 : LogVarError("Error : %s", ex.what());
143 : 0 : }
144 : 0 : return false;
145 : 3 : }
146 : :
147 : : /*
148 : : * Check if a DEM filename is valid.
149 : : * Expected format: N00E000 (7 chars total)
150 : : * - N|S followed by 2 digits
151 : : * - E|W followed by 3 digits
152 : : * Returns true if valid, also fills voLat and voLon with longitude/latitude codes.
153 : : */
154 : 0 : inline bool checkDemFileName(const std::string& vFileName, int16_t& voLat, int16_t& voLon) {
155 : 0 : std::string lat;
156 : 0 : std::string lon;
157 : 0 : if (checkDemFileName(vFileName, lat, lon)) {
158 : 0 : voLat = parseDemCoordinate(lat);
159 : 0 : voLon = parseDemCoordinate(lon);
160 : 0 : return true;
161 : 0 : }
162 : 0 : return false;
163 : 0 : }
164 : :
165 : : /*
166 : : * DMS class:
167 : : * Stores coordinates in Degrees, Minutes, Seconds and a cardinal letter.
168 : : * Examples of valid inputs: "48°51'24\"N", "45 30 15 S"
169 : : */
170 : : struct dms {
171 : : float deg{};
172 : : float min{};
173 : : float sec{};
174 : : char letter{' '};
175 : : bool valid{false};
176 : :
177 : 36 : dms() = default;
178 : :
179 : : // Construct from a DMS string like "45°30'15\"N" or "45 30 15 N"
180 : 5 : explicit dms(const std::string& vDms) { valid = parse(vDms); }
181 : :
182 : : // Construct from a decimal angle, with flag isLat (true for latitude, false for longitude)
183 : 154 : explicit dms(double vAngle, bool vIsLat) {
184 : 154 : fromAngle(vAngle, vIsLat);
185 : 154 : valid = true;
186 : 154 : }
187 : :
188 : : // Basic setters
189 : 9 : dms& setDeg(float vDeg) {
190 : 9 : deg = vDeg;
191 : 9 : return *this;
192 : 9 : }
193 : 9 : dms& setMin(float vMin) {
194 : 9 : min = vMin;
195 : 9 : return *this;
196 : 9 : }
197 : 9 : dms& setSec(float vSec) {
198 : 9 : sec = vSec;
199 : 9 : return *this;
200 : 9 : }
201 : 10 : dms& setLetter(char vLetter) {
202 : 10 : letter = vLetter;
203 : 10 : return *this;
204 : 10 : }
205 : :
206 : : // offsets
207 : : // Adds degrees, pushes fractional part to minutes, then wraps if lat/lon.
208 : 25 : dms& offsetDeg(float vDeg) {
209 : : // 0) First, purge any existing fractional degrees -> push into minutes
210 : 25 : const float dIntExisting = ez::floor(deg);
211 : 25 : const float dFracExisting = deg - dIntExisting;
212 : 25 : deg = dIntExisting;
213 [ + + ]: 25 : if (dFracExisting != 0.0f) {
214 : : // push to minutes as fractional part
215 : 1 : sec += (dFracExisting * 60.0f) * 60.0f; // deg.frac * 60 = min.frac; *60 again -> seconds
216 : 1 : }
217 : :
218 : : // 1) Split input degrees: integer to deg, fractional cascades to min/sec
219 : 25 : const float dInt = ez::floor(vDeg);
220 : 25 : const float dFrac = vDeg - dInt;
221 : 25 : deg += dInt;
222 : :
223 [ + + ]: 25 : if (dFrac != 0.0f) {
224 : : // dFrac degrees -> minutes.frac = dFrac * 60 -> seconds addition
225 : 2 : sec += (dFrac * 60.0f) * 60.0f;
226 : 2 : }
227 : :
228 : : // 2) Normalize seconds -> minutes
229 : 25 : const float carryMin = ez::floor(sec / 60.0f);
230 : 25 : sec -= carryMin * 60.0f;
231 : 25 : min += carryMin;
232 : :
233 : : // 3) Ensure 'min' is integer: push fractional residue into seconds
234 : 25 : const float minFrac = min - ez::floor(min);
235 [ - + ]: 25 : if (minFrac != 0.0f) {
236 : 0 : min = ez::floor(min);
237 : 0 : sec += minFrac * 60.0f;
238 : :
239 : : // Re-normalize seconds
240 : 0 : const float carryMin2 = ez::floor(sec / 60.0f);
241 : 0 : sec -= carryMin2 * 60.0f;
242 : 0 : min += carryMin2;
243 : 0 : }
244 : :
245 : : // 4) Carry minutes into degrees
246 : 25 : const float carryDeg = ez::floor(min / 60.0f);
247 : 25 : min -= carryDeg * 60.0f;
248 : 25 : deg += carryDeg;
249 : :
250 : : // 5) Optional wrap
251 [ + + ]: 25 : if (isLat()) {
252 : 14 : deg = ez::mod(deg, 90.0f);
253 [ + - ]: 14 : } else if (isLon()) {
254 : 11 : deg = ez::mod(deg, 180.0f);
255 : 11 : }
256 : 25 : return *this;
257 : 25 : }
258 : :
259 : : // Adds minutes, normalizes into [0,60), carries into degrees, wraps if needed.
260 : 3 : dms& offsetMin(float vMin) {
261 : : // 1) Split input minutes: integer part goes to minutes, fractional to seconds
262 : 3 : const float mInt = ez::floor(vMin);
263 : 3 : const float mFrac = vMin - mInt;
264 : :
265 : : // 2) Add integer minutes
266 : 3 : min += mInt;
267 : :
268 : : // 3) Push fractional minutes into seconds
269 [ + - ]: 3 : if (mFrac != 0.0f) {
270 : 3 : sec += mFrac * 60.0f;
271 : 3 : }
272 : :
273 : : // 4) Normalize seconds -> minutes
274 : 3 : const float carryMin = ez::floor(sec / 60.0f);
275 : 3 : sec -= carryMin * 60.0f;
276 : 3 : min += carryMin;
277 : :
278 : : // 5) Ensure 'min' is integer: push any fractional residue into seconds
279 : 3 : const float minFrac = min - ez::floor(min);
280 [ - + ]: 3 : if (minFrac != 0.0f) {
281 : 0 : min = ez::floor(min);
282 : 0 : sec += minFrac * 60.0f;
283 : :
284 : : // Re-normalize seconds since we just added some
285 : 0 : const float carryMin2 = ez::floor(sec / 60.0f);
286 : 0 : sec -= carryMin2 * 60.0f;
287 : 0 : min += carryMin2;
288 : 0 : }
289 : :
290 : : // 6) Carry minutes into degrees (keep min in [0,60) and integer)
291 : 3 : const float carryDeg = ez::floor(min / 60.0f);
292 : 3 : min -= carryDeg * 60.0f;
293 : 3 : deg += carryDeg;
294 : :
295 : : // 7) Optional wrap
296 [ + - ]: 3 : if (isLat()) {
297 : 3 : deg = ez::mod(deg, 90.0f);
298 [ # # ]: 3 : } else if (isLon()) {
299 : 0 : deg = ez::mod(deg, 180.0f);
300 : 0 : }
301 : 3 : return *this;
302 : 3 : }
303 : :
304 : : // Adds seconds, normalizes into [0,60), carries into minutes (then degrees via offsetMin).
305 : 3 : dms& offsetSec(float vSec) {
306 : : // 1) Add seconds
307 : 3 : sec += vSec;
308 : :
309 : : // 2) Normalize seconds into [0,60) and carry into minutes (floor works for negatives)
310 : 3 : const float carryMin = ez::floor(sec / 60.0f);
311 : 3 : sec -= carryMin * 60.0f;
312 : 3 : min += carryMin;
313 : :
314 : : // 3) Ensure 'min' is an integer: push any fractional minutes into seconds
315 : 3 : const float minFrac = min - ez::floor(min);
316 [ - + ]: 3 : if (minFrac != 0.0f) {
317 : 0 : min = ez::floor(min);
318 : 0 : sec += minFrac * 60.0f;
319 : :
320 : : // Re-normalize seconds because we just added some
321 : 0 : const float carryMin2 = ez::floor(sec / 60.0f);
322 : 0 : sec -= carryMin2 * 60.0f;
323 : 0 : min += carryMin2;
324 : 0 : }
325 : :
326 : : // 4) Carry minutes into degrees, keep min in [0,60) and integer
327 : 3 : const float carryDeg = ez::floor(min / 60.0f);
328 : 3 : min -= carryDeg * 60.0f;
329 : 3 : deg += carryDeg;
330 : :
331 : : // 5) Optional wrap on degrees (no hemisphere flip)
332 [ + - ]: 3 : if (isLat()) {
333 : 3 : deg = ez::mod(deg, 90.0f);
334 [ # # ]: 3 : } else if (isLon()) {
335 : 0 : deg = ez::mod(deg, 180.0f);
336 : 0 : }
337 : 3 : return *this;
338 : 3 : }
339 : :
340 : : // Type checks
341 [ + + ][ - + ]: 33 : bool isLat() const { return (letter == 'N') || (letter == 'S'); }
342 [ + + ][ - + ]: 13 : bool isLon() const { return (letter == 'E') || (letter == 'W'); }
343 : 0 : bool isValid() const { return valid; }
344 : :
345 : : // Normalize DMS values by carrying overflow of seconds and minutes
346 : 7 : static void normalizeDMS(double& D, double& M, double& S) {
347 [ - + ]: 7 : if (S >= 59.9999995) {
348 : 0 : S = 0.0;
349 : 0 : M += 1.0;
350 : 0 : }
351 [ - + ]: 7 : if (M >= 60.0) {
352 : 0 : M = 0.0;
353 : 0 : D += 1.0;
354 : 0 : }
355 : 7 : }
356 : :
357 : : // Parse a DMS string into components
358 : 12 : bool parse(const std::string& vDms) {
359 : 12 : std::string str;
360 [ + + ]: 88 : for (char c : vDms) {
361 : 88 : unsigned char uc = static_cast<unsigned char>(c);
362 [ + + ][ + + ]: 88 : if (std::isdigit(uc) || uc == '.' || uc == '-' || uc == '+' || uc == ' ') {
[ + + ][ - + ]
[ + + ]
363 : 76 : str += static_cast<char>(uc);
364 [ + + ]: 76 : } else if (std::isalpha(uc)) {
365 : 5 : letter = static_cast<char>(std::toupper(uc));
366 : 7 : } else {
367 : 7 : str += ' ';
368 : 7 : }
369 : 88 : }
370 : 12 : std::istringstream iss(str);
371 : 12 : std::vector<double> vals;
372 [ + + ]: 32 : for (double x; iss >> x;) {
373 : 28 : vals.push_back(x);
374 [ + + ]: 28 : if (vals.size() == 3) {
375 : 8 : break;
376 : 8 : }
377 : 28 : }
378 [ + + ]: 12 : if (vals.empty()) {
379 : 1 : return false;
380 : 1 : }
381 : 11 : double d = 0.0, m = 0.0, s = 0.0;
382 [ + + ]: 11 : const double sign = (vals[0] < 0.0) ? -1.0 : 1.0;
383 [ + + ][ + + ]: 11 : if (letter == ' ' && sign < 0.0){
384 : 1 : letter = '-'; // special case for parse("-12.5")
385 : 1 : }
386 [ + + ]: 11 : if (vals.size() == 1) {
387 : : // Decimal degrees -> convert to DMS
388 : 2 : const double ad = std::fabs(vals[0]);
389 : 2 : d = std::floor(ad + 1e-12);
390 : 2 : const double rem_min = (ad - d) * 60.0;
391 : 2 : m = std::floor(rem_min + 1e-12);
392 : 2 : s = (rem_min - m) * 60.0;
393 : 2 : normalizeDMS(d, m, s);
394 [ + + ]: 9 : } else if (vals.size() == 2) {
395 [ - + ]: 1 : if (vals[1] < 0.0) {
396 : 0 : return false;
397 : 0 : }
398 : : // Degrees + minutes
399 : 1 : d = std::floor(std::fabs(vals[0]) + 1e-12);
400 : 1 : m = std::fabs(vals[1]);
401 [ + - ][ + - ]: 1 : if (!(m >= 0.0 && m < 60.0)) {
402 : 0 : return false;
403 : 0 : }
404 : 1 : const double mi = std::floor(m + 1e-12);
405 : 1 : s = (m - mi) * 60.0;
406 : 1 : m = mi;
407 : 1 : normalizeDMS(d, m, s);
408 : 8 : } else {
409 [ + + ]: 8 : if (vals[1] < 0.0) {
410 : 1 : return false;
411 : 1 : }
412 [ + + ]: 7 : if (vals[2] < 0.0) {
413 : 1 : return false;
414 : 1 : }
415 : : // Degrees + minutes + seconds
416 : 6 : d = std::floor(std::fabs(vals[0]) + 1e-12);
417 : 6 : m = std::fabs(vals[1]);
418 : 6 : s = std::fabs(vals[2]);
419 [ + - ][ + + ]: 6 : if (!(m >= 0.0 && m < 60.0)) {
420 : 1 : return false;
421 : 1 : }
422 [ + - ][ + + ]: 5 : if (!(s >= 0.0 && s < 60.0)) {
423 : 1 : return false;
424 : 1 : }
425 : 4 : normalizeDMS(d, m, s);
426 : 4 : }
427 : 7 : deg = static_cast<float>(d);
428 : 7 : min = static_cast<float>(m);
429 : 7 : sec = static_cast<float>(s);
430 : 7 : return true;
431 : 11 : }
432 : :
433 : : // Convert DMS -> decimal angle
434 : 206 : double toAngle() const {
435 : 206 : return getSign() * (deg + min / 60.0 + sec / 3600.0);
436 : 206 : }
437 : :
438 : : // get the sign
439 [ + + ][ - + ]: 206 : double getSign() const {return (letter == 'S' || letter == 'W' || letter == '-') ? -1.0 : 1.0;}
[ + + ]
440 : :
441 : : // Convert DMS -> string
442 : 0 : std::string toDmsString() const {
443 : 0 : std::stringstream ss;
444 : 0 : ss << deg << "°" << min << "'" << sec << "\"" << letter;
445 : 0 : return ss.str();
446 : 0 : }
447 : :
448 : : // Convert decimal angle -> DMS + letter
449 : 154 : void fromAngle(const double vAngle, const bool vIsLat) {
450 [ + + ]: 154 : const double sign = (vAngle < 0.0) ? -1.0 : 1.0;
451 : 154 : const double a = std::fabs(vAngle);
452 : 154 : deg = static_cast<float>(std::floor(a));
453 : 154 : const double frac = (a - deg) * 60.0;
454 : 154 : min = static_cast<float>(std::floor(frac));
455 : 154 : sec = static_cast<float>((frac - min) * 60.0);
456 [ + + ]: 154 : if (vIsLat) {
457 [ + + ]: 77 : letter = (sign < 0.0) ? 'S' : 'N';
458 : 77 : } else {
459 [ - + ]: 77 : letter = (sign < 0.0) ? 'W' : 'E';
460 : 77 : }
461 : 154 : }
462 : : };
463 : :
464 : : } // namespace geo
465 : : } // namespace ez
|
