Marble

GeoDataCoordinates.cpp
1// SPDX-License-Identifier: LGPL-2.1-or-later
2//
3// SPDX-FileCopyrightText: 2004-2007 Torsten Rahn <tackat@kde.org>
4// SPDX-FileCopyrightText: 2007-2008 Inge Wallin <ingwa@kde.org>
5// SPDX-FileCopyrightText: 2008 Patrick Spendrin <ps_ml@gmx.de>
6// SPDX-FileCopyrightText: 2011 Friedrich W. H. Kossebau <kossebau@kde.org>
7// SPDX-FileCopyrightText: 2011 Bernhard Beschow <bbeschow@cs.tu-berlin.de>
8// SPDX-FileCopyrightText: 2015 Alejandro Garcia Montoro <alejandro.garciamontoro@gmail.com>
9//
10
11#include "GeoDataCoordinates.h"
12#include "GeoDataCoordinates_p.h"
13#include "LonLatParser_p.h"
14
15#include <QDataStream>
16#include <QPointF>
17#include <qmath.h>
18
19#include "MarbleDebug.h"
20#include "MarbleMath.h"
21
22#include "Quaternion.h"
23
24namespace Marble
25{
26
27const qreal GeoDataCoordinatesPrivate::sm_semiMajorAxis = 6378137.0;
28const qreal GeoDataCoordinatesPrivate::sm_semiMinorAxis = 6356752.314;
29const qreal GeoDataCoordinatesPrivate::sm_eccentricitySquared = 6.69437999013e-03;
30const qreal GeoDataCoordinatesPrivate::sm_utmScaleFactor = 0.9996;
31GeoDataCoordinates::Notation GeoDataCoordinates::s_notation = GeoDataCoordinates::DMS;
32
33const GeoDataCoordinates GeoDataCoordinates::null = GeoDataCoordinates(0, 0, 0); // don't use default constructor!
34
35GeoDataCoordinates::GeoDataCoordinates(qreal _lon, qreal _lat, qreal _alt, GeoDataCoordinates::Unit unit, int _detail)
36 : d(new GeoDataCoordinatesPrivate(_lon, _lat, _alt, unit, _detail))
37{
38 d->ref.ref();
39}
40
41/* simply copy the d pointer
42 * it will be replaced in the detach function instead
43 */
45 : d(other.d)
46{
47 d->ref.ref();
48}
49
50/* simply copy null's d pointer
51 * it will be replaced in the detach function
52 */
54 : d(null.d)
55{
56 d->ref.ref();
57}
58
59/*
60 * only delete the private d pointer if the number of references is 0
61 * remember that all copies share the same d pointer!
62 */
63GeoDataCoordinates::~GeoDataCoordinates()
64{
65 delete d->m_q;
66 d->m_q = nullptr;
67
68 if (!d->ref.deref())
69 delete d;
70#ifdef DEBUG_GEODATA
71// mDebug() << "delete coordinates";
72#endif
73}
74
76{
77 return d != null.d;
78}
79
80/*
81 * if only one copy exists, return
82 * else make a new private d pointer object and assign the values of the current
83 * one to it
84 * at the end, if the number of references thus reaches 0 delete it
85 * this state shouldn't happen, but if it does, we have to clean up behind us.
86 */
87void GeoDataCoordinates::detach()
88{
89 delete d->m_q;
90 d->m_q = nullptr;
91
92 if (d->ref.fetchAndAddRelaxed(0) == 1) {
93 return;
94 }
95
96 auto new_d = new GeoDataCoordinatesPrivate(*d);
97
98 if (!d->ref.deref()) {
99 delete d;
100 }
101
102 d = new_d;
103 d->ref.ref();
104}
105
106/*
107 * call detach() at the start of all non-static, non-const functions
108 */
109void GeoDataCoordinates::set(qreal _lon, qreal _lat, qreal _alt, GeoDataCoordinates::Unit unit)
110{
111 detach();
112 d->m_altitude = _alt;
113 switch (unit) {
114 default:
115 case Radian:
116 d->m_lon = _lon;
117 d->m_lat = _lat;
118 break;
119 case Degree:
120 d->m_lon = _lon * DEG2RAD;
121 d->m_lat = _lat * DEG2RAD;
122 break;
123 }
124}
125
126/*
127 * call detach() at the start of all non-static, non-const functions
128 */
130{
131 detach();
132 switch (unit) {
133 default:
134 case Radian:
135 d->m_lon = _lon;
136 break;
137 case Degree:
138 d->m_lon = _lon * DEG2RAD;
139 break;
140 }
141}
142
143/*
144 * call detach() at the start of all non-static, non-const functions
145 */
147{
148 detach();
149 switch (unit) {
150 case Radian:
151 d->m_lat = _lat;
152 break;
153 case Degree:
154 d->m_lat = _lat * DEG2RAD;
155 break;
156 }
157}
158
159void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat, GeoDataCoordinates::Unit unit) const
160{
161 switch (unit) {
162 default:
163 case Radian:
164 lon = d->m_lon;
165 lat = d->m_lat;
166 break;
167 case Degree:
168 lon = d->m_lon * RAD2DEG;
169 lat = d->m_lat * RAD2DEG;
170 break;
171 }
172}
173
174void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat) const
175{
176 lon = d->m_lon;
177 lat = d->m_lat;
178}
179
180void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat, qreal &alt, GeoDataCoordinates::Unit unit) const
181{
182 geoCoordinates(lon, lat, unit);
183 alt = d->m_altitude;
184}
185
186void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat, qreal &alt) const
187{
188 lon = d->m_lon;
189 lat = d->m_lat;
190 alt = d->m_altitude;
191}
192
193qreal GeoDataCoordinates::longitude(GeoDataCoordinates::Unit unit) const
194{
195 switch (unit) {
196 default:
197 case Radian:
198 return d->m_lon;
199 case Degree:
200 return d->m_lon * RAD2DEG;
201 }
202}
203
204qreal GeoDataCoordinates::longitude() const
205{
206 return d->m_lon;
207}
208
209qreal GeoDataCoordinates::latitude(GeoDataCoordinates::Unit unit) const
210{
211 switch (unit) {
212 default:
213 case Radian:
214 return d->m_lat;
215 case Degree:
216 return d->m_lat * RAD2DEG;
217 }
218}
219
220qreal GeoDataCoordinates::latitude() const
221{
222 return d->m_lat;
223}
224
225// static
230
231// static
233{
234 s_notation = notation;
235}
236
237// static
239{
240 qreal halfCircle;
241 if (unit == GeoDataCoordinates::Radian) {
242 halfCircle = M_PI;
243 } else {
244 halfCircle = 180;
245 }
246
247 if (lon > halfCircle) {
248 int cycles = (int)((lon + halfCircle) / (2 * halfCircle));
249 return lon - (cycles * 2 * halfCircle);
250 }
251 if (lon < -halfCircle) {
252 int cycles = (int)((lon - halfCircle) / (2 * halfCircle));
253 return lon - (cycles * 2 * halfCircle);
254 }
255
256 return lon;
257}
258
259// static
261{
262 qreal halfCircle;
263 if (unit == GeoDataCoordinates::Radian) {
264 halfCircle = M_PI;
265 } else {
266 halfCircle = 180;
267 }
268
269 if (lat > (halfCircle / 2.0)) {
270 int cycles = (int)((lat + halfCircle) / (2 * halfCircle));
271 qreal temp;
272 if (cycles == 0) { // pi/2 < lat < pi
273 temp = halfCircle - lat;
274 } else {
275 temp = lat - (cycles * 2 * halfCircle);
276 }
277 if (temp > (halfCircle / 2.0)) {
278 return (halfCircle - temp);
279 }
280 if (temp < (-halfCircle / 2.0)) {
281 return (-halfCircle - temp);
282 }
283 return temp;
284 }
285 if (lat < (-halfCircle / 2.0)) {
286 int cycles = (int)((lat - halfCircle) / (2 * halfCircle));
287 qreal temp;
288 if (cycles == 0) {
289 temp = -halfCircle - lat;
290 } else {
291 temp = lat - (cycles * 2 * halfCircle);
292 }
293 if (temp > (+halfCircle / 2.0)) {
294 return (+halfCircle - temp);
295 }
296 if (temp < (-halfCircle / 2.0)) {
297 return (-halfCircle - temp);
298 }
299 return temp;
300 }
301 return lat;
302}
303
304// static
306{
307 qreal halfCircle;
308 if (unit == GeoDataCoordinates::Radian) {
309 halfCircle = M_PI;
310 } else {
311 halfCircle = 180;
312 }
313
314 if (lon > +halfCircle) {
315 int cycles = (int)((lon + halfCircle) / (2 * halfCircle));
316 lon = lon - (cycles * 2 * halfCircle);
317 }
318 if (lon < -halfCircle) {
319 int cycles = (int)((lon - halfCircle) / (2 * halfCircle));
320 lon = lon - (cycles * 2 * halfCircle);
321 }
322
323 if (lat > (+halfCircle / 2.0)) {
324 int cycles = (int)((lat + halfCircle) / (2 * halfCircle));
325 qreal temp;
326 if (cycles == 0) { // pi/2 < lat < pi
327 temp = halfCircle - lat;
328 } else {
329 temp = lat - (cycles * 2 * halfCircle);
330 }
331 if (temp > (+halfCircle / 2.0)) {
332 lat = +halfCircle - temp;
333 }
334 if (temp < (-halfCircle / 2.0)) {
335 lat = -halfCircle - temp;
336 }
337 lat = temp;
338 if (lon > 0) {
339 lon = -halfCircle + lon;
340 } else {
341 lon = halfCircle + lon;
342 }
343 }
344 if (lat < (-halfCircle / 2.0)) {
345 int cycles = (int)((lat - halfCircle) / (2 * halfCircle));
346 qreal temp;
347 if (cycles == 0) {
348 temp = -halfCircle - lat;
349 } else {
350 temp = lat - (cycles * 2 * halfCircle);
351 }
352 if (temp > (+halfCircle / 2.0)) {
353 lat = +halfCircle - temp;
354 }
355 if (temp < (-halfCircle / 2.0)) {
356 lat = -halfCircle - temp;
357 }
358 lat = temp;
359 if (lon > 0) {
360 lon = -halfCircle + lon;
361 } else {
362 lon = halfCircle + lon;
363 }
364 }
365 return;
366}
367
369{
370 LonLatParser parser;
371 successful = parser.parse(string);
372 if (successful) {
373 return {parser.lon(), parser.lat(), 0, GeoDataCoordinates::Degree};
374 } else {
375 return {};
376 }
377}
378
380{
381 return GeoDataCoordinates::toString(s_notation);
382}
383
385{
386 QString coordString;
387
388 if (notation == GeoDataCoordinates::UTM) {
389 int zoneNumber = GeoDataCoordinatesPrivate::lonLatToZone(d->m_lon, d->m_lat);
390
391 // Handle lack of UTM zone number in the poles
392 const QString zoneString = (zoneNumber > 0) ? QString::number(zoneNumber) : QString();
393
394 QString bandString = GeoDataCoordinatesPrivate::lonLatToLatitudeBand(d->m_lon, d->m_lat);
395
396 QString eastingString = QString::number(GeoDataCoordinatesPrivate::lonLatToEasting(d->m_lon, d->m_lat), 'f', 2);
397 QString northingString = QString::number(GeoDataCoordinatesPrivate::lonLatToNorthing(d->m_lon, d->m_lat), 'f', 2);
398
399 return QStringLiteral("%1%2 %3 m E, %4 m N").arg(zoneString, bandString, eastingString, northingString);
400 } else {
401 coordString = lonToString(d->m_lon, notation, Radian, precision) + QLatin1StringView(", ") + latToString(d->m_lat, notation, Radian, precision);
402 }
403
404 return coordString;
405}
406
408{
409 if (notation == GeoDataCoordinates::UTM) {
410 /**
411 * @FIXME: UTM needs lon + lat to know zone number and easting
412 * By now, this code returns the zone+easting of the point
413 * (lon, equator), but this can differ a lot at different locations
414 * See bug 347536 https://bugs.kde.org/show_bug.cgi?id=347536
415 */
416
417 qreal lonRad = (unit == Radian) ? lon : lon * DEG2RAD;
418
419 int zoneNumber = GeoDataCoordinatesPrivate::lonLatToZone(lonRad, 0);
420
421 // Handle lack of UTM zone number in the poles
422 QString result = (zoneNumber > 0) ? QString::number(zoneNumber) : QString();
423
424 if (precision > 0) {
425 QString eastingString = QString::number(GeoDataCoordinatesPrivate::lonLatToEasting(lonRad, 0), 'f', 2);
426 result += QStringLiteral(" %1 m E").arg(eastingString);
427 }
428
429 return result;
430 }
431
432 QString weString = (lon < 0) ? tr("W") : tr("E");
433
434 QString lonString;
435
436 qreal lonDegF = (unit == Degree) ? fabs(lon) : fabs((qreal)(lon)*RAD2DEG);
437
438 // Take care of -1 case
439 precision = (precision < 0) ? 5 : precision;
440
441 if (notation == DMS || notation == DM) {
442 int lonDeg = (int)lonDegF;
443 qreal lonMinF = 60 * (lonDegF - lonDeg);
444 int lonMin = (int)lonMinF;
445 qreal lonSecF = 60 * (lonMinF - lonMin);
446 int lonSec = (int)lonSecF;
447
448 // Adjustment for fuzziness (like 49.999999999999999999999)
449 if (precision == 0) {
450 lonDeg = qRound(lonDegF);
451 } else if (precision <= 2) {
452 lonMin = qRound(lonMinF);
453 } else if (precision <= 4 && notation == DMS) {
454 lonSec = qRound(lonSecF);
455 } else {
456 if (notation == DMS) {
457 lonSec = lonSecF = qRound(lonSecF * qPow(10, precision - 4)) / qPow(10, precision - 4);
458 } else {
459 lonMin = lonMinF = qRound(lonMinF * qPow(10, precision - 2)) / qPow(10, precision - 2);
460 }
461 }
462
463 if (lonSec > 59 && notation == DMS) {
464 lonSec = lonSecF = 0;
465 lonMin = lonMinF = lonMinF + 1;
466 }
467 if (lonMin > 59) {
468 lonMin = lonMinF = 0;
469 lonDeg = lonDegF = lonDegF + 1;
470 }
471
472 // Evaluate the string
473 lonString = QString::fromUtf8("%1\xc2\xb0").arg(lonDeg, 3, 10, QLatin1Char(' '));
474
475 if (precision == 0) {
476 return lonString + weString;
477 }
478
479 if (notation == DMS || precision < 3) {
480 lonString += QStringLiteral(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0'));
481 }
482
483 if (precision < 3) {
484 return lonString + weString;
485 }
486
487 if (notation == DMS) {
488 // Includes -1 case!
489 if (precision < 5) {
490 lonString += QStringLiteral(" %3\"").arg(lonSec, 2, 'f', 0, QLatin1Char('0'));
491 return lonString + weString;
492 }
493
494 lonString += QStringLiteral(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0'));
495 } else {
496 lonString += QStringLiteral(" %L3'").arg(lonMinF, precision + 1, 'f', precision - 2, QLatin1Char('0'));
497 }
498 } else if (notation == GeoDataCoordinates::Decimal) {
499 lonString = QString::fromUtf8("%L1\xc2\xb0").arg(lonDegF, 4 + precision, format, precision, QLatin1Char(' '));
500 } else if (notation == GeoDataCoordinates::Astro) {
501 if (lon < 0) {
502 lon += (unit == Degree) ? 360 : 2 * M_PI;
503 }
504
505 qreal lonHourF = (unit == Degree) ? fabs(lon / 15.0) : fabs((qreal)(lon / 15.0) * RAD2DEG);
506 int lonHour = (int)lonHourF;
507 qreal lonMinF = 60 * (lonHourF - lonHour);
508 int lonMin = (int)lonMinF;
509 qreal lonSecF = 60 * (lonMinF - lonMin);
510 int lonSec = (int)lonSecF;
511
512 // Adjustment for fuzziness (like 49.999999999999999999999)
513 if (precision == 0) {
514 lonHour = qRound(lonHourF);
515 } else if (precision <= 2) {
516 lonMin = qRound(lonMinF);
517 } else if (precision <= 4) {
518 lonSec = qRound(lonSecF);
519 } else {
520 lonSec = lonSecF = qRound(lonSecF * qPow(10, precision - 4)) / qPow(10, precision - 4);
521 }
522
523 if (lonSec > 59) {
524 lonSec = lonSecF = 0;
525 lonMin = lonMinF = lonMinF + 1;
526 }
527 if (lonMin > 59) {
528 lonMin = lonMinF = 0;
529 lonHour = lonHourF = lonHourF + 1;
530 }
531
532 // Evaluate the string
533 lonString = QString::fromUtf8("%1h").arg(lonHour, 3, 10, QLatin1Char(' '));
534
535 if (precision == 0) {
536 return lonString;
537 }
538
539 lonString += QStringLiteral(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0'));
540
541 if (precision < 3) {
542 return lonString;
543 }
544
545 // Includes -1 case!
546 if (precision < 5) {
547 lonString += QStringLiteral(" %3\"").arg(lonSec, 2, 'f', 0, QLatin1Char('0'));
548 return lonString;
549 }
550
551 lonString += QStringLiteral(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0'));
552 return lonString;
553 }
554
555 return lonString + weString;
556}
557
559{
560 return GeoDataCoordinates::lonToString(d->m_lon, s_notation);
561}
562
564{
565 if (notation == GeoDataCoordinates::UTM) {
566 /**
567 * @FIXME: UTM needs lon + lat to know latitude band and northing
568 * By now, this code returns the band+northing of the point
569 * (meridian, lat), but this can differ a lot at different locations
570 * See bug 347536 https://bugs.kde.org/show_bug.cgi?id=347536
571 */
572
573 qreal latRad = (unit == Radian) ? lat : lat * DEG2RAD;
574
575 QString result = GeoDataCoordinatesPrivate::lonLatToLatitudeBand(0, latRad);
576
577 if (precision > 0) {
578 QString northingString = QString::number(GeoDataCoordinatesPrivate::lonLatToNorthing(0, latRad), 'f', 2);
579 result += QStringLiteral(" %1 m N").arg(northingString);
580 }
581
582 return result;
583 }
584
585 QString pmString;
586 QString nsString;
587
588 if (notation == GeoDataCoordinates::Astro) {
589 pmString = (lat > 0) ? QStringLiteral("+") : QStringLiteral("-");
590 } else {
591 nsString = (lat > 0) ? tr("N") : tr("S");
592 }
593
594 QString latString;
595
596 qreal latDegF = (unit == Degree) ? fabs(lat) : fabs((qreal)(lat)*RAD2DEG);
597
598 // Take care of -1 case
599 precision = (precision < 0) ? 5 : precision;
600
601 if (notation == DMS || notation == DM || notation == Astro) {
602 int latDeg = (int)latDegF;
603 qreal latMinF = 60 * (latDegF - latDeg);
604 int latMin = (int)latMinF;
605 qreal latSecF = 60 * (latMinF - latMin);
606 int latSec = (int)latSecF;
607
608 // Adjustment for fuzziness (like 49.999999999999999999999)
609 if (precision == 0) {
610 latDeg = qRound(latDegF);
611 } else if (precision <= 2) {
612 latMin = qRound(latMinF);
613 } else if (precision <= 4 && notation == DMS) {
614 latSec = qRound(latSecF);
615 } else {
616 if (notation == DMS || notation == Astro) {
617 latSec = latSecF = qRound(latSecF * qPow(10, precision - 4)) / qPow(10, precision - 4);
618 } else {
619 latMin = latMinF = qRound(latMinF * qPow(10, precision - 2)) / qPow(10, precision - 2);
620 }
621 }
622
623 if (latSec > 59 && (notation == DMS || notation == Astro)) {
624 latSecF = 0;
625 latSec = latSecF;
626 latMin = latMin + 1;
627 }
628 if (latMin > 59) {
629 latMinF = 0;
630 latMin = latMinF;
631 latDeg = latDeg + 1;
632 }
633
634 // Evaluate the string
635 latString = QString::fromUtf8("%1\xc2\xb0").arg(latDeg, 3, 10, QLatin1Char(' '));
636
637 if (precision == 0) {
638 return pmString + latString + nsString;
639 }
640
641 if (notation == DMS || notation == Astro || precision < 3) {
642 latString += QStringLiteral(" %2\'").arg(latMin, 2, 10, QLatin1Char('0'));
643 }
644
645 if (precision < 3) {
646 return pmString + latString + nsString;
647 }
648
649 if (notation == DMS || notation == Astro) {
650 // Includes -1 case!
651 if (precision < 5) {
652 latString += QStringLiteral(" %3\"").arg(latSec, 2, 'f', 0, QLatin1Char('0'));
653 return latString + nsString;
654 }
655
656 latString += QStringLiteral(" %L3\"").arg(latSecF, precision - 1, 'f', precision - 4, QLatin1Char('0'));
657 } else {
658 latString += QStringLiteral(" %L3'").arg(latMinF, precision + 1, 'f', precision - 2, QLatin1Char('0'));
659 }
660 } else // notation = GeoDataCoordinates::Decimal
661 {
662 latString = QString::fromUtf8("%L1\xc2\xb0").arg(latDegF, 4 + precision, format, precision, QLatin1Char(' '));
663 }
664 return pmString + latString + nsString;
665}
666
668{
669 return GeoDataCoordinates::latToString(d->m_lat, s_notation);
670}
671
672bool GeoDataCoordinates::operator==(const GeoDataCoordinates &rhs) const
673{
674 return *d == *rhs.d;
675}
676
677bool GeoDataCoordinates::operator!=(const GeoDataCoordinates &rhs) const
678{
679 return *d != *rhs.d;
680}
681
682void GeoDataCoordinates::setAltitude(const qreal altitude)
683{
684 detach();
685 d->m_altitude = altitude;
686}
687
689{
690 return d->m_altitude;
691}
692
694{
695 return GeoDataCoordinatesPrivate::lonLatToZone(d->m_lon, d->m_lat);
696}
697
699{
700 return GeoDataCoordinatesPrivate::lonLatToEasting(d->m_lon, d->m_lat);
701}
702
704{
705 return GeoDataCoordinatesPrivate::lonLatToLatitudeBand(d->m_lon, d->m_lat);
706}
707
709{
710 return GeoDataCoordinatesPrivate::lonLatToNorthing(d->m_lon, d->m_lat);
711}
712
714{
715 return d->m_detail;
716}
717
719{
720 detach();
721 d->m_detail = detail;
722}
723
725{
726 const Quaternion quatAxis = Quaternion::fromEuler(-axis.latitude(), axis.longitude(), 0);
727 const Quaternion rotationAmount = Quaternion::fromEuler(0, 0, unit == Radian ? angle : angle * DEG2RAD);
728 const Quaternion resultAxis = quatAxis * rotationAmount * quatAxis.inverse();
729
730 return rotateAround(resultAxis);
731}
732
733GeoDataCoordinates GeoDataCoordinates::rotateAround(const Quaternion &rotAxis) const
734{
735 Quaternion rotatedQuat = quaternion();
736 rotatedQuat.rotateAroundAxis(rotAxis);
737 qreal rotatedLon, rotatedLat;
738 rotatedQuat.getSpherical(rotatedLon, rotatedLat);
739 return {rotatedLon, rotatedLat, altitude()};
740}
741
743{
744 if (type == FinalBearing) {
745 double const offset = unit == Degree ? 180.0 : M_PI;
746 return offset + other.bearing(*this, unit, InitialBearing);
747 }
748
749 qreal const delta = other.d->m_lon - d->m_lon;
750 double const bearing = atan2(sin(delta) * cos(other.d->m_lat), cos(d->m_lat) * sin(other.d->m_lat) - sin(d->m_lat) * cos(other.d->m_lat) * cos(delta));
751 return unit == Radian ? bearing : bearing * RAD2DEG;
752}
753
754GeoDataCoordinates GeoDataCoordinates::moveByBearing(qreal bearing, qreal distance) const
755{
756 qreal newLat = asin(sin(d->m_lat) * cos(distance) + cos(d->m_lat) * sin(distance) * cos(bearing));
757 qreal newLon = d->m_lon + atan2(sin(bearing) * sin(distance) * cos(d->m_lat), cos(distance) - sin(d->m_lat) * sin(newLat));
758
759 return {newLon, newLat};
760}
761
762const Quaternion &GeoDataCoordinates::quaternion() const
763{
764 if (d->m_q == nullptr) {
765 d->m_q = new Quaternion(Quaternion::fromSpherical(d->m_lon, d->m_lat));
766 }
767 return *d->m_q;
768}
769
771{
772 double const t = qBound(0.0, t_, 1.0);
773 Quaternion const quat = Quaternion::slerp(quaternion(), target.quaternion(), t);
774 qreal lon, lat;
775 quat.getSpherical(lon, lat);
776 double const alt = (1.0 - t) * d->m_altitude + t * target.d->m_altitude;
777 return {lon, lat, alt};
778}
779
781{
782 qreal lon = 0.0;
783 qreal lat = 0.0;
784
785 const Quaternion itpos = Quaternion::nlerp(quaternion(), target.quaternion(), t);
786 itpos.getSpherical(lon, lat);
787
788 const qreal altitude = 0.5 * (d->m_altitude + target.altitude());
789
790 return {lon, lat, altitude};
791}
792
794GeoDataCoordinates::interpolate(const GeoDataCoordinates &before, const GeoDataCoordinates &target, const GeoDataCoordinates &after, double t_) const
795{
796 double const t = qBound(0.0, t_, 1.0);
797 Quaternion const b1 = GeoDataCoordinatesPrivate::basePoint(before.quaternion(), quaternion(), target.quaternion());
798 Quaternion const a2 = GeoDataCoordinatesPrivate::basePoint(quaternion(), target.quaternion(), after.quaternion());
799 Quaternion const a = Quaternion::slerp(quaternion(), target.quaternion(), t);
800 Quaternion const b = Quaternion::slerp(b1, a2, t);
801 Quaternion c = Quaternion::slerp(a, b, 2 * t * (1.0 - t));
802 qreal lon, lat;
803 c.getSpherical(lon, lat);
804 // @todo spline interpolation of altitude?
805 double const alt = (1.0 - t) * d->m_altitude + t * target.d->m_altitude;
806 return {lon, lat, alt};
807}
808
810{
811 // Evaluate the most likely case first:
812 // The case where we haven't hit the pole and where our latitude is normalized
813 // to the range of 90 deg S ... 90 deg N
814 if (fabs((qreal)2.0 * d->m_lat) < M_PI) {
815 return false;
816 } else {
817 if (fabs((qreal)2.0 * d->m_lat) == M_PI) {
818 // Ok, we have hit a pole. Now let's check whether it's the one we've asked for:
819 if (pole == AnyPole) {
820 return true;
821 } else {
822 if (pole == NorthPole && 2.0 * d->m_lat == +M_PI) {
823 return true;
824 }
825 if (pole == SouthPole && 2.0 * d->m_lat == -M_PI) {
826 return true;
827 }
828 return false;
829 }
830 }
831 //
832 else {
833 // FIXME: Should we just normalize latitude and longitude and be done?
834 // While this might work well for persistent data it would create some
835 // possible overhead for temporary data, so this needs careful thinking.
836 mDebug() << "GeoDataCoordinates not normalized!";
837
838 // Only as a last resort we cover the unlikely case where
839 // the latitude is not normalized to the range of
840 // 90 deg S ... 90 deg N
841 if (fabs((qreal)2.0 * normalizeLat(d->m_lat)) < M_PI) {
842 return false;
843 } else {
844 // Ok, we have hit a pole. Now let's check whether it's the one we've asked for:
845 if (pole == AnyPole) {
846 return true;
847 } else {
848 if (pole == NorthPole && 2.0 * d->m_lat == +M_PI) {
849 return true;
850 }
851 if (pole == SouthPole && 2.0 * d->m_lat == -M_PI) {
852 return true;
853 }
854 return false;
855 }
856 }
857 }
858 }
859}
860
862{
863 qreal lon2, lat2;
864 other.geoCoordinates(lon2, lat2);
865
866 // FIXME: Take the altitude into account!
867
868 return distanceSphere(d->m_lon, d->m_lat, lon2, lat2);
869}
870
871GeoDataCoordinates &GeoDataCoordinates::operator=(const GeoDataCoordinates &other)
872{
873 qAtomicAssign(d, other.d);
874 return *this;
875}
876
878{
879 stream << d->m_lon;
880 stream << d->m_lat;
881 stream << d->m_altitude;
882}
883
885{
886 // call detach even though it shouldn't be needed - one never knows
887 detach();
888 stream >> d->m_lon;
889 stream >> d->m_lat;
890 stream >> d->m_altitude;
891}
892
893Quaternion GeoDataCoordinatesPrivate::basePoint(const Quaternion &q1, const Quaternion &q2, const Quaternion &q3)
894{
895 Quaternion const a = (q2.inverse() * q3).log();
896 Quaternion const b = (q2.inverse() * q1).log();
897 return q2 * ((a + b) * -0.25).exp();
898}
899
900qreal GeoDataCoordinatesPrivate::arcLengthOfMeridian(qreal phi)
901{
902 // Precalculate n
903 qreal const n = (GeoDataCoordinatesPrivate::sm_semiMajorAxis - GeoDataCoordinatesPrivate::sm_semiMinorAxis)
904 / (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis);
905
906 // Precalculate alpha
907 qreal const alpha = ((GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis) / 2.0)
908 * (1.0 + (qPow(n, 2.0) / 4.0) + (qPow(n, 4.0) / 64.0));
909
910 // Precalculate beta
911 qreal const beta = (-3.0 * n / 2.0) + (9.0 * qPow(n, 3.0) / 16.0) + (-3.0 * qPow(n, 5.0) / 32.0);
912
913 // Precalculate gamma
914 qreal const gamma = (15.0 * qPow(n, 2.0) / 16.0) + (-15.0 * qPow(n, 4.0) / 32.0);
915
916 // Precalculate delta
917 qreal const delta = (-35.0 * qPow(n, 3.0) / 48.0) + (105.0 * qPow(n, 5.0) / 256.0);
918
919 // Precalculate epsilon
920 qreal const epsilon = (315.0 * qPow(n, 4.0) / 512.0);
921
922 // Now calculate the sum of the series and return
923 qreal const result = alpha * (phi + (beta * qSin(2.0 * phi)) + (gamma * qSin(4.0 * phi)) + (delta * qSin(6.0 * phi)) + (epsilon * qSin(8.0 * phi)));
924
925 return result;
926}
927
928qreal GeoDataCoordinatesPrivate::centralMeridianUTM(qreal zone)
929{
930 return DEG2RAD * (-183.0 + (zone * 6.0));
931}
932
933qreal GeoDataCoordinatesPrivate::footpointLatitude(qreal northing)
934{
935 // Precalculate n (Eq. 10.18)
936 qreal const n = (GeoDataCoordinatesPrivate::sm_semiMajorAxis - GeoDataCoordinatesPrivate::sm_semiMinorAxis)
937 / (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis);
938
939 // Precalculate alpha (Eq. 10.22)
940 // (Same as alpha in Eq. 10.17)
941 qreal const alpha =
942 ((GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis) / 2.0) * (1 + (qPow(n, 2.0) / 4) + (qPow(n, 4.0) / 64));
943
944 // Precalculate y (Eq. 10.23)
945 qreal const y = northing / alpha;
946
947 // Precalculate beta (Eq. 10.22)
948 qreal const beta = (3.0 * n / 2.0) + (-27.0 * qPow(n, 3.0) / 32.0) + (269.0 * qPow(n, 5.0) / 512.0);
949
950 // Precalculate gamma (Eq. 10.22)
951 qreal const gamma = (21.0 * qPow(n, 2.0) / 16.0) + (-55.0 * qPow(n, 4.0) / 32.0);
952
953 // Precalculate delta (Eq. 10.22)
954 qreal const delta = (151.0 * qPow(n, 3.0) / 96.0) + (-417.0 * qPow(n, 5.0) / 128.0);
955
956 // Precalculate epsilon (Eq. 10.22)
957 qreal const epsilon = (1097.0 * qPow(n, 4.0) / 512.0);
958
959 // Now calculate the sum of the series (Eq. 10.21)
960 qreal const result = y + (beta * qSin(2.0 * y)) + (gamma * qSin(4.0 * y)) + (delta * qSin(6.0 * y)) + (epsilon * qSin(8.0 * y));
961
962 return result;
963}
964
965QPointF GeoDataCoordinatesPrivate::mapLonLatToXY(qreal lambda, qreal phi, qreal lambda0)
966{
967 // Equation (10.15)
968
969 // Precalculate second numerical eccentricity
970 const qreal ep2 = (qPow(GeoDataCoordinatesPrivate::sm_semiMajorAxis, 2.0) - qPow(GeoDataCoordinatesPrivate::sm_semiMinorAxis, 2.0))
971 / qPow(GeoDataCoordinatesPrivate::sm_semiMinorAxis, 2.0);
972
973 // Precalculate the square of nu, just an auxilar quantity
974 const qreal nu2 = ep2 * qPow(qCos(phi), 2.0);
975
976 // Precalculate the radius of curvature in prime vertical
977 const qreal N = qPow(GeoDataCoordinatesPrivate::sm_semiMajorAxis, 2.0) / (GeoDataCoordinatesPrivate::sm_semiMinorAxis * qSqrt(1 + nu2));
978
979 // Precalculate the tangent of phi and its square
980 const qreal t = qTan(phi);
981 const qreal t2 = t * t;
982
983 // Precalculate longitude difference
984 const qreal l = lambda - lambda0;
985
986 /*
987 * Precalculate coefficients for l**n in the equations below
988 * so a normal human being can read the expressions for easting
989 * and northing
990 * -- l**1 and l**2 have coefficients of 1.0
991 *
992 * The actual used coefficients starts at coef[1], just to
993 * follow the meaningful nomenclature in equation 10.15
994 * (coef[n] corresponds to qPow(l,n) factor)
995 */
996
997 const qreal coef1 = 1;
998 const qreal coef2 = 1;
999 const qreal coef3 = 1.0 - t2 + nu2;
1000 const qreal coef4 = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2);
1001 const qreal coef5 = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 - 58.0 * t2 * nu2;
1002 const qreal coef6 = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 - 330.0 * t2 * nu2;
1003 const qreal coef7 = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2);
1004 const qreal coef8 = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
1005
1006 // Calculate easting (x)
1007 const qreal easting = N * qCos(phi) * coef1 * l + (N / 6.0 * qPow(qCos(phi), 3.0) * coef3 * qPow(l, 3.0))
1008 + (N / 120.0 * qPow(qCos(phi), 5.0) * coef5 * qPow(l, 5.0)) + (N / 5040.0 * qPow(qCos(phi), 7.0) * coef7 * qPow(l, 7.0));
1009
1010 // Calculate northing (y)
1011 const qreal northing = arcLengthOfMeridian(phi) + (t / 2.0 * N * qPow(qCos(phi), 2.0) * coef2 * qPow(l, 2.0))
1012 + (t / 24.0 * N * qPow(qCos(phi), 4.0) * coef4 * qPow(l, 4.0)) + (t / 720.0 * N * qPow(qCos(phi), 6.0) * coef6 * qPow(l, 6.0))
1013 + (t / 40320.0 * N * qPow(qCos(phi), 8.0) * coef8 * qPow(l, 8.0));
1014
1015 return {easting, northing};
1016}
1017
1018int GeoDataCoordinatesPrivate::lonLatToZone(qreal lon, qreal lat)
1019{
1020 // Converts lon and lat to degrees
1021 qreal lonDeg = lon * RAD2DEG;
1022 qreal latDeg = lat * RAD2DEG;
1023
1024 /* Round the value of the longitude when the distance to the nearest integer
1025 * is less than 0.0000001. This avoids fuzzy values such as -114.0000000001, which
1026 * can produce a misbehaviour when calculating the zone associated at the borders
1027 * of the zone intervals (for example, the interval [-114, -108[ is associated with
1028 * zone number 12; if the following rounding is not done, the value returned by
1029 * lonLatToZone(114,0) is 11 instead of 12, as the function actually receives
1030 * -114.0000000001, which is in the interval [-120,-114[, associated to zone 11
1031 */
1032 qreal precision = 0.0000001;
1033
1034 if (qAbs(lonDeg - qFloor(lonDeg)) < precision || qAbs(lonDeg - qCeil(lonDeg)) < precision) {
1035 lonDeg = qRound(lonDeg);
1036 }
1037
1038 // There is no numbering associated to the poles, special value 0 is returned.
1039 if (latDeg < -80 || latDeg > 84) {
1040 return 0;
1041 }
1042
1043 // Obtains the zone number handling all the so called "exceptions"
1044 // See problem: https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions
1045 // See solution: https://gis.stackexchange.com/questions/13291/computing-utm-zone-from-lat-long-point
1046
1047 // General
1048 int zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1;
1049
1050 // Southwest Norway
1051 if (latDeg >= 56 && latDeg < 64 && lonDeg >= 3 && lonDeg < 12) {
1052 zoneNumber = 32;
1053 }
1054
1055 // Svalbard
1056 if (latDeg >= 72 && latDeg < 84) {
1057 if (lonDeg >= 0 && lonDeg < 9) {
1058 zoneNumber = 31;
1059 } else if (lonDeg >= 9 && lonDeg < 21) {
1060 zoneNumber = 33;
1061 } else if (lonDeg >= 21 && lonDeg < 33) {
1062 zoneNumber = 35;
1063 } else if (lonDeg >= 33 && lonDeg < 42) {
1064 zoneNumber = 37;
1065 }
1066 }
1067
1068 return zoneNumber;
1069}
1070
1071qreal GeoDataCoordinatesPrivate::lonLatToEasting(qreal lon, qreal lat)
1072{
1073 int zoneNumber = lonLatToZone(lon, lat);
1074
1075 if (zoneNumber == 0) {
1076 qreal lonDeg = lon * RAD2DEG;
1077 zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1;
1078 }
1079
1080 QPointF coordinates = GeoDataCoordinatesPrivate::mapLonLatToXY(lon, lat, GeoDataCoordinatesPrivate::centralMeridianUTM(zoneNumber));
1081
1082 // Adjust easting and northing for UTM system
1083 qreal easting = coordinates.x() * GeoDataCoordinatesPrivate::sm_utmScaleFactor + 500000.0;
1084
1085 return easting;
1086}
1087
1088QString GeoDataCoordinatesPrivate::lonLatToLatitudeBand(qreal lon, qreal lat)
1089{
1090 // Obtains the latitude bands handling all the so called "exceptions"
1091
1092 // Converts lon and lat to degrees
1093 qreal lonDeg = lon * RAD2DEG;
1094 qreal latDeg = lat * RAD2DEG;
1095
1096 // Regular latitude bands between 80 S and 80 N (that is, between 10 and 170 in the [0,180] interval)
1097 int bandLetterIndex = 24; // Avoids "may be used uninitialized" warning
1098
1099 if (latDeg < -80) {
1100 // South pole (A for zones 1-30, B for zones 31-60)
1101 bandLetterIndex = ((lonDeg + 180) < 6 * 31) ? 0 : 1;
1102 } else if (latDeg >= -80 && latDeg <= 80) {
1103 // General (+2 because the general lettering starts in C)
1104 bandLetterIndex = static_cast<int>((latDeg + 80.0) / 8.0) + 2;
1105 } else if (latDeg >= 80 && latDeg < 84) {
1106 // Band X is extended 4 more degrees
1107 bandLetterIndex = 21;
1108 } else if (latDeg >= 84) {
1109 // North pole (Y for zones 1-30, Z for zones 31-60)
1110 bandLetterIndex = ((lonDeg + 180) < 6 * 31) ? 22 : 23;
1111 }
1112
1113 return QStringLiteral("ABCDEFGHJKLMNPQRSTUVWXYZ?").at(bandLetterIndex);
1114}
1115
1116qreal GeoDataCoordinatesPrivate::lonLatToNorthing(qreal lon, qreal lat)
1117{
1118 int zoneNumber = lonLatToZone(lon, lat);
1119
1120 if (zoneNumber == 0) {
1121 qreal lonDeg = lon * RAD2DEG;
1122 zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1;
1123 }
1124
1125 QPointF coordinates = GeoDataCoordinatesPrivate::mapLonLatToXY(lon, lat, GeoDataCoordinatesPrivate::centralMeridianUTM(zoneNumber));
1126
1127 qreal northing = coordinates.y() * GeoDataCoordinatesPrivate::sm_utmScaleFactor;
1128
1129 if (northing < 0.0) {
1130 northing += 10000000.0;
1131 }
1132
1133 return northing;
1134}
1135
1136size_t qHash(const GeoDataCoordinates &coordinates)
1137{
1138 uint seed = ::qHash(coordinates.altitude());
1139 seed = ::qHash(coordinates.latitude(), seed);
1140
1141 return ::qHash(coordinates.longitude(), seed);
1142}
1143
1144}
A 3d point representation.
bool isPole(Pole=AnyPole) const
return whether our coordinates represent a pole This method can be used to check whether the coordina...
void unpack(QDataStream &stream)
Unserialize the contents of the feature from stream.
void set(qreal lon, qreal lat, qreal alt=0, GeoDataCoordinates::Unit unit=GeoDataCoordinates::Radian)
(re)set the coordinates in a GeoDataCoordinates object
static void normalizeLonLat(qreal &lon, qreal &lat, GeoDataCoordinates::Unit=GeoDataCoordinates::Radian)
normalize both longitude and latitude at the same time This method normalizes both latitude and longi...
static qreal normalizeLon(qreal lon, GeoDataCoordinates::Unit=GeoDataCoordinates::Radian)
normalize the longitude to always be -M_PI <= lon <= +M_PI (Radian).
GeoDataCoordinates moveByBearing(qreal bearing, qreal distance) const
Returns the coordinates of the resulting point after moving this point according to the distance and ...
void setLatitude(qreal lat, GeoDataCoordinates::Unit unit=GeoDataCoordinates::Radian)
set the longitude in a GeoDataCoordinates object
qreal utmNorthing() const
retrieves the UTM northing of the GeoDataCoordinates object, in meters
GeoDataCoordinates()
constructs an invalid instance
qreal bearing(const GeoDataCoordinates &other, Unit unit=Radian, BearingType type=InitialBearing) const
Returns the bearing (true bearing, the angle between the line defined by this point and the other and...
QString latToString() const
return a string representation of latitude of the coordinate convenience function that uses the defau...
void pack(QDataStream &stream) const
Serialize the contents of the feature to stream.
GeoDataCoordinates nlerp(const GeoDataCoordinates &target, double t) const
nlerp (normalized linear interpolation) between this coordinates and the given target coordinates
void setAltitude(const qreal altitude)
set the altitude of the Point in meters
GeoDataCoordinates rotateAround(const GeoDataCoordinates &axis, qreal angle, Unit unit=Radian) const
Rotates one coordinate around another.
Notation
enum used to specify the notation / numerical system
@ Astro
< "RA and DEC" notation (used for astronomical sky coordinates)
@ DM
"Sexagesimal DM" notation (base-60)
@ DMS
"Sexagesimal DMS" notation (base-60)
@ Decimal
"Decimal" notation (base-10)
BearingType
The BearingType enum specifies where to measure the bearing along great circle arcs.
qreal altitude() const
return the altitude of the Point in meters
static void setDefaultNotation(GeoDataCoordinates::Notation notation)
set the Notation of the string representation
void setDetail(quint8 detail)
set the detail flag
static GeoDataCoordinates::Notation defaultNotation()
return Notation of string representation
qreal utmEasting() const
retrieves the UTM easting of the GeoDataCoordinates object, in meters.
void setLongitude(qreal lon, GeoDataCoordinates::Unit unit=GeoDataCoordinates::Radian)
set the longitude in a GeoDataCoordinates object
Unit
enum used constructor to specify the units used
QString lonToString() const
return a string representation of longitude of the coordinate convenience function that uses the defa...
int utmZone() const
retrieves the UTM zone of the GeoDataCoordinates object.
static GeoDataCoordinates fromString(const QString &string, bool &successful)
try to parse the string into a coordinate pair
GeoDataCoordinates interpolate(const GeoDataCoordinates &target, double t) const
slerp (spherical linear) interpolation between this coordinate and the given target coordinate
static qreal normalizeLat(qreal lat, GeoDataCoordinates::Unit=GeoDataCoordinates::Radian)
normalize latitude to always be in -M_PI / 2.
qreal sphericalDistanceTo(const GeoDataCoordinates &other) const
This method calculates the shortest distance between two points on a sphere.
quint8 detail() const
return the detail flag detail range: 0 for most important points, 5 for least important
QString utmLatitudeBand() const
retrieves the UTM latitude band of the GeoDataCoordinates object
void geoCoordinates(qreal &lon, qreal &lat, GeoDataCoordinates::Unit unit) const
use this function to get the longitude and latitude with one call - use the unit parameter to switch ...
QString toString() const
return a string representation of the coordinate this is a convenience function which uses the defaul...
const Quaternion & quaternion() const
return a Quaternion with the used coordinates
KTEXTEDITOR_EXPORT size_t qHash(KTextEditor::Cursor cursor, size_t seed=0) noexcept
Binds a QML item to a specific geodetic location in screen coordinates.
@ SouthPole
Only South Pole.
@ NorthPole
Only North Pole.
@ AnyPole
Any pole.
qreal distanceSphere(qreal lon1, qreal lat1, qreal lon2, qreal lat2)
This method calculates the shortest distance between two points on a sphere.
Definition MarbleMath.h:45
qreal x() const const
qreal y() const const
QString arg(Args &&... args) const const
QString fromUtf8(QByteArrayView str)
QString number(double n, char format, int precision)
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri Dec 6 2024 12:08:52 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.