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
QgsScaleCalculator(double dpi,QgsUnitTypes::DistanceUnit mapUnits)24 QgsScaleCalculator::QgsScaleCalculator( double dpi, QgsUnitTypes::DistanceUnit mapUnits )
25 : mDpi( dpi )
26 , mMapUnits( mapUnits )
27 {}
28
setDpi(double dpi)29 void QgsScaleCalculator::setDpi( double dpi )
30 {
31 mDpi = dpi;
32 }
dpi()33 double QgsScaleCalculator::dpi()
34 {
35 return mDpi;
36 }
37
setMapUnits(QgsUnitTypes::DistanceUnit mapUnits)38 void QgsScaleCalculator::setMapUnits( QgsUnitTypes::DistanceUnit mapUnits )
39 {
40 QgsDebugMsgLevel( QStringLiteral( "Map units set to %1" ).arg( QString::number( mapUnits ) ), 3 );
41 mMapUnits = mapUnits;
42 }
43
mapUnits() const44 QgsUnitTypes::DistanceUnit QgsScaleCalculator::mapUnits() const
45 {
46 QgsDebugMsgLevel( QStringLiteral( "Map units returned as %1" ).arg( QString::number( mMapUnits ) ), 4 );
47 return mMapUnits;
48 }
49
calculate(const QgsRectangle & mapExtent,double canvasWidth) const50 double QgsScaleCalculator::calculate( const QgsRectangle &mapExtent, double canvasWidth ) const
51 {
52 double conversionFactor = 0;
53 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 switch ( mMapUnits )
57 {
58 case QgsUnitTypes::DistanceMeters:
59 // convert meters to inches
60 conversionFactor = 39.3700787;
61 delta = mapExtent.xMaximum() - mapExtent.xMinimum();
62 break;
63 case QgsUnitTypes::DistanceFeet:
64 conversionFactor = 12.0;
65 delta = mapExtent.xMaximum() - mapExtent.xMinimum();
66 break;
67 case QgsUnitTypes::DistanceNauticalMiles:
68 // convert nautical miles to inches
69 conversionFactor = 72913.4;
70 delta = mapExtent.xMaximum() - mapExtent.xMinimum();
71 break;
72
73 default:
74 case QgsUnitTypes::DistanceDegrees:
75 // degrees require conversion to meters first
76 conversionFactor = 39.3700787;
77 delta = calculateGeographicDistance( mapExtent );
78 break;
79 }
80 if ( qgsDoubleNear( canvasWidth, 0. ) || qgsDoubleNear( mDpi, 0.0 ) )
81 {
82 QgsDebugMsg( QStringLiteral( "Can't calculate scale from the input values" ) );
83 return 0;
84 }
85 const double scale = ( delta * conversionFactor ) / ( static_cast< double >( canvasWidth ) / mDpi );
86 QgsDebugMsgLevel( QStringLiteral( "scale = %1 conversionFactor = %2" ).arg( scale ).arg( conversionFactor ), 4 );
87 return scale;
88 }
89
90
calculateGeographicDistance(const QgsRectangle & mapExtent) const91 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 const double lat = ( mapExtent.yMaximum() + mapExtent.yMinimum() ) * 0.5;
118 static const double RADS = ( 4.0 * std::atan( 1.0 ) ) / 180.0;
119 const double a = std::pow( std::cos( lat * RADS ), 2 );
120 const 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 const double radius = RA * ( 1.0 - E * E ) /
126 std::pow( 1.0 - E * E * std::sin( lat * RADS ) * std::sin( lat * RADS ), 1.5 );
127 const double meters = ( mapExtent.xMaximum() - mapExtent.xMinimum() ) / 180.0 * radius * c;
128
129 QgsDebugMsgLevel( "Distance across map extent (m): " + QString::number( meters ), 4 );
130
131 return meters;
132 }
133