QGIS API Documentation  2.17.0-Master (973e4b0)
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  QgsDebugMsg( QString( "Map units set to %1" ).arg( QString::number( mapUnits ) ) );
41  mMapUnits = mapUnits;
42 }
43 
45 {
46  QgsDebugMsgLevel( QString( "Map units returned as %1" ).arg( QString::number( mMapUnits ) ), 4 );
47  return mMapUnits;
48 }
49 
50 double QgsScaleCalculator::calculate( const QgsRectangle &mapExtent, int canvasWidth )
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 QGis::Meters:
59  // convert meters to inches
60  conversionFactor = 39.3700787;
61  delta = mapExtent.xMaximum() - mapExtent.xMinimum();
62  break;
63  case QGis::Feet:
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:
74  case QGis::Degrees:
75  // degrees require conversion to meters first
76  conversionFactor = 39.3700787;
77  delta = calculateGeographicDistance( mapExtent );
78  break;
79  }
80  if ( canvasWidth == 0 || qgsDoubleNear( mDpi, 0.0 ) )
81  {
82  QgsDebugMsg( "Can't calculate scale from the input values" );
83  return 0;
84  }
85  double scale = ( delta * conversionFactor ) / ( static_cast< double >( canvasWidth ) / mDpi );
86  QgsDebugMsg( QString( "scale = %1 conversionFactor = %2" ).arg( scale ).arg( conversionFactor ) );
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  const static double rads = ( 4.0 * atan( 1.0 ) ) / 180.0;
119  double a = pow( cos( lat * rads ), 2 );
120  double c = 2.0 * atan2( sqrt( a ), sqrt( 1.0 - a ) );
121  const static 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  const static double e = 0.0810820288;
125  double radius = ra * ( 1.0 - e * e ) /
126  pow( 1.0 - e * e * sin( lat * rads ) * sin( lat * rads ), 1.5 );
127  double meters = ( mapExtent.xMaximum() - mapExtent.xMinimum() ) / 180.0 * radius * c;
128 
129  QgsDebugMsg( "Distance across map extent (m): " + QString::number( meters ) );
130 
131  return meters;
132 }
void setMapUnits(QGis::UnitType mapUnits)
Set the map units.
double calculateGeographicDistance(const QgsRectangle &mapExtent)
Calculate the distance between two points in geographic coordinates.
A rectangle specified with double values.
Definition: qgsrectangle.h:35
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:197
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
void setDpi(double dpi)
Set the dpi to be used in scale calculations.
QGis::UnitType mapUnits() const
Returns current map units.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:353
QString number(int n, int base)
double calculate(const QgsRectangle &mapExtent, int canvasWidth)
Calculate the scale denominator.
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:202
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:187
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
double dpi()
Accessor for dpi used in scale calculations.
QgsScaleCalculator(double dpi=0, QGis::UnitType mapUnits=QGis::Meters)
Constructor.
UnitType
Map units that qgis supports.
Definition: qgis.h:159
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:192