QGIS API Documentation  2.8.2-Wien
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 *
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 ), mMapUnits( mapUnits )
26 {}
27
29 {}
30
31 void QgsScaleCalculator::setDpi( double dpi )
32 {
33  mDpi = dpi;
34 }
36 {
37  return mDpi;
38 }
39
41 {
42  QgsDebugMsg( QString( "Map units set to %1" ).arg( QString::number( mapUnits ) ) );
43  mMapUnits = mapUnits;
44 }
45
47 {
48  QgsDebugMsgLevel( QString( "Map units returned as %1" ).arg( QString::number( mMapUnits ) ), 4 );
49  return mMapUnits;
50 }
51
52 double QgsScaleCalculator::calculate( const QgsRectangle &mapExtent, int canvasWidth )
53 {
54  double conversionFactor = 0;
55  double delta = 0;
56  // calculation is based on the map units and extent, the dpi of the
57  // users display, and the canvas width
58  switch ( mMapUnits )
59  {
60  case QGis::Meters:
61  // convert meters to inches
62  conversionFactor = 39.3700787;
63  delta = mapExtent.xMaximum() - mapExtent.xMinimum();
64  break;
65  case QGis::Feet:
66  conversionFactor = 12.0;
67  delta = mapExtent.xMaximum() - mapExtent.xMinimum();
68  break;
70  // convert nautical miles to inches
71  conversionFactor = 72913.4;
72  delta = mapExtent.xMaximum() - mapExtent.xMinimum();
73  break;
74
75  default:
76  case QGis::Degrees:
77  // degrees require conversion to meters first
78  conversionFactor = 39.3700787;
79  delta = calculateGeographicDistance( mapExtent );
80  break;
81  }
82  if ( canvasWidth == 0 || mDpi == 0 )
83  {
84  QgsDebugMsg( "Can't calculate scale from the input values" );
85  return 0;
86  }
87  double scale = ( delta * conversionFactor ) / (( double )canvasWidth / mDpi );
88  QgsDebugMsg( QString( "scale = %1 conversionFactor = %2" ).arg( scale ).arg( conversionFactor ) );
89  return scale;
90 }
91
92
94 {
95  // need to calculate the x distance in meters
96  // We'll use the middle latitude for the calculation
97  // Note this is an approximation (although very close) but calculating scale
98  // for geographic data over large extents is quasi-meaningless
99
100  // The distance between two points on a sphere can be estimated
101  // using the Haversine formula. This gives the shortest distance
102  // between two points on the sphere. However, what we're after is
103  // the distance from the left of the given extent and the right of
104  // it. This is not necessarily the shortest distance between two
105  // points on a sphere.
106  //
107  // The code below uses the Haversine formula, but with some changes
108  // to cope with the above problem, and also to deal with the extent
109  // possibly extending beyond +/-180 degrees:
110  //
111  // - Use the Halversine formula to calculate the distance from -90 to
112  // +90 degrees at the mean latitude.
113  // - Scale this distance by the number of degrees between
114  // mapExtent.xMinimum() and mapExtent.xMaximum();
115  // - For a slight improvemnt, allow for the ellipsoid shape of earth.
116
117
118  // For a longitude change of 180 degrees
119  double lat = ( mapExtent.yMaximum() + mapExtent.yMinimum() ) * 0.5;
120  const static double rads = ( 4.0 * atan( 1.0 ) ) / 180.0;
121  double a = pow( cos( lat * rads ), 2 );
122  double c = 2.0 * atan2( sqrt( a ), sqrt( 1.0 - a ) );
123  const static double ra = 6378000; // [m]
124  // The eccentricity. This comes from sqrt(1.0 - rb*rb/(ra*ra)) with rb set
125  // to 6357000 m.
126  const static double e = 0.0810820288;
127  double radius = ra * ( 1.0 - e * e ) /
128  pow( 1.0 - e * e * sin( lat * rads ) * sin( lat * rads ), 1.5 );
129  double meters = ( mapExtent.xMaximum() - mapExtent.xMinimum() ) / 180.0 * radius * c;
130
131  QgsDebugMsg( "Distance across map extent (m): " + QString::number( meters ) );
132
133  return meters;
134 }