File indexing completed on 2025-01-05 03:59:33

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