QGIS API Documentation  3.9.0-Master (224899f119)
qgsscalecalculator.cpp
Go to the documentation of this file.
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 
25  : mDpi( dpi )
26  , mMapUnits( mapUnits )
27 {}
28 
30 {
31  mDpi = dpi;
32 }
34 {
35  return mDpi;
36 }
37 
39 {
40  QgsDebugMsgLevel( QStringLiteral( "Map units set to %1" ).arg( QString::number( mapUnits ) ), 3 );
41  mMapUnits = mapUnits;
42 }
43 
45 {
46  QgsDebugMsgLevel( QStringLiteral( "Map units returned as %1" ).arg( QString::number( mMapUnits ) ), 4 );
47  return mMapUnits;
48 }
49 
50 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  {
59  // convert meters to inches
60  conversionFactor = 39.3700787;
61  delta = mapExtent.xMaximum() - mapExtent.xMinimum();
62  break;
64  conversionFactor = 12.0;
65  delta = mapExtent.xMaximum() - mapExtent.xMinimum();
66  break;
68  // convert nautical miles to inches
69  conversionFactor = 72913.4;
70  delta = mapExtent.xMaximum() - mapExtent.xMinimum();
71  break;
72 
73  default:
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  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 
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  double lat = ( mapExtent.yMaximum() + mapExtent.yMinimum() ) * 0.5;
118  static const double RADS = ( 4.0 * std::atan( 1.0 ) ) / 180.0;
119  double a = std::pow( std::cos( lat * RADS ), 2 );
120  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  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  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 }
A rectangle specified with double values.
Definition: qgsrectangle.h:41
QgsScaleCalculator(double dpi=0, QgsUnitTypes::DistanceUnit mapUnits=QgsUnitTypes::DistanceMeters)
Constructor.
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
double calculateGeographicDistance(const QgsRectangle &mapExtent) const
Calculate the distance between two points in geographic coordinates.
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:265
void setDpi(double dpi)
Sets the dpi (dots per inch) for the output resolution, to be used in scale calculations.
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
void setMapUnits(QgsUnitTypes::DistanceUnit mapUnits)
Set the map units.
Degrees, for planar geographic CRS distance measurements.
Definition: qgsunittypes.h:62
double dpi()
Returns the DPI (dots per inch) used in scale calculations.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:177
DistanceUnit
Units of distance.
Definition: qgsunittypes.h:54
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:162
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:167
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:172
QgsUnitTypes::DistanceUnit mapUnits() const
Returns current map units.
double calculate(const QgsRectangle &mapExtent, double canvasWidth) const
Calculate the scale denominator.