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