Quantum GIS API Documentation
1.7.4
|
00001 /*************************************************************************** 00002 qgsscalecalculator.h 00003 Calculates scale based on map extent and units 00004 ------------------- 00005 begin : May 18, 2004 00006 copyright : (C) 2004 by Gary E.Sherman 00007 email : sherman at mrcc.com 00008 ***************************************************************************/ 00009 00010 /*************************************************************************** 00011 * * 00012 * This program is free software; you can redistribute it and/or modify * 00013 * it under the terms of the GNU General Public License as published by * 00014 * the Free Software Foundation; either version 2 of the License, or * 00015 * (at your option) any later version. * 00016 * * 00017 ***************************************************************************/ 00018 /* $Id$ */ 00019 00020 #include <cmath> 00021 #include "qgslogger.h" 00022 #include "qgsrectangle.h" 00023 #include "qgsscalecalculator.h" 00024 00025 QgsScaleCalculator::QgsScaleCalculator( double dpi, QGis::UnitType mapUnits ) 00026 : mDpi( dpi ), mMapUnits( mapUnits ) 00027 {} 00028 00029 QgsScaleCalculator::~QgsScaleCalculator() 00030 {} 00031 00032 void QgsScaleCalculator::setDpi( double dpi ) 00033 { 00034 mDpi = dpi; 00035 } 00036 double QgsScaleCalculator::dpi() 00037 { 00038 return mDpi; 00039 } 00040 00041 void QgsScaleCalculator::setMapUnits( QGis::UnitType mapUnits ) 00042 { 00043 QgsDebugMsg( QString( "Map units set to %1" ).arg( QString::number( mapUnits ) ) ); 00044 mMapUnits = mapUnits; 00045 } 00046 00047 QGis::UnitType QgsScaleCalculator::mapUnits() const 00048 { 00049 QgsDebugMsgLevel( QString( "Map units returned as %1" ).arg( QString::number( mMapUnits ) ), 3 ); 00050 return mMapUnits; 00051 } 00052 00053 double QgsScaleCalculator::calculate( const QgsRectangle &mapExtent, int canvasWidth ) 00054 { 00055 double conversionFactor = 0; 00056 double delta = 0; 00057 // calculation is based on the map units and extent, the dpi of the 00058 // users display, and the canvas width 00059 switch ( mMapUnits ) 00060 { 00061 case QGis::Meters: 00062 // convert meters to inches 00063 conversionFactor = 39.3700787; 00064 delta = mapExtent.xMaximum() - mapExtent.xMinimum(); 00065 break; 00066 case QGis::Feet: 00067 conversionFactor = 12.0; 00068 delta = mapExtent.xMaximum() - mapExtent.xMinimum(); 00069 break; 00070 case QGis::DecimalDegrees: 00071 // degrees require conversion to meters first 00072 conversionFactor = 39.3700787; 00073 delta = calculateGeographicDistance( mapExtent ); 00074 break; 00075 case QGis::DegreesMinutesSeconds: 00076 // degrees require conversion to meters first 00077 conversionFactor = 39.3700787; 00078 delta = calculateGeographicDistance( mapExtent ); 00079 break; 00080 case QGis::DegreesDecimalMinutes: 00081 // degrees require conversion to meters first 00082 conversionFactor = 39.3700787; 00083 delta = calculateGeographicDistance( mapExtent ); 00084 break; 00085 default: 00086 Q_ASSERT( "bad map units" ); 00087 break; 00088 } 00089 QgsDebugMsg( "Using conversionFactor of " + QString::number( conversionFactor ) ); 00090 if ( canvasWidth == 0 || mDpi == 0 ) 00091 { 00092 QgsDebugMsg( "Can't calculate scale from the input values" ); 00093 return 0; 00094 } 00095 double scale = ( delta * conversionFactor ) / (( double )canvasWidth / mDpi ); 00096 return scale; 00097 } 00098 00099 00100 double QgsScaleCalculator::calculateGeographicDistance( const QgsRectangle &mapExtent ) 00101 { 00102 // need to calculate the x distance in meters 00103 // We'll use the middle latitude for the calculation 00104 // Note this is an approximation (although very close) but calculating scale 00105 // for geographic data over large extents is quasi-meaningless 00106 00107 // The distance between two points on a sphere can be estimated 00108 // using the Haversine formula. This gives the shortest distance 00109 // between two points on the sphere. However, what we're after is 00110 // the distance from the left of the given extent and the right of 00111 // it. This is not necessarily the shortest distance between two 00112 // points on a sphere. 00113 // 00114 // The code below uses the Haversine formula, but with some changes 00115 // to cope with the above problem, and also to deal with the extent 00116 // possibly extending beyond +/-180 degrees: 00117 // 00118 // - Use the Halversine formula to calculate the distance from -90 to 00119 // +90 degrees at the mean latitude. 00120 // - Scale this distance by the number of degrees between 00121 // mapExtent.xMinimum() and mapExtent.xMaximum(); 00122 // - For a slight improvemnt, allow for the ellipsoid shape of earth. 00123 00124 00125 // For a longitude change of 180 degrees 00126 double lat = ( mapExtent.yMaximum() + mapExtent.yMinimum() ) * 0.5; 00127 const static double rads = ( 4.0 * atan( 1.0 ) ) / 180.0; 00128 double a = pow( cos( lat * rads ), 2 ); 00129 double c = 2.0 * atan2( sqrt( a ), sqrt( 1.0 - a ) ); 00130 const static double ra = 6378000; // [m] 00131 // The eccentricity. This comes from sqrt(1.0 - rb*rb/(ra*ra)) with rb set 00132 // to 6357000 m. 00133 const static double e = 0.0810820288; 00134 double radius = ra * ( 1.0 - e * e ) / 00135 pow( 1.0 - e * e * sin( lat * rads ) * sin( lat * rads ), 1.5 ); 00136 double meters = ( mapExtent.xMaximum() - mapExtent.xMinimum() ) / 180.0 * radius * c; 00137 00138 QgsDebugMsg( "Distance across map extent (m): " + QString::number( meters ) ); 00139 00140 return meters; 00141 }