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