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