Branch data     Line data    Source code

       1                 :            : /***************************************************************************
       2                 :            :                               qgsscalecalculator.h
       3                 :            :                  Calculates scale based on map extent and units
       4                 :            :                               -------------------
       5                 :            :   begin                : May 18, 2004
       6                 :            :   copyright            : (C) 2004 by Gary E.Sherman
       7                 :            :   email                : sherman at mrcc.com
       8                 :            :  ***************************************************************************/
       9                 :            : 
      10                 :            : /***************************************************************************
      11                 :            :  *                                                                         *
      12                 :            :  *   This program is free software; you can redistribute it and/or modify  *
      13                 :            :  *   it under the terms of the GNU General Public License as published by  *
      14                 :            :  *   the Free Software Foundation; either version 2 of the License, or     *
      15                 :            :  *   (at your option) any later version.                                   *
      16                 :            :  *                                                                         *
      17                 :            :  ***************************************************************************/
      18                 :            : 
      19                 :            : #include <cmath>
      20                 :            : #include "qgslogger.h"
      21                 :            : #include "qgsrectangle.h"
      22                 :            : #include "qgsscalecalculator.h"
      23                 :            : 
      24                 :          6 : QgsScaleCalculator::QgsScaleCalculator( double dpi, QgsUnitTypes::DistanceUnit mapUnits )
      25                 :          6 :   : mDpi( dpi )
      26                 :          6 :   , mMapUnits( mapUnits )
      27                 :          6 : {}
      28                 :            : 
      29                 :          0 : void QgsScaleCalculator::setDpi( double dpi )
      30                 :            : {
      31                 :          0 :   mDpi = dpi;
      32                 :          0 : }
      33                 :          0 : double QgsScaleCalculator::dpi()
      34                 :            : {
      35                 :          0 :   return mDpi;
      36                 :            : }
      37                 :            : 
      38                 :          6 : void QgsScaleCalculator::setMapUnits( QgsUnitTypes::DistanceUnit mapUnits )
      39                 :            : {
      40                 :          6 :   QgsDebugMsgLevel( QStringLiteral( "Map units set to %1" ).arg( QString::number( mapUnits ) ), 3 );
      41                 :          6 :   mMapUnits = mapUnits;
      42                 :          6 : }
      43                 :            : 
      44                 :          0 : QgsUnitTypes::DistanceUnit QgsScaleCalculator::mapUnits() const
      45                 :            : {
      46                 :          0 :   QgsDebugMsgLevel( QStringLiteral( "Map units returned as %1" ).arg( QString::number( mMapUnits ) ), 4 );
      47                 :          0 :   return mMapUnits;
      48                 :            : }
      49                 :            : 
      50                 :          0 : double QgsScaleCalculator::calculate( const QgsRectangle &mapExtent, double canvasWidth )  const
      51                 :            : {
      52                 :          0 :   double conversionFactor = 0;
      53                 :          0 :   double delta = 0;
      54                 :            :   // calculation is based on the map units and extent, the dpi of the
      55                 :            :   // users display, and the canvas width
      56                 :          0 :   switch ( mMapUnits )
      57                 :            :   {
      58                 :            :     case QgsUnitTypes::DistanceMeters:
      59                 :            :       // convert meters to inches
      60                 :          0 :       conversionFactor = 39.3700787;
      61                 :          0 :       delta = mapExtent.xMaximum() - mapExtent.xMinimum();
      62                 :          0 :       break;
      63                 :            :     case QgsUnitTypes::DistanceFeet:
      64                 :          0 :       conversionFactor = 12.0;
      65                 :          0 :       delta = mapExtent.xMaximum() - mapExtent.xMinimum();
      66                 :          0 :       break;
      67                 :            :     case QgsUnitTypes::DistanceNauticalMiles:
      68                 :            :       // convert nautical miles to inches
      69                 :          0 :       conversionFactor = 72913.4;
      70                 :          0 :       delta = mapExtent.xMaximum() - mapExtent.xMinimum();
      71                 :          0 :       break;
      72                 :            : 
      73                 :            :     default:
      74                 :            :     case QgsUnitTypes::DistanceDegrees:
      75                 :            :       // degrees require conversion to meters first
      76                 :          0 :       conversionFactor = 39.3700787;
      77                 :          0 :       delta = calculateGeographicDistance( mapExtent );
      78                 :          0 :       break;
      79                 :            :   }
      80                 :          0 :   if ( qgsDoubleNear( canvasWidth, 0. ) || qgsDoubleNear( mDpi, 0.0 ) )
      81                 :            :   {
      82                 :          0 :     QgsDebugMsg( QStringLiteral( "Can't calculate scale from the input values" ) );
      83                 :          0 :     return 0;
      84                 :            :   }
      85                 :          0 :   double scale = ( delta * conversionFactor ) / ( static_cast< double >( canvasWidth ) / mDpi );
      86                 :          0 :   QgsDebugMsgLevel( QStringLiteral( "scale = %1 conversionFactor = %2" ).arg( scale ).arg( conversionFactor ), 4 );
      87                 :          0 :   return scale;
      88                 :          0 : }
      89                 :            : 
      90                 :            : 
      91                 :          0 : double QgsScaleCalculator::calculateGeographicDistance( const QgsRectangle &mapExtent ) const
      92                 :            : {
      93                 :            :   // need to calculate the x distance in meters
      94                 :            :   // We'll use the middle latitude for the calculation
      95                 :            :   // Note this is an approximation (although very close) but calculating scale
      96                 :            :   // for geographic data over large extents is quasi-meaningless
      97                 :            : 
      98                 :            :   // The distance between two points on a sphere can be estimated
      99                 :            :   // using the Haversine formula. This gives the shortest distance
     100                 :            :   // between two points on the sphere. However, what we're after is
     101                 :            :   // the distance from the left of the given extent and the right of
     102                 :            :   // it. This is not necessarily the shortest distance between two
     103                 :            :   // points on a sphere.
     104                 :            :   //
     105                 :            :   // The code below uses the Haversine formula, but with some changes
     106                 :            :   // to cope with the above problem, and also to deal with the extent
     107                 :            :   // possibly extending beyond +/-180 degrees:
     108                 :            :   //
     109                 :            :   // - Use the Halversine formula to calculate the distance from -90 to
     110                 :            :   //   +90 degrees at the mean latitude.
     111                 :            :   // - Scale this distance by the number of degrees between
     112                 :            :   //   mapExtent.xMinimum() and mapExtent.xMaximum();
     113                 :            :   // - For a slight improvemnt, allow for the ellipsoid shape of earth.
     114                 :            : 
     115                 :            : 
     116                 :            :   // For a longitude change of 180 degrees
     117                 :          0 :   double lat = ( mapExtent.yMaximum() + mapExtent.yMinimum() ) * 0.5;
     118                 :          0 :   static const double RADS = ( 4.0 * std::atan( 1.0 ) ) / 180.0;
     119                 :          0 :   double a = std::pow( std::cos( lat * RADS ), 2 );
     120                 :          0 :   double c = 2.0 * std::atan2( std::sqrt( a ), std::sqrt( 1.0 - a ) );
     121                 :            :   static const double RA = 6378000; // [m]
     122                 :            :   // The eccentricity. This comes from sqrt(1.0 - rb*rb/(ra*ra)) with rb set
     123                 :            :   // to 6357000 m.
     124                 :            :   static const double E = 0.0810820288;
     125                 :          0 :   double radius = RA * ( 1.0 - E * E ) /
     126                 :          0 :                   std::pow( 1.0 - E * E * std::sin( lat * RADS ) * std::sin( lat * RADS ), 1.5 );
     127                 :          0 :   double meters = ( mapExtent.xMaximum() - mapExtent.xMinimum() ) / 180.0 * radius * c;
     128                 :            : 
     129                 :          0 :   QgsDebugMsgLevel( "Distance across map extent (m): " + QString::number( meters ), 4 );
     130                 :            : 
     131                 :          0 :   return meters;
     132                 :            : }