File: | /home/gilles/devel/GIT/8.x/core/utilities/geolocation/engine/geodata/data/GeoDataCoordinates.cpp |
Warning: | line 706, column 22 Although the value stored to 'lonMinF' is used in the enclosing expression, the value is never actually read from 'lonMinF' |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
1 | /* ============================================================ |
2 | * |
3 | * This file is a part of digiKam project |
4 | * https://www.digikam.org |
5 | * |
6 | * Date : 2023-05-15 |
7 | * Description : geolocation engine based on Marble. |
8 | * (c) 2007-2022 Marble Team |
9 | * https://invent.kde.org/education/marble/-/raw/master/data/credits_authors.html |
10 | * |
11 | * SPDX-FileCopyrightText: 2023-2024 by Gilles Caulier <caulier dot gilles at gmail dot com> |
12 | * |
13 | * SPDX-License-Identifier: LGPL-2.1-or-later |
14 | * |
15 | * ============================================================ */ |
16 | |
17 | #include "GeoDataCoordinates.h" |
18 | #include "GeoDataCoordinates_p.h" |
19 | |
20 | // Qt includes |
21 | |
22 | #include <qmath.h> |
23 | #include <QDataStream> |
24 | #include <QPointF> |
25 | |
26 | // KDE includes |
27 | |
28 | #include <klocalizedstring.h> |
29 | |
30 | // Local includes |
31 | |
32 | #include "LonLatParser_p.h" |
33 | #include "MarbleMath.h" |
34 | #include "Quaternion.h" |
35 | #include "digikam_debug.h" |
36 | |
37 | namespace Marble |
38 | { |
39 | |
40 | const qreal GeoDataCoordinatesPrivate::sm_semiMajorAxis = 6378137.0; |
41 | const qreal GeoDataCoordinatesPrivate::sm_semiMinorAxis = 6356752.314; |
42 | const qreal GeoDataCoordinatesPrivate::sm_eccentricitySquared = 6.69437999013e-03; |
43 | const qreal GeoDataCoordinatesPrivate::sm_utmScaleFactor = 0.9996; |
44 | GeoDataCoordinates::Notation GeoDataCoordinates::s_notation = GeoDataCoordinates::DMS; |
45 | |
46 | const GeoDataCoordinates GeoDataCoordinates::null = GeoDataCoordinates(0, 0, 0); // don't use default constructor! |
47 | |
48 | GeoDataCoordinates::GeoDataCoordinates(qreal _lon, qreal _lat, qreal _alt, GeoDataCoordinates::Unit unit, int _detail) |
49 | : d(new GeoDataCoordinatesPrivate(_lon, _lat, _alt, unit, _detail)) |
50 | { |
51 | d->ref.ref(); |
52 | } |
53 | |
54 | /* simply copy the d pointer |
55 | * it will be replaced in the detach function instead |
56 | */ |
57 | GeoDataCoordinates::GeoDataCoordinates(const GeoDataCoordinates& other) |
58 | : d(other.d) |
59 | { |
60 | d->ref.ref(); |
61 | } |
62 | |
63 | /* simply copy null's d pointer |
64 | * it will be replaced in the detach function |
65 | */ |
66 | GeoDataCoordinates::GeoDataCoordinates() |
67 | : d(null.d) |
68 | { |
69 | d->ref.ref(); |
70 | } |
71 | |
72 | /* |
73 | * only delete the private d pointer if the number of references is 0 |
74 | * remember that all copies share the same d pointer! |
75 | */ |
76 | GeoDataCoordinates::~GeoDataCoordinates() |
77 | { |
78 | delete d->m_q; |
79 | d->m_q = nullptr; |
80 | |
81 | if (!d->ref.deref()) |
82 | { |
83 | delete d; |
84 | } |
85 | |
86 | #ifdef DEBUG_GEODATA |
87 | // qCDebug(DIGIKAM_MARBLE_LOG) << "delete coordinates"; |
88 | #endif |
89 | } |
90 | |
91 | bool GeoDataCoordinates::isValid() const |
92 | { |
93 | return d != null.d; |
94 | } |
95 | |
96 | /* |
97 | * if only one copy exists, return |
98 | * else make a new private d pointer object and assign the values of the current |
99 | * one to it |
100 | * at the end, if the number of references thus reaches 0 delete it |
101 | * this state shouldn't happen, but if it does, we have to clean up behind us. |
102 | */ |
103 | void GeoDataCoordinates::detach() |
104 | { |
105 | delete d->m_q; |
106 | d->m_q = nullptr; |
107 | |
108 | if (d->ref.fetchAndAddRelaxed(0) == 1) |
109 | { |
110 | return; |
111 | } |
112 | |
113 | GeoDataCoordinatesPrivate* new_d = new GeoDataCoordinatesPrivate(*d); |
114 | |
115 | if (!d->ref.deref()) |
116 | { |
117 | delete d; |
118 | } |
119 | |
120 | d = new_d; |
121 | d->ref.ref(); |
122 | } |
123 | |
124 | /* |
125 | * call detach() at the start of all non-static, non-const functions |
126 | */ |
127 | void GeoDataCoordinates::set(qreal _lon, qreal _lat, qreal _alt, GeoDataCoordinates::Unit unit) |
128 | { |
129 | detach(); |
130 | d->m_altitude = _alt; |
131 | |
132 | switch (unit) |
133 | { |
134 | default: |
135 | case Radian: |
136 | d->m_lon = _lon; |
137 | d->m_lat = _lat; |
138 | break; |
139 | |
140 | case Degree: |
141 | d->m_lon = _lon * DEG2RAD; |
142 | d->m_lat = _lat * DEG2RAD; |
143 | break; |
144 | } |
145 | } |
146 | |
147 | /* |
148 | * call detach() at the start of all non-static, non-const functions |
149 | */ |
150 | void GeoDataCoordinates::setLongitude(qreal _lon, GeoDataCoordinates::Unit unit) |
151 | { |
152 | detach(); |
153 | |
154 | switch (unit) |
155 | { |
156 | default: |
157 | case Radian: |
158 | d->m_lon = _lon; |
159 | break; |
160 | |
161 | case Degree: |
162 | d->m_lon = _lon * DEG2RAD; |
163 | break; |
164 | } |
165 | } |
166 | |
167 | |
168 | /* |
169 | * call detach() at the start of all non-static, non-const functions |
170 | */ |
171 | void GeoDataCoordinates::setLatitude(qreal _lat, GeoDataCoordinates::Unit unit) |
172 | { |
173 | detach(); |
174 | |
175 | switch (unit) |
176 | { |
177 | case Radian: |
178 | d->m_lat = _lat; |
179 | break; |
180 | |
181 | case Degree: |
182 | d->m_lat = _lat * DEG2RAD; |
183 | break; |
184 | } |
185 | } |
186 | |
187 | void GeoDataCoordinates::geoCoordinates(qreal& lon, qreal& lat, |
188 | GeoDataCoordinates::Unit unit) const |
189 | { |
190 | switch (unit) |
191 | { |
192 | default: |
193 | case Radian: |
194 | lon = d->m_lon; |
195 | lat = d->m_lat; |
196 | break; |
197 | |
198 | case Degree: |
199 | lon = d->m_lon * RAD2DEG; |
200 | lat = d->m_lat * RAD2DEG; |
201 | break; |
202 | } |
203 | } |
204 | |
205 | void GeoDataCoordinates::geoCoordinates(qreal& lon, qreal& lat) const |
206 | { |
207 | lon = d->m_lon; |
208 | lat = d->m_lat; |
209 | } |
210 | |
211 | void GeoDataCoordinates::geoCoordinates(qreal& lon, qreal& lat, qreal& alt, |
212 | GeoDataCoordinates::Unit unit) const |
213 | { |
214 | geoCoordinates(lon, lat, unit); |
215 | alt = d->m_altitude; |
216 | } |
217 | |
218 | void GeoDataCoordinates::geoCoordinates(qreal& lon, qreal& lat, qreal& alt) const |
219 | { |
220 | lon = d->m_lon; |
221 | lat = d->m_lat; |
222 | alt = d->m_altitude; |
223 | } |
224 | |
225 | qreal GeoDataCoordinates::longitude(GeoDataCoordinates::Unit unit) const |
226 | { |
227 | switch (unit) |
228 | { |
229 | default: |
230 | case Radian: |
231 | return d->m_lon; |
232 | |
233 | case Degree: |
234 | return d->m_lon * RAD2DEG; |
235 | } |
236 | } |
237 | |
238 | qreal GeoDataCoordinates::longitude() const |
239 | { |
240 | return d->m_lon; |
241 | } |
242 | |
243 | qreal GeoDataCoordinates::latitude(GeoDataCoordinates::Unit unit) const |
244 | { |
245 | switch (unit) |
246 | { |
247 | default: |
248 | case Radian: |
249 | return d->m_lat; |
250 | |
251 | case Degree: |
252 | return d->m_lat * RAD2DEG; |
253 | } |
254 | } |
255 | |
256 | qreal GeoDataCoordinates::latitude() const |
257 | { |
258 | return d->m_lat; |
259 | } |
260 | |
261 | //static |
262 | GeoDataCoordinates::Notation GeoDataCoordinates::defaultNotation() |
263 | { |
264 | return s_notation; |
265 | } |
266 | |
267 | //static |
268 | void GeoDataCoordinates::setDefaultNotation(GeoDataCoordinates::Notation notation) |
269 | { |
270 | s_notation = notation; |
271 | } |
272 | |
273 | //static |
274 | qreal GeoDataCoordinates::normalizeLon(qreal lon, GeoDataCoordinates::Unit unit) |
275 | { |
276 | qreal halfCircle; |
277 | |
278 | if (unit == GeoDataCoordinates::Radian) |
279 | { |
280 | halfCircle = M_PI3.14159265358979323846; |
281 | } |
282 | |
283 | else |
284 | { |
285 | halfCircle = 180; |
286 | } |
287 | |
288 | if (lon > halfCircle) |
289 | { |
290 | int cycles = (int)((lon + halfCircle) / (2 * halfCircle)); |
291 | return lon - (cycles * 2 * halfCircle); |
292 | } |
293 | |
294 | if (lon < -halfCircle) |
295 | { |
296 | int cycles = (int)((lon - halfCircle) / (2 * halfCircle)); |
297 | return lon - (cycles * 2 * halfCircle); |
298 | } |
299 | |
300 | return lon; |
301 | } |
302 | |
303 | //static |
304 | qreal GeoDataCoordinates::normalizeLat(qreal lat, GeoDataCoordinates::Unit unit) |
305 | { |
306 | qreal halfCircle; |
307 | |
308 | if (unit == GeoDataCoordinates::Radian) |
309 | { |
310 | halfCircle = M_PI3.14159265358979323846; |
311 | } |
312 | |
313 | else |
314 | { |
315 | halfCircle = 180; |
316 | } |
317 | |
318 | if (lat > (halfCircle / 2.0)) |
319 | { |
320 | int cycles = (int)((lat + halfCircle) / (2 * halfCircle)); |
321 | qreal temp; |
322 | |
323 | if (cycles == 0) // pi/2 < lat < pi |
324 | { |
325 | temp = halfCircle - lat; |
326 | } |
327 | |
328 | else |
329 | { |
330 | temp = lat - (cycles * 2 * halfCircle); |
331 | } |
332 | |
333 | if (temp > (halfCircle / 2.0)) |
334 | { |
335 | return (halfCircle - temp); |
336 | } |
337 | |
338 | if (temp < (-halfCircle / 2.0)) |
339 | { |
340 | return (-halfCircle - temp); |
341 | } |
342 | |
343 | return temp; |
344 | } |
345 | |
346 | if (lat < (-halfCircle / 2.0)) |
347 | { |
348 | int cycles = (int)((lat - halfCircle) / (2 * halfCircle)); |
349 | qreal temp; |
350 | |
351 | if (cycles == 0) |
352 | { |
353 | temp = -halfCircle - lat; |
354 | } |
355 | |
356 | else |
357 | { |
358 | temp = lat - (cycles * 2 * halfCircle); |
359 | } |
360 | |
361 | if (temp > (+halfCircle / 2.0)) |
362 | { |
363 | return (+halfCircle - temp); |
364 | } |
365 | |
366 | if (temp < (-halfCircle / 2.0)) |
367 | { |
368 | return (-halfCircle - temp); |
369 | } |
370 | |
371 | return temp; |
372 | } |
373 | |
374 | return lat; |
375 | } |
376 | |
377 | //static |
378 | void GeoDataCoordinates::normalizeLonLat(qreal& lon, qreal& lat, GeoDataCoordinates::Unit unit) |
379 | { |
380 | qreal halfCircle; |
381 | |
382 | if (unit == GeoDataCoordinates::Radian) |
383 | { |
384 | halfCircle = M_PI3.14159265358979323846; |
385 | } |
386 | |
387 | else |
388 | { |
389 | halfCircle = 180; |
390 | } |
391 | |
392 | if (lon > +halfCircle) |
393 | { |
394 | int cycles = (int)((lon + halfCircle) / (2 * halfCircle)); |
395 | lon = lon - (cycles * 2 * halfCircle); |
396 | } |
397 | |
398 | if (lon < -halfCircle) |
399 | { |
400 | int cycles = (int)((lon - halfCircle) / (2 * halfCircle)); |
401 | lon = lon - (cycles * 2 * halfCircle); |
402 | } |
403 | |
404 | if (lat > (+halfCircle / 2.0)) |
405 | { |
406 | int cycles = (int)((lat + halfCircle) / (2 * halfCircle)); |
407 | qreal temp; |
408 | |
409 | if (cycles == 0) // pi/2 < lat < pi |
410 | { |
411 | temp = halfCircle - lat; |
412 | } |
413 | |
414 | else |
415 | { |
416 | temp = lat - (cycles * 2 * halfCircle); |
417 | } |
418 | |
419 | if (temp > (+halfCircle / 2.0)) |
420 | { |
421 | lat = +halfCircle - temp; |
422 | } |
423 | |
424 | if (temp < (-halfCircle / 2.0)) |
425 | { |
426 | lat = -halfCircle - temp; |
427 | } |
428 | |
429 | lat = temp; |
430 | |
431 | if (lon > 0) |
432 | { |
433 | lon = -halfCircle + lon; |
434 | } |
435 | |
436 | else |
437 | { |
438 | lon = halfCircle + lon; |
439 | } |
440 | } |
441 | |
442 | if (lat < (-halfCircle / 2.0)) |
443 | { |
444 | int cycles = (int)((lat - halfCircle) / (2 * halfCircle)); |
445 | qreal temp; |
446 | |
447 | if (cycles == 0) |
448 | { |
449 | temp = -halfCircle - lat; |
450 | } |
451 | |
452 | else |
453 | { |
454 | temp = lat - (cycles * 2 * halfCircle); |
455 | } |
456 | |
457 | if (temp > (+halfCircle / 2.0)) |
458 | { |
459 | lat = +halfCircle - temp; |
460 | } |
461 | |
462 | if (temp < (-halfCircle / 2.0)) |
463 | { |
464 | lat = -halfCircle - temp; |
465 | } |
466 | |
467 | lat = temp; |
468 | |
469 | if (lon > 0) |
470 | { |
471 | lon = -halfCircle + lon; |
472 | } |
473 | |
474 | else |
475 | { |
476 | lon = halfCircle + lon; |
477 | } |
478 | } |
479 | |
480 | return; |
481 | } |
482 | |
483 | GeoDataCoordinates GeoDataCoordinates::fromString(const QString& string, bool& successful) |
484 | { |
485 | LonLatParser parser; |
486 | successful = parser.parse(string); |
487 | |
488 | if (successful) |
489 | { |
490 | return GeoDataCoordinates(parser.lon(), parser.lat(), 0, GeoDataCoordinates::Degree); |
491 | } |
492 | |
493 | else |
494 | { |
495 | return GeoDataCoordinates(); |
496 | } |
497 | } |
498 | |
499 | |
500 | QString GeoDataCoordinates::toString() const |
501 | { |
502 | return GeoDataCoordinates::toString(s_notation); |
503 | } |
504 | |
505 | QString GeoDataCoordinates::toString(GeoDataCoordinates::Notation notation, int precision) const |
506 | { |
507 | QString coordString; |
508 | |
509 | if (notation == GeoDataCoordinates::UTM) |
510 | { |
511 | int zoneNumber = GeoDataCoordinatesPrivate::lonLatToZone(d->m_lon, d->m_lat); |
512 | |
513 | // Handle lack of UTM zone number in the poles |
514 | const QString zoneString = (zoneNumber > 0) ? QString::number(zoneNumber) : QString(); |
515 | |
516 | QString bandString = GeoDataCoordinatesPrivate::lonLatToLatitudeBand(d->m_lon, d->m_lat); |
517 | |
518 | QString eastingString = QString::number(GeoDataCoordinatesPrivate::lonLatToEasting(d->m_lon, d->m_lat), 'f', 2); |
519 | QString northingString = QString::number(GeoDataCoordinatesPrivate::lonLatToNorthing(d->m_lon, d->m_lat), 'f', 2); |
520 | |
521 | return QString::fromUtf8("%1%2 %3 m E, %4 m N").arg(zoneString, bandString, eastingString, northingString); |
522 | } |
523 | |
524 | else |
525 | { |
526 | coordString = lonToString(d->m_lon, notation, Radian, precision) |
527 | + QLatin1String(", ") |
528 | + latToString(d->m_lat, notation, Radian, precision); |
529 | } |
530 | |
531 | return coordString; |
532 | } |
533 | |
534 | QString GeoDataCoordinates::lonToString(qreal lon, GeoDataCoordinates::Notation notation, |
535 | GeoDataCoordinates::Unit unit, |
536 | int precision, |
537 | char format) |
538 | { |
539 | if (notation == GeoDataCoordinates::UTM) |
540 | { |
541 | /** |
542 | * @FIXME: UTM needs lon + lat to know zone number and easting |
543 | * By now, this code returns the zone+easting of the point |
544 | * (lon, equator), but this can differ a lot at different locations |
545 | * See bug 347536 https://bugs.kde.org/show_bug.cgi?id=347536 |
546 | */ |
547 | |
548 | qreal lonRad = (unit == Radian) ? lon : lon * DEG2RAD; |
549 | |
550 | int zoneNumber = GeoDataCoordinatesPrivate::lonLatToZone(lonRad, 0); |
551 | |
552 | // Handle lack of UTM zone number in the poles |
553 | QString result = (zoneNumber > 0) ? QString::number(zoneNumber) : QString(); |
554 | |
555 | if (precision > 0) |
556 | { |
557 | QString eastingString = QString::number(GeoDataCoordinatesPrivate::lonLatToEasting(lonRad, 0), 'f', 2); |
558 | result += QString::fromUtf8(" %1 m E").arg(eastingString); |
559 | } |
560 | |
561 | return result; |
562 | } |
563 | |
564 | QString weString = (lon < 0) ? i18n("W")i18nd("digikam", "W") : i18n("E")i18nd("digikam", "E"); |
565 | |
566 | QString lonString; |
567 | |
568 | qreal lonDegF = (unit == Degree) ? fabs(lon) : fabs((qreal)(lon) * RAD2DEG); |
569 | |
570 | // Take care of -1 case |
571 | precision = (precision < 0) ? 5 : precision; |
572 | |
573 | if (notation == DMS || notation == DM) |
574 | { |
575 | int lonDeg = (int) lonDegF; |
576 | qreal lonMinF = 60 * (lonDegF - lonDeg); |
577 | int lonMin = (int) lonMinF; |
578 | qreal lonSecF = 60 * (lonMinF - lonMin); |
579 | int lonSec = (int) lonSecF; |
580 | |
581 | // Adjustment for fuzziness (like 49.999999999999999999999) |
582 | if (precision == 0) |
583 | { |
584 | lonDeg = qRound(lonDegF); |
585 | } |
586 | |
587 | else if (precision <= 2) |
588 | { |
589 | lonMin = qRound(lonMinF); |
590 | } |
591 | |
592 | else if (precision <= 4 && notation == DMS) |
593 | { |
594 | lonSec = qRound(lonSecF); |
595 | } |
596 | |
597 | else |
598 | { |
599 | if (notation == DMS) |
600 | { |
601 | lonSec = lonSecF = qRound(lonSecF * qPow(10, precision - 4)) / qPow(10, precision - 4); |
602 | } |
603 | |
604 | else |
605 | { |
606 | lonMin = lonMinF = qRound(lonMinF * qPow(10, precision - 2)) / qPow(10, precision - 2); |
607 | } |
608 | } |
609 | |
610 | if (lonSec > 59 && notation == DMS) |
611 | { |
612 | lonSec = lonSecF = 0; |
613 | lonMin = lonMinF = lonMinF + 1; |
614 | } |
615 | |
616 | if (lonMin > 59) |
617 | { |
618 | lonMin = lonMinF = 0; |
619 | lonDeg = lonDegF = lonDegF + 1; |
620 | } |
621 | |
622 | // Evaluate the string |
623 | lonString = QString::fromUtf8("%1\xc2\xb0").arg(lonDeg, 3, 10, QLatin1Char(' ')); |
624 | |
625 | if (precision == 0) |
626 | { |
627 | return lonString + weString; |
628 | } |
629 | |
630 | if (notation == DMS || precision < 3) |
631 | { |
632 | lonString += QString::fromUtf8(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0')); |
633 | } |
634 | |
635 | if (precision < 3) |
636 | { |
637 | return lonString + weString; |
638 | } |
639 | |
640 | if (notation == DMS) |
641 | { |
642 | // Includes -1 case! |
643 | if (precision < 5) |
644 | { |
645 | lonString += QString::fromUtf8(" %3\"").arg(lonSec, 2, 'f', 0, QLatin1Char('0')); |
646 | return lonString + weString; |
647 | } |
648 | |
649 | lonString += QString::fromUtf8(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0')); |
650 | } |
651 | |
652 | else |
653 | { |
654 | lonString += QString::fromUtf8(" %L3'").arg(lonMinF, precision + 1, 'f', precision - 2, QLatin1Char('0')); |
655 | } |
656 | } |
657 | |
658 | else if (notation == GeoDataCoordinates::Decimal) |
659 | { |
660 | lonString = QString::fromUtf8("%L1\xc2\xb0").arg(lonDegF, 4 + precision, format, precision, QLatin1Char(' ')); |
661 | } |
662 | |
663 | else if (notation == GeoDataCoordinates::Astro) |
664 | { |
665 | if (lon < 0) |
666 | { |
667 | lon += (unit == Degree) ? 360 : 2 * M_PI3.14159265358979323846; |
668 | } |
669 | |
670 | qreal lonHourF = (unit == Degree) ? fabs(lon / 15.0) : fabs((qreal)(lon / 15.0) * RAD2DEG); |
671 | int lonHour = (int) lonHourF; |
672 | qreal lonMinF = 60 * (lonHourF - lonHour); |
673 | int lonMin = (int) lonMinF; |
674 | qreal lonSecF = 60 * (lonMinF - lonMin); |
675 | int lonSec = (int) lonSecF; |
676 | |
677 | // Adjustment for fuzziness (like 49.999999999999999999999) |
678 | if (precision == 0) |
679 | { |
680 | lonHour = qRound(lonHourF); |
681 | } |
682 | |
683 | else if (precision <= 2) |
684 | { |
685 | lonMin = qRound(lonMinF); |
686 | } |
687 | |
688 | else if (precision <= 4) |
689 | { |
690 | lonSec = qRound(lonSecF); |
691 | } |
692 | |
693 | else |
694 | { |
695 | lonSec = lonSecF = qRound(lonSecF * qPow(10, precision - 4)) / qPow(10, precision - 4); |
696 | } |
697 | |
698 | if (lonSec > 59) |
699 | { |
700 | lonSec = lonSecF = 0; |
701 | lonMin = lonMinF = lonMinF + 1; |
702 | } |
703 | |
704 | if (lonMin > 59) |
705 | { |
706 | lonMin = lonMinF = 0; |
Although the value stored to 'lonMinF' is used in the enclosing expression, the value is never actually read from 'lonMinF' | |
707 | lonHour = lonHourF = lonHourF + 1; |
708 | } |
709 | |
710 | // Evaluate the string |
711 | lonString = QString::fromUtf8("%1h").arg(lonHour, 3, 10, QLatin1Char(' ')); |
712 | |
713 | if (precision == 0) |
714 | { |
715 | return lonString; |
716 | } |
717 | |
718 | lonString += QString::fromUtf8(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0')); |
719 | |
720 | if (precision < 3) |
721 | { |
722 | return lonString; |
723 | } |
724 | |
725 | // Includes -1 case! |
726 | if (precision < 5) |
727 | { |
728 | lonString += QString::fromUtf8(" %3\"").arg(lonSec, 2, 'f', 0, QLatin1Char('0')); |
729 | return lonString; |
730 | } |
731 | |
732 | lonString += QString::fromUtf8(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0')); |
733 | return lonString; |
734 | } |
735 | |
736 | return lonString + weString; |
737 | } |
738 | |
739 | QString GeoDataCoordinates::lonToString() const |
740 | { |
741 | return GeoDataCoordinates::lonToString(d->m_lon, s_notation); |
742 | } |
743 | |
744 | QString GeoDataCoordinates::latToString(qreal lat, GeoDataCoordinates::Notation notation, |
745 | GeoDataCoordinates::Unit unit, |
746 | int precision, |
747 | char format) |
748 | { |
749 | if (notation == GeoDataCoordinates::UTM) |
750 | { |
751 | /** |
752 | * @FIXME: UTM needs lon + lat to know latitude band and northing |
753 | * By now, this code returns the band+northing of the point |
754 | * (meridian, lat), but this can differ a lot at different locations |
755 | * See bug 347536 https://bugs.kde.org/show_bug.cgi?id=347536 |
756 | */ |
757 | |
758 | qreal latRad = (unit == Radian) ? lat : lat * DEG2RAD; |
759 | |
760 | QString result = GeoDataCoordinatesPrivate::lonLatToLatitudeBand(0, latRad); |
761 | |
762 | if (precision > 0) |
763 | { |
764 | QString northingString = QString::number(GeoDataCoordinatesPrivate::lonLatToNorthing(0, latRad), 'f', 2); |
765 | result += QString::fromUtf8(" %1 m N").arg(northingString); |
766 | } |
767 | |
768 | return result; |
769 | } |
770 | |
771 | QString pmString; |
772 | QString nsString; |
773 | |
774 | if (notation == GeoDataCoordinates::Astro) |
775 | { |
776 | pmString = (lat > 0) ? QLatin1String("+") : QLatin1String("-"); |
777 | } |
778 | |
779 | else |
780 | { |
781 | nsString = (lat > 0) ? i18n("N")i18nd("digikam", "N") : i18n("S")i18nd("digikam", "S"); |
782 | } |
783 | |
784 | QString latString; |
785 | |
786 | qreal latDegF = (unit == Degree) ? fabs(lat) : fabs((qreal)(lat) * RAD2DEG); |
787 | |
788 | // Take care of -1 case |
789 | precision = (precision < 0) ? 5 : precision; |
790 | |
791 | if (notation == DMS || notation == DM || notation == Astro) |
792 | { |
793 | int latDeg = (int) latDegF; |
794 | qreal latMinF = 60 * (latDegF - latDeg); |
795 | int latMin = (int) latMinF; |
796 | qreal latSecF = 60 * (latMinF - latMin); |
797 | int latSec = (int) latSecF; |
798 | |
799 | // Adjustment for fuzziness (like 49.999999999999999999999) |
800 | if (precision == 0) |
801 | { |
802 | latDeg = qRound(latDegF); |
803 | } |
804 | |
805 | else if (precision <= 2) |
806 | { |
807 | latMin = qRound(latMinF); |
808 | } |
809 | |
810 | else if (precision <= 4 && notation == DMS) |
811 | { |
812 | latSec = qRound(latSecF); |
813 | } |
814 | |
815 | else |
816 | { |
817 | if (notation == DMS || notation == Astro) |
818 | { |
819 | latSec = latSecF = qRound(latSecF * qPow(10, precision - 4)) / qPow(10, precision - 4); |
820 | } |
821 | |
822 | else |
823 | { |
824 | latMin = latMinF = qRound(latMinF * qPow(10, precision - 2)) / qPow(10, precision - 2); |
825 | } |
826 | } |
827 | |
828 | if (latSec > 59 && (notation == DMS || notation == Astro)) |
829 | { |
830 | latSecF = 0; |
831 | latSec = latSecF; |
832 | latMin = latMin + 1; |
833 | } |
834 | |
835 | if (latMin > 59) |
836 | { |
837 | latMinF = 0; |
838 | latMin = latMinF; |
839 | latDeg = latDeg + 1; |
840 | } |
841 | |
842 | // Evaluate the string |
843 | latString = QString::fromUtf8("%1\xc2\xb0").arg(latDeg, 3, 10, QLatin1Char(' ')); |
844 | |
845 | if (precision == 0) |
846 | { |
847 | return pmString + latString + nsString; |
848 | } |
849 | |
850 | if (notation == DMS || notation == Astro || precision < 3) |
851 | { |
852 | latString += QString::fromUtf8(" %2\'").arg(latMin, 2, 10, QLatin1Char('0')); |
853 | } |
854 | |
855 | if (precision < 3) |
856 | { |
857 | return pmString + latString + nsString; |
858 | } |
859 | |
860 | if (notation == DMS || notation == Astro) |
861 | { |
862 | // Includes -1 case! |
863 | if (precision < 5) |
864 | { |
865 | latString += QString::fromUtf8(" %3\"").arg(latSec, 2, 'f', 0, QLatin1Char('0')); |
866 | return latString + nsString; |
867 | } |
868 | |
869 | latString += QString::fromUtf8(" %L3\"").arg(latSecF, precision - 1, 'f', precision - 4, QLatin1Char('0')); |
870 | } |
871 | |
872 | else |
873 | { |
874 | latString += QString::fromUtf8(" %L3'").arg(latMinF, precision + 1, 'f', precision - 2, QLatin1Char('0')); |
875 | } |
876 | } |
877 | |
878 | else // notation = GeoDataCoordinates::Decimal |
879 | { |
880 | latString = QString::fromUtf8("%L1\xc2\xb0").arg(latDegF, 4 + precision, format, precision, QLatin1Char(' ')); |
881 | } |
882 | |
883 | return pmString + latString + nsString; |
884 | } |
885 | |
886 | QString GeoDataCoordinates::latToString() const |
887 | { |
888 | return GeoDataCoordinates::latToString(d->m_lat, s_notation); |
889 | } |
890 | |
891 | bool GeoDataCoordinates::operator==(const GeoDataCoordinates& rhs) const |
892 | { |
893 | return *d == *rhs.d; |
894 | } |
895 | |
896 | bool GeoDataCoordinates::operator!=(const GeoDataCoordinates& rhs) const |
897 | { |
898 | return *d != *rhs.d; |
899 | } |
900 | |
901 | void GeoDataCoordinates::setAltitude(const qreal altitude) |
902 | { |
903 | detach(); |
904 | d->m_altitude = altitude; |
905 | } |
906 | |
907 | qreal GeoDataCoordinates::altitude() const |
908 | { |
909 | return d->m_altitude; |
910 | } |
911 | |
912 | int GeoDataCoordinates::utmZone() const |
913 | { |
914 | return GeoDataCoordinatesPrivate::lonLatToZone(d->m_lon, d->m_lat); |
915 | } |
916 | |
917 | qreal GeoDataCoordinates::utmEasting() const |
918 | { |
919 | return GeoDataCoordinatesPrivate::lonLatToEasting(d->m_lon, d->m_lat); |
920 | } |
921 | |
922 | QString GeoDataCoordinates::utmLatitudeBand() const |
923 | { |
924 | return GeoDataCoordinatesPrivate::lonLatToLatitudeBand(d->m_lon, d->m_lat); |
925 | } |
926 | |
927 | qreal GeoDataCoordinates::utmNorthing() const |
928 | { |
929 | return GeoDataCoordinatesPrivate::lonLatToNorthing(d->m_lon, d->m_lat); |
930 | } |
931 | |
932 | quint8 GeoDataCoordinates::detail() const |
933 | { |
934 | return d->m_detail; |
935 | } |
936 | |
937 | void GeoDataCoordinates::setDetail(quint8 detail) |
938 | { |
939 | detach(); |
940 | d->m_detail = detail; |
941 | } |
942 | |
943 | GeoDataCoordinates GeoDataCoordinates::rotateAround(const GeoDataCoordinates& axis, qreal angle, Unit unit) const |
944 | { |
945 | const Quaternion quatAxis = Quaternion::fromEuler(-axis.latitude(), axis.longitude(), 0); |
946 | const Quaternion rotationAmount = Quaternion::fromEuler(0, 0, unit == Radian ? angle : angle * DEG2RAD); |
947 | const Quaternion resultAxis = quatAxis * rotationAmount * quatAxis.inverse(); |
948 | |
949 | return rotateAround(resultAxis); |
950 | } |
951 | |
952 | GeoDataCoordinates GeoDataCoordinates::rotateAround(const Quaternion& rotAxis) const |
953 | { |
954 | Quaternion rotatedQuat = quaternion(); |
955 | rotatedQuat.rotateAroundAxis(rotAxis); |
956 | qreal rotatedLon, rotatedLat; |
957 | rotatedQuat.getSpherical(rotatedLon, rotatedLat); |
958 | return GeoDataCoordinates(rotatedLon, rotatedLat, altitude()); |
959 | } |
960 | |
961 | qreal GeoDataCoordinates::bearing(const GeoDataCoordinates& other, Unit unit, BearingType type) const |
962 | { |
963 | if (type == FinalBearing) |
964 | { |
965 | double const offset = unit == Degree ? 180.0 : M_PI3.14159265358979323846; |
966 | return offset + other.bearing(*this, unit, InitialBearing); |
967 | } |
968 | |
969 | qreal const delta = other.d->m_lon - d->m_lon; |
970 | double const bearing = atan2(sin(delta) * cos(other.d->m_lat), |
971 | cos(d->m_lat) * sin(other.d->m_lat) - sin(d->m_lat) * cos(other.d->m_lat) * cos(delta)); |
972 | return unit == Radian ? bearing : bearing * RAD2DEG; |
973 | } |
974 | |
975 | GeoDataCoordinates GeoDataCoordinates::moveByBearing(qreal bearing, qreal distance) const |
976 | { |
977 | qreal newLat = asin(sin(d->m_lat) * cos(distance) + |
978 | cos(d->m_lat) * sin(distance) * cos(bearing)); |
979 | qreal newLon = d->m_lon + atan2(sin(bearing) * sin(distance) * cos(d->m_lat), |
980 | cos(distance) - sin(d->m_lat) * sin(newLat)); |
981 | |
982 | return GeoDataCoordinates(newLon, newLat); |
983 | } |
984 | |
985 | const Quaternion& GeoDataCoordinates::quaternion() const |
986 | { |
987 | if (d->m_q == nullptr) |
988 | { |
989 | d->m_q = new Quaternion(Quaternion::fromSpherical(d->m_lon, d->m_lat)); |
990 | } |
991 | |
992 | return *d->m_q; |
993 | } |
994 | |
995 | GeoDataCoordinates GeoDataCoordinates::interpolate(const GeoDataCoordinates& target, double t_) const |
996 | { |
997 | double const t = qBound(0.0, t_, 1.0); |
998 | Quaternion const quat = Quaternion::slerp(quaternion(), target.quaternion(), t); |
999 | qreal lon, lat; |
1000 | quat.getSpherical(lon, lat); |
1001 | double const alt = (1.0 - t) * d->m_altitude + t * target.d->m_altitude; |
1002 | return GeoDataCoordinates(lon, lat, alt); |
1003 | } |
1004 | |
1005 | GeoDataCoordinates GeoDataCoordinates::nlerp(const GeoDataCoordinates& target, double t) const |
1006 | { |
1007 | qreal lon = 0.0; |
1008 | qreal lat = 0.0; |
1009 | |
1010 | const Quaternion itpos = Quaternion::nlerp(quaternion(), target.quaternion(), t); |
1011 | itpos.getSpherical(lon, lat); |
1012 | |
1013 | const qreal altitude = 0.5 * (d->m_altitude + target.altitude()); |
1014 | |
1015 | return GeoDataCoordinates(lon, lat, altitude); |
1016 | } |
1017 | |
1018 | GeoDataCoordinates GeoDataCoordinates::interpolate(const GeoDataCoordinates& before, const GeoDataCoordinates& target, const GeoDataCoordinates& after, double t_) const |
1019 | { |
1020 | double const t = qBound(0.0, t_, 1.0); |
1021 | Quaternion const b1 = GeoDataCoordinatesPrivate::basePoint(before.quaternion(), quaternion(), target.quaternion()); |
1022 | Quaternion const a2 = GeoDataCoordinatesPrivate::basePoint(quaternion(), target.quaternion(), after.quaternion()); |
1023 | Quaternion const a = Quaternion::slerp(quaternion(), target.quaternion(), t); |
1024 | Quaternion const b = Quaternion::slerp(b1, a2, t); |
1025 | Quaternion c = Quaternion::slerp(a, b, 2 * t * (1.0 - t)); |
1026 | qreal lon, lat; |
1027 | c.getSpherical(lon, lat); |
1028 | // @todo spline interpolation of altitude? |
1029 | double const alt = (1.0 - t) * d->m_altitude + t * target.d->m_altitude; |
1030 | return GeoDataCoordinates(lon, lat, alt); |
1031 | } |
1032 | |
1033 | bool GeoDataCoordinates::isPole(Pole pole) const |
1034 | { |
1035 | // Evaluate the most likely case first: |
1036 | // The case where we haven't hit the pole and where our latitude is normalized |
1037 | // to the range of 90 deg S ... 90 deg N |
1038 | if (fabs((qreal) 2.0 * d->m_lat) < M_PI3.14159265358979323846) |
1039 | { |
1040 | return false; |
1041 | } |
1042 | |
1043 | else |
1044 | { |
1045 | if (fabs((qreal) 2.0 * d->m_lat) == M_PI3.14159265358979323846) |
1046 | { |
1047 | // Ok, we have hit a pole. Now let's check whether it's the one we've asked for: |
1048 | if (pole == AnyPole) |
1049 | { |
1050 | return true; |
1051 | } |
1052 | |
1053 | else |
1054 | { |
1055 | if (pole == NorthPole && 2.0 * d->m_lat == +M_PI3.14159265358979323846) |
1056 | { |
1057 | return true; |
1058 | } |
1059 | |
1060 | if (pole == SouthPole && 2.0 * d->m_lat == -M_PI3.14159265358979323846) |
1061 | { |
1062 | return true; |
1063 | } |
1064 | |
1065 | return false; |
1066 | } |
1067 | } |
1068 | |
1069 | // |
1070 | else |
1071 | { |
1072 | // FIXME: Should we just normalize latitude and longitude and be done? |
1073 | // While this might work well for persistent data it would create some |
1074 | // possible overhead for temporary data, so this needs careful thinking. |
1075 | qCDebug(DIGIKAM_MARBLE_LOG)for (QLoggingCategoryMacroHolder<QtDebugMsg> qt_category ((DIGIKAM_MARBLE_LOG)()); qt_category; qt_category.control = false ) QMessageLogger(static_cast<const char *>("/home/gilles/devel/GIT/8.x/core/utilities/geolocation/engine/geodata/data/GeoDataCoordinates.cpp" ), 1075, static_cast<const char *>(__PRETTY_FUNCTION__) , qt_category.name()).debug() << "GeoDataCoordinates not normalized!"; |
1076 | |
1077 | // Only as a last resort we cover the unlikely case where |
1078 | // the latitude is not normalized to the range of |
1079 | // 90 deg S ... 90 deg N |
1080 | if (fabs((qreal) 2.0 * normalizeLat(d->m_lat)) < M_PI3.14159265358979323846) |
1081 | { |
1082 | return false; |
1083 | } |
1084 | |
1085 | else |
1086 | { |
1087 | // Ok, we have hit a pole. Now let's check whether it's the one we've asked for: |
1088 | if (pole == AnyPole) |
1089 | { |
1090 | return true; |
1091 | } |
1092 | |
1093 | else |
1094 | { |
1095 | if (pole == NorthPole && 2.0 * d->m_lat == +M_PI3.14159265358979323846) |
1096 | { |
1097 | return true; |
1098 | } |
1099 | |
1100 | if (pole == SouthPole && 2.0 * d->m_lat == -M_PI3.14159265358979323846) |
1101 | { |
1102 | return true; |
1103 | } |
1104 | |
1105 | return false; |
1106 | } |
1107 | } |
1108 | } |
1109 | } |
1110 | } |
1111 | |
1112 | qreal GeoDataCoordinates::sphericalDistanceTo(const GeoDataCoordinates& other) const |
1113 | { |
1114 | qreal lon2, lat2; |
1115 | other.geoCoordinates(lon2, lat2); |
1116 | |
1117 | // FIXME: Take the altitude into account! |
1118 | |
1119 | return distanceSphere(d->m_lon, d->m_lat, lon2, lat2); |
1120 | } |
1121 | |
1122 | GeoDataCoordinates& GeoDataCoordinates::operator=(const GeoDataCoordinates& other) |
1123 | { |
1124 | qAtomicAssign(d, other.d); |
1125 | return *this; |
1126 | } |
1127 | |
1128 | void GeoDataCoordinates::pack(QDataStream& stream) const |
1129 | { |
1130 | stream << d->m_lon; |
1131 | stream << d->m_lat; |
1132 | stream << d->m_altitude; |
1133 | } |
1134 | |
1135 | void GeoDataCoordinates::unpack(QDataStream& stream) |
1136 | { |
1137 | // call detach even though it shouldn't be needed - one never knows |
1138 | detach(); |
1139 | stream >> d->m_lon; |
1140 | stream >> d->m_lat; |
1141 | stream >> d->m_altitude; |
1142 | } |
1143 | |
1144 | Quaternion GeoDataCoordinatesPrivate::basePoint(const Quaternion& q1, const Quaternion& q2, const Quaternion& q3) |
1145 | { |
1146 | Quaternion const a = (q2.inverse() * q3).log(); |
1147 | Quaternion const b = (q2.inverse() * q1).log(); |
1148 | return q2 * ((a + b) * -0.25).exp(); |
1149 | } |
1150 | |
1151 | |
1152 | |
1153 | qreal GeoDataCoordinatesPrivate::arcLengthOfMeridian(qreal phi) |
1154 | { |
1155 | // Precalculate n |
1156 | qreal const n = (GeoDataCoordinatesPrivate::sm_semiMajorAxis - GeoDataCoordinatesPrivate::sm_semiMinorAxis) |
1157 | / (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis); |
1158 | |
1159 | // Precalculate alpha |
1160 | qreal const alpha = ((GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis) / 2.0) |
1161 | * (1.0 + (qPow(n, 2.0) / 4.0) + (qPow(n, 4.0) / 64.0)); |
1162 | |
1163 | // Precalculate beta |
1164 | qreal const beta = (-3.0 * n / 2.0) |
1165 | + (9.0 * qPow(n, 3.0) / 16.0) |
1166 | + (-3.0 * qPow(n, 5.0) / 32.0); |
1167 | |
1168 | // Precalculate gamma |
1169 | qreal const gamma = (15.0 * qPow(n, 2.0) / 16.0) |
1170 | + (-15.0 * qPow(n, 4.0) / 32.0); |
1171 | |
1172 | // Precalculate delta |
1173 | qreal const delta = (-35.0 * qPow(n, 3.0) / 48.0) |
1174 | + (105.0 * qPow(n, 5.0) / 256.0); |
1175 | |
1176 | // Precalculate epsilon |
1177 | qreal const epsilon = (315.0 * qPow(n, 4.0) / 512.0); |
1178 | |
1179 | // Now calculate the sum of the series and return |
1180 | qreal const result = alpha * (phi + (beta * qSin(2.0 * phi)) |
1181 | + (gamma * qSin(4.0 * phi)) |
1182 | + (delta * qSin(6.0 * phi)) |
1183 | + (epsilon * qSin(8.0 * phi))); |
1184 | |
1185 | return result; |
1186 | } |
1187 | |
1188 | qreal GeoDataCoordinatesPrivate::centralMeridianUTM(qreal zone) |
1189 | { |
1190 | return DEG2RAD * (-183.0 + (zone * 6.0)); |
1191 | } |
1192 | |
1193 | qreal GeoDataCoordinatesPrivate::footpointLatitude(qreal northing) |
1194 | { |
1195 | // Precalculate n (Eq. 10.18) |
1196 | qreal const n = (GeoDataCoordinatesPrivate::sm_semiMajorAxis - GeoDataCoordinatesPrivate::sm_semiMinorAxis) |
1197 | / (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis); |
1198 | |
1199 | // Precalculate alpha (Eq. 10.22) |
1200 | // (Same as alpha in Eq. 10.17) |
1201 | qreal const alpha = ((GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis) / 2.0) |
1202 | * (1 + (qPow(n, 2.0) / 4) + (qPow(n, 4.0) / 64)); |
1203 | |
1204 | // Precalculate y (Eq. 10.23) |
1205 | qreal const y = northing / alpha; |
1206 | |
1207 | // Precalculate beta (Eq. 10.22) |
1208 | qreal const beta = (3.0 * n / 2.0) + (-27.0 * qPow(n, 3.0) / 32.0) |
1209 | + (269.0 * qPow(n, 5.0) / 512.0); |
1210 | |
1211 | // Precalculate gamma (Eq. 10.22) |
1212 | qreal const gamma = (21.0 * qPow(n, 2.0) / 16.0) |
1213 | + (-55.0 * qPow(n, 4.0) / 32.0); |
1214 | |
1215 | // Precalculate delta (Eq. 10.22) |
1216 | qreal const delta = (151.0 * qPow(n, 3.0) / 96.0) |
1217 | + (-417.0 * qPow(n, 5.0) / 128.0); |
1218 | |
1219 | // Precalculate epsilon (Eq. 10.22) |
1220 | qreal const epsilon = (1097.0 * qPow(n, 4.0) / 512.0); |
1221 | |
1222 | // Now calculate the sum of the series (Eq. 10.21) |
1223 | qreal const result = y + (beta * qSin(2.0 * y)) |
1224 | + (gamma * qSin(4.0 * y)) |
1225 | + (delta * qSin(6.0 * y)) |
1226 | + (epsilon * qSin(8.0 * y)); |
1227 | |
1228 | return result; |
1229 | } |
1230 | |
1231 | QPointF GeoDataCoordinatesPrivate::mapLonLatToXY(qreal lambda, qreal phi, qreal lambda0) |
1232 | { |
1233 | // Equation (10.15) |
1234 | |
1235 | // Precalculate second numerical eccentricity |
1236 | const qreal ep2 = (qPow(GeoDataCoordinatesPrivate::sm_semiMajorAxis, 2.0) - qPow(GeoDataCoordinatesPrivate::sm_semiMinorAxis, 2.0)) |
1237 | / qPow(GeoDataCoordinatesPrivate::sm_semiMinorAxis, 2.0); |
1238 | |
1239 | // Precalculate the square of nu, just an auxilar quantity |
1240 | const qreal nu2 = ep2 * qPow(qCos(phi), 2.0); |
1241 | |
1242 | // Precalculate the radius of curvature in prime vertical |
1243 | const qreal N = qPow(GeoDataCoordinatesPrivate::sm_semiMajorAxis, 2.0) / (GeoDataCoordinatesPrivate::sm_semiMinorAxis * qSqrt(1 + nu2)); |
1244 | |
1245 | // Precalculate the tangent of phi and its square |
1246 | const qreal t = qTan(phi); |
1247 | const qreal t2 = t * t; |
1248 | |
1249 | // Precalculate longitude difference |
1250 | const qreal l = lambda - lambda0; |
1251 | |
1252 | /* |
1253 | * Precalculate coefficients for l**n in the equations below |
1254 | * so a normal human being can read the expressions for easting |
1255 | * and northing |
1256 | * -- l**1 and l**2 have coefficients of 1.0 |
1257 | * |
1258 | * The actual used coefficients starts at coef[1], just to |
1259 | * follow the meaningful nomenclature in equation 10.15 |
1260 | * (coef[n] corresponds to qPow(l,n) factor) |
1261 | */ |
1262 | |
1263 | const qreal coef1 = 1; |
1264 | const qreal coef2 = 1; |
1265 | const qreal coef3 = 1.0 - t2 + nu2; |
1266 | const qreal coef4 = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2); |
1267 | const qreal coef5 = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 - 58.0 * t2 * nu2; |
1268 | const qreal coef6 = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 - 330.0 * t2 * nu2; |
1269 | const qreal coef7 = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2); |
1270 | const qreal coef8 = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2); |
1271 | |
1272 | // Calculate easting (x) |
1273 | const qreal easting = N * qCos(phi) * coef1 * l |
1274 | + (N / 6.0 * qPow(qCos(phi), 3.0) * coef3 * qPow(l, 3.0)) |
1275 | + (N / 120.0 * qPow(qCos(phi), 5.0) * coef5 * qPow(l, 5.0)) |
1276 | + (N / 5040.0 * qPow(qCos(phi), 7.0) * coef7 * qPow(l, 7.0)); |
1277 | |
1278 | // Calculate northing (y) |
1279 | const qreal northing = arcLengthOfMeridian(phi) |
1280 | + (t / 2.0 * N * qPow(qCos(phi), 2.0) * coef2 * qPow(l, 2.0)) |
1281 | + (t / 24.0 * N * qPow(qCos(phi), 4.0) * coef4 * qPow(l, 4.0)) |
1282 | + (t / 720.0 * N * qPow(qCos(phi), 6.0) * coef6 * qPow(l, 6.0)) |
1283 | + (t / 40320.0 * N * qPow(qCos(phi), 8.0) * coef8 * qPow(l, 8.0)); |
1284 | |
1285 | return QPointF(easting, northing); |
1286 | } |
1287 | |
1288 | int GeoDataCoordinatesPrivate::lonLatToZone(qreal lon, qreal lat) |
1289 | { |
1290 | // Converts lon and lat to degrees |
1291 | qreal lonDeg = lon * RAD2DEG; |
1292 | qreal latDeg = lat * RAD2DEG; |
1293 | |
1294 | /* Round the value of the longitude when the distance to the nearest integer |
1295 | * is less than 0.0000001. This avoids fuzzy values such as -114.0000000001, which |
1296 | * can produce a misbehaviour when calculating the zone associated at the borders |
1297 | * of the zone intervals (for example, the interval [-114, -108[ is associated with |
1298 | * zone number 12; if the following rounding is not done, the value returned by |
1299 | * lonLatToZone(114,0) is 11 instead of 12, as the function actually receives |
1300 | * -114.0000000001, which is in the interval [-120,-114[, associated to zone 11 |
1301 | */ |
1302 | qreal precision = 0.0000001; |
1303 | |
1304 | if (qAbs(lonDeg - qFloor(lonDeg)) < precision || qAbs(lonDeg - qCeil(lonDeg)) < precision) |
1305 | { |
1306 | lonDeg = qRound(lonDeg); |
1307 | } |
1308 | |
1309 | // There is no numbering associated to the poles, special value 0 is returned. |
1310 | if (latDeg < -80 || latDeg > 84) |
1311 | { |
1312 | return 0; |
1313 | } |
1314 | |
1315 | // Obtains the zone number handling all the so called "exceptions" |
1316 | // See problem: https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions |
1317 | // See solution: https://gis.stackexchange.com/questions/13291/computing-utm-zone-from-lat-long-point |
1318 | |
1319 | // General |
1320 | int zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1; |
1321 | |
1322 | // Southwest Norway |
1323 | if (latDeg >= 56 && latDeg < 64 && lonDeg >= 3 && lonDeg < 12) |
1324 | { |
1325 | zoneNumber = 32; |
1326 | } |
1327 | |
1328 | // Svalbard |
1329 | if (latDeg >= 72 && latDeg < 84) |
1330 | { |
1331 | if (lonDeg >= 0 && lonDeg < 9) |
1332 | { |
1333 | zoneNumber = 31; |
1334 | } |
1335 | |
1336 | else if (lonDeg >= 9 && lonDeg < 21) |
1337 | { |
1338 | zoneNumber = 33; |
1339 | } |
1340 | |
1341 | else if (lonDeg >= 21 && lonDeg < 33) |
1342 | { |
1343 | zoneNumber = 35; |
1344 | } |
1345 | |
1346 | else if (lonDeg >= 33 && lonDeg < 42) |
1347 | { |
1348 | zoneNumber = 37; |
1349 | } |
1350 | } |
1351 | |
1352 | return zoneNumber; |
1353 | } |
1354 | |
1355 | qreal GeoDataCoordinatesPrivate::lonLatToEasting(qreal lon, qreal lat) |
1356 | { |
1357 | int zoneNumber = lonLatToZone(lon, lat); |
1358 | |
1359 | if (zoneNumber == 0) |
1360 | { |
1361 | qreal lonDeg = lon * RAD2DEG; |
1362 | zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1; |
1363 | } |
1364 | |
1365 | QPointF coordinates = GeoDataCoordinatesPrivate::mapLonLatToXY(lon, lat, GeoDataCoordinatesPrivate::centralMeridianUTM(zoneNumber)); |
1366 | |
1367 | // Adjust easting and northing for UTM system |
1368 | qreal easting = coordinates.x() * GeoDataCoordinatesPrivate::sm_utmScaleFactor + 500000.0; |
1369 | |
1370 | return easting; |
1371 | } |
1372 | |
1373 | QString GeoDataCoordinatesPrivate::lonLatToLatitudeBand(qreal lon, qreal lat) |
1374 | { |
1375 | // Obtains the latitude bands handling all the so called "exceptions" |
1376 | |
1377 | // Converts lon and lat to degrees |
1378 | qreal lonDeg = lon * RAD2DEG; |
1379 | qreal latDeg = lat * RAD2DEG; |
1380 | |
1381 | // Regular latitude bands between 80 S and 80 N (that is, between 10 and 170 in the [0,180] interval) |
1382 | int bandLetterIndex = 24; //Avoids "may be used uninitialized" warning |
1383 | |
1384 | if (latDeg < -80) |
1385 | { |
1386 | // South pole (A for zones 1-30, B for zones 31-60) |
1387 | bandLetterIndex = ((lonDeg + 180) < 6 * 31) ? 0 : 1; |
1388 | } |
1389 | |
1390 | else if (latDeg >= -80 && latDeg <= 80) |
1391 | { |
1392 | // General (+2 because the general lettering starts in C) |
1393 | bandLetterIndex = static_cast<int>((latDeg + 80.0) / 8.0) + 2; |
1394 | } |
1395 | |
1396 | else if (latDeg >= 80 && latDeg < 84) |
1397 | { |
1398 | // Band X is extended 4 more degrees |
1399 | bandLetterIndex = 21; |
1400 | } |
1401 | |
1402 | else if (latDeg >= 84) |
1403 | { |
1404 | // North pole (Y for zones 1-30, Z for zones 31-60) |
1405 | bandLetterIndex = ((lonDeg + 180) < 6 * 31) ? 22 : 23; |
1406 | } |
1407 | |
1408 | return QStringLiteral("ABCDEFGHJKLMNPQRSTUVWXYZ?")(QString(QtPrivate::qMakeStringPrivate(u"" "ABCDEFGHJKLMNPQRSTUVWXYZ?" ))).at(bandLetterIndex); |
1409 | } |
1410 | |
1411 | qreal GeoDataCoordinatesPrivate::lonLatToNorthing(qreal lon, qreal lat) |
1412 | { |
1413 | int zoneNumber = lonLatToZone(lon, lat); |
1414 | |
1415 | if (zoneNumber == 0) |
1416 | { |
1417 | qreal lonDeg = lon * RAD2DEG; |
1418 | zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1; |
1419 | } |
1420 | |
1421 | QPointF coordinates = GeoDataCoordinatesPrivate::mapLonLatToXY(lon, lat, GeoDataCoordinatesPrivate::centralMeridianUTM(zoneNumber)); |
1422 | |
1423 | qreal northing = coordinates.y() * GeoDataCoordinatesPrivate::sm_utmScaleFactor; |
1424 | |
1425 | if (northing < 0.0) |
1426 | { |
1427 | northing += 10000000.0; |
1428 | } |
1429 | |
1430 | return northing; |
1431 | } |
1432 | |
1433 | size_t qHash(const GeoDataCoordinates& coordinates) |
1434 | { |
1435 | uint seed = ::qHash(coordinates.altitude()); |
1436 | seed = ::qHash(coordinates.latitude(), seed); |
1437 | |
1438 | return ::qHash(coordinates.longitude(), seed); |
1439 | } |
1440 | |
1441 | } // namespace Marble |