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

KDE's Doxygen guidelines are available online.