File indexing completed on 2024-05-05 03:49:52

0001 // SPDX-License-Identifier: LGPL-2.1-or-later
0002 //
0003 // SPDX-FileCopyrightText: 2014 Torsten Rahn <rahn@kde.org>
0004 //
0005 
0006 // Local
0007 #include "LambertAzimuthalProjection.h"
0008 #include "AbstractProjection_p.h"
0009 
0010 #include "MarbleDebug.h"
0011 
0012 // Marble
0013 #include "ViewportParams.h"
0014 #include "GeoDataPoint.h"
0015 #include "GeoDataLineString.h"
0016 #include "GeoDataCoordinates.h"
0017 #include "MarbleGlobal.h"
0018 #include "AzimuthalProjection_p.h"
0019 
0020 #include <QIcon>
0021 #include <qmath.h>
0022 
0023 #define SAFE_DISTANCE
0024 
0025 namespace Marble
0026 {
0027 
0028 class LambertAzimuthalProjectionPrivate : public AzimuthalProjectionPrivate
0029 {
0030   public:
0031     explicit LambertAzimuthalProjectionPrivate( LambertAzimuthalProjection * parent );
0032 
0033     Q_DECLARE_PUBLIC( LambertAzimuthalProjection )
0034 };
0035 
0036 LambertAzimuthalProjection::LambertAzimuthalProjection()
0037     : AzimuthalProjection( new LambertAzimuthalProjectionPrivate( this ) )
0038 {
0039     setMinLat( minValidLat() );
0040     setMaxLat( maxValidLat() );
0041 }
0042 
0043 LambertAzimuthalProjection::LambertAzimuthalProjection( LambertAzimuthalProjectionPrivate *dd )
0044         : AzimuthalProjection( dd )
0045 {
0046     setMinLat( minValidLat() );
0047     setMaxLat( maxValidLat() );
0048 }
0049 
0050 LambertAzimuthalProjection::~LambertAzimuthalProjection()
0051 {
0052 }
0053 
0054 
0055 LambertAzimuthalProjectionPrivate::LambertAzimuthalProjectionPrivate( LambertAzimuthalProjection * parent )
0056         : AzimuthalProjectionPrivate( parent )
0057 {
0058 }
0059 
0060 QString LambertAzimuthalProjection::name() const
0061 {
0062     return QObject::tr( "Lambert Azimuthal Equal-Area" );
0063 }
0064 
0065 QString LambertAzimuthalProjection::description() const
0066 {
0067     return QObject::tr( "<p><b>Lambert Azimuthal Equal-Area Projection</b></p><p>Applications: Used in structural geology to plot directional data.</p>" );
0068 }
0069 
0070 QIcon LambertAzimuthalProjection::icon() const
0071 {
0072     return QIcon(QStringLiteral(":/icons/map-globe.png"));
0073 }
0074 
0075 qreal LambertAzimuthalProjection::clippingRadius() const
0076 {
0077     return 1;
0078 }
0079 
0080 bool LambertAzimuthalProjection::screenCoordinates( const GeoDataCoordinates &coordinates,
0081                                              const ViewportParams *viewport,
0082                                              qreal &x, qreal &y, bool &globeHidesPoint ) const
0083 {
0084     const qreal lambda = coordinates.longitude();
0085     const qreal phi = coordinates.latitude();
0086     const qreal lambdaPrime = viewport->centerLongitude();
0087     const qreal phi1 = viewport->centerLatitude();
0088 
0089     qreal cosC = qSin( phi1 ) * qSin( phi ) + qCos( phi1 ) * qCos( phi ) * qCos( lambda - lambdaPrime );
0090     // Prevent division by zero
0091     if (cosC <= 0) {
0092         globeHidesPoint = true;
0093         return false;
0094     }
0095 
0096     qreal k = qSqrt(2 / (1 + cosC));
0097 
0098     // Let (x, y) be the position on the screen of the placemark..
0099     x = ( qCos( phi ) * qSin( lambda - lambdaPrime ) ) * k;
0100     y = ( qCos( phi1 ) * qSin( phi ) - qSin( phi1 ) * qCos( phi ) * qCos( lambda - lambdaPrime ) ) * k;
0101 
0102     x *= viewport->radius() / qSqrt(2);
0103     y *= viewport->radius() / qSqrt(2);
0104 
0105     const qint64  radius  = clippingRadius() * viewport->radius();
0106 
0107     if (x*x + y*y > radius * radius) {
0108         globeHidesPoint = true;
0109         return false;
0110     }
0111 
0112     globeHidesPoint = false;
0113 
0114     x += viewport->width() / 2;
0115     y = viewport->height() / 2 - y;
0116 
0117     // Skip placemarks that are outside the screen area
0118     return !(x < 0 || x >= viewport->width() || y < 0 || y >= viewport->height());
0119 }
0120 
0121 bool LambertAzimuthalProjection::screenCoordinates( const GeoDataCoordinates &coordinates,
0122                                              const ViewportParams *viewport,
0123                                              qreal *x, qreal &y,
0124                                              int &pointRepeatNum,
0125                                              const QSizeF& size,
0126                                              bool &globeHidesPoint ) const
0127 {
0128     pointRepeatNum = 0;
0129     globeHidesPoint = false;
0130 
0131     bool visible = screenCoordinates( coordinates, viewport, *x, y, globeHidesPoint );
0132 
0133     // Skip placemarks that are outside the screen area
0134     if ( *x + size.width() / 2.0 < 0.0 || *x >= viewport->width() + size.width() / 2.0
0135          || y + size.height() / 2.0 < 0.0 || y >= viewport->height() + size.height() / 2.0 )
0136     {
0137         return false;
0138     }
0139 
0140     // This projection doesn't have any repetitions,
0141     // so the number of screen points referring to the geopoint is one.
0142     pointRepeatNum = 1;
0143     return visible;
0144 }
0145 
0146 
0147 bool LambertAzimuthalProjection::geoCoordinates( const int x, const int y,
0148                                           const ViewportParams *viewport,
0149                                           qreal& lon, qreal& lat,
0150                                           GeoDataCoordinates::Unit unit ) const
0151 {
0152     const qint64  radius  = viewport->radius();
0153     // Calculate how many degrees are being represented per pixel.
0154     const qreal centerLon = viewport->centerLongitude();
0155     const qreal centerLat = viewport->centerLatitude();
0156     const qreal rx = ( - viewport->width()  / 2 + x );
0157     const qreal ry = (   viewport->height() / 2 - y );
0158     const qreal p = qMax( qSqrt( rx*rx + ry*ry ), qreal(0.0001) ); // ensure we don't divide by zero
0159 
0160     // exclude area too far away from map, as it may cause undefined arcsin result
0161     if ( p > (qSqrt(2) * radius) ) {
0162         return false;
0163     }
0164     const qreal c = 2 * qAsin( p / (qSqrt(2) * radius)  );
0165     const qreal sinc = qSin(c);
0166 
0167     lon = centerLon + qAtan2( rx*sinc , ( p*qCos( centerLat )*qCos( c ) - ry*qSin( centerLat )*sinc  ) );
0168 
0169     while ( lon < -M_PI ) lon += 2 * M_PI;
0170     while ( lon >  M_PI ) lon -= 2 * M_PI;
0171 
0172     lat = qAsin( qCos(c)*qSin(centerLat) + (ry*sinc*qCos(centerLat))/p );
0173 
0174     if ( unit == GeoDataCoordinates::Degree ) {
0175         lon *= RAD2DEG;
0176         lat *= RAD2DEG;
0177     }
0178 
0179     return true;
0180 }
0181 
0182 }