Quantum GIS API Documentation  1.8
src/core/qgsdistancearea.cpp
Go to the documentation of this file.
00001 /***************************************************************************
00002   qgsdistancearea.cpp - Distance and area calculations on the ellipsoid
00003  ---------------------------------------------------------------------------
00004   Date                 : September 2005
00005   Copyright            : (C) 2005 by Martin Dobias
00006   email                : won.der at centrum.sk
00007  ***************************************************************************
00008  *                                                                         *
00009  *   This program is free software; you can redistribute it and/or modify  *
00010  *   it under the terms of the GNU General Public License as published by  *
00011  *   the Free Software Foundation; either version 2 of the License, or     *
00012  *   (at your option) any later version.                                   *
00013  *                                                                         *
00014  ***************************************************************************/
00015 
00016 #include <cmath>
00017 #include <sqlite3.h>
00018 #include <QDir>
00019 #include <QString>
00020 #include <QLocale>
00021 #include <QObject>
00022 
00023 #include "qgis.h"
00024 #include "qgspoint.h"
00025 #include "qgscoordinatetransform.h"
00026 #include "qgscoordinatereferencesystem.h"
00027 #include "qgsgeometry.h"
00028 #include "qgsdistancearea.h"
00029 #include "qgsapplication.h"
00030 #include "qgslogger.h"
00031 #include "qgsmessagelog.h"
00032 
00033 // MSVC compiler doesn't have defined M_PI in math.h
00034 #ifndef M_PI
00035 #define M_PI          3.14159265358979323846
00036 #endif
00037 
00038 #define DEG2RAD(x)    ((x)*M_PI/180)
00039 
00040 
00041 QgsDistanceArea::QgsDistanceArea()
00042 {
00043   // init with default settings
00044   mProjectionsEnabled = false;
00045   mCoordTransform = new QgsCoordinateTransform;
00046   setSourceCrs( GEOCRS_ID ); // WGS 84
00047   setEllipsoid( "WGS84" );
00048 }
00049 
00050 
00051 QgsDistanceArea::~QgsDistanceArea()
00052 {
00053   delete mCoordTransform;
00054 }
00055 
00056 
00057 void QgsDistanceArea::setProjectionsEnabled( bool flag )
00058 {
00059   mProjectionsEnabled = flag;
00060 }
00061 
00062 void QgsDistanceArea::setSourceCrs( long srsid )
00063 {
00064   QgsCoordinateReferenceSystem srcCRS;
00065   srcCRS.createFromSrsId( srsid );
00066   mCoordTransform->setSourceCrs( srcCRS );
00067 }
00068 
00069 void QgsDistanceArea::setSourceEpsgCrsId( long epsgId )
00070 {
00071   QgsCoordinateReferenceSystem srcCRS;
00072   srcCRS.createFromOgcWmsCrs( QString( "EPSG:%1" ).arg( epsgId ) );
00073   mCoordTransform->setSourceCrs( srcCRS );
00074 }
00075 
00076 void QgsDistanceArea::setSourceAuthId( QString authId )
00077 {
00078   QgsCoordinateReferenceSystem srcCRS;
00079   srcCRS.createFromOgcWmsCrs( authId );
00080   mCoordTransform->setSourceCrs( srcCRS );
00081 }
00082 
00083 bool QgsDistanceArea::setEllipsoid( const QString& ellipsoid )
00084 {
00085   QString radius, parameter2;
00086   //
00087   // SQLITE3 stuff - get parameters for selected ellipsoid
00088   //
00089   sqlite3      *myDatabase;
00090   const char   *myTail;
00091   sqlite3_stmt *myPreparedStatement;
00092   int           myResult;
00093 
00094   // Shortcut if ellipsoid is none.
00095   if ( ellipsoid == "NONE" )
00096   {
00097     mEllipsoid = "NONE";
00098     return true;
00099   }
00100 
00101   //check the db is available
00102   myResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().data(), &myDatabase, SQLITE_OPEN_READONLY, NULL );
00103   if ( myResult )
00104   {
00105     QgsMessageLog::logMessage( QObject::tr( "Can't open database: %1" ).arg( sqlite3_errmsg( myDatabase ) ) );
00106     // XXX This will likely never happen since on open, sqlite creates the
00107     //     database if it does not exist.
00108     return false;
00109   }
00110   // Set up the query to retrieve the projection information needed to populate the ELLIPSOID list
00111   QString mySql = "select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid + "'";
00112   myResult = sqlite3_prepare( myDatabase, mySql.toUtf8(), mySql.toUtf8().length(), &myPreparedStatement, &myTail );
00113   // XXX Need to free memory from the error msg if one is set
00114   if ( myResult == SQLITE_OK )
00115   {
00116     if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
00117     {
00118       radius = QString(( char * )sqlite3_column_text( myPreparedStatement, 0 ) );
00119       parameter2 = QString(( char * )sqlite3_column_text( myPreparedStatement, 1 ) );
00120     }
00121   }
00122   // close the sqlite3 statement
00123   sqlite3_finalize( myPreparedStatement );
00124   sqlite3_close( myDatabase );
00125 
00126   // row for this ellipsoid wasn't found?
00127   if ( radius.isEmpty() || parameter2.isEmpty() )
00128   {
00129     QgsDebugMsg( QString( "setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
00130     return false;
00131   }
00132 
00133   // get major semiaxis
00134   if ( radius.left( 2 ) == "a=" )
00135     mSemiMajor = radius.mid( 2 ).toDouble();
00136   else
00137   {
00138     QgsDebugMsg( QString( "setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
00139     return false;
00140   }
00141 
00142   // get second parameter
00143   // one of values 'b' or 'f' is in field parameter2
00144   // second one must be computed using formula: invf = a/(a-b)
00145   if ( parameter2.left( 2 ) == "b=" )
00146   {
00147     mSemiMinor = parameter2.mid( 2 ).toDouble();
00148     mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
00149   }
00150   else if ( parameter2.left( 3 ) == "rf=" )
00151   {
00152     mInvFlattening = parameter2.mid( 3 ).toDouble();
00153     mSemiMinor = mSemiMajor - ( mSemiMajor / mInvFlattening );
00154   }
00155   else
00156   {
00157     QgsDebugMsg( QString( "setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
00158     return false;
00159   }
00160 
00161   QgsDebugMsg( QString( "setEllipsoid: a=%1, b=%2, 1/f=%3" ).arg( mSemiMajor ).arg( mSemiMinor ).arg( mInvFlattening ) );
00162 
00163 
00164   // get spatial ref system for ellipsoid
00165   QString proj4 = "+proj=longlat +ellps=" + ellipsoid + " +no_defs";
00166   QgsCoordinateReferenceSystem destCRS;
00167   destCRS.createFromProj4( proj4 );
00168 
00169   // set transformation from project CRS to ellipsoid coordinates
00170   mCoordTransform->setDestCRS( destCRS );
00171 
00172   // precalculate some values for area calculations
00173   computeAreaInit();
00174 
00175   mEllipsoid = ellipsoid;
00176   return true;
00177 }
00178 
00179 
00180 double QgsDistanceArea::measure( QgsGeometry* geometry )
00181 {
00182   if ( !geometry )
00183     return 0.0;
00184 
00185   unsigned char* wkb = geometry->asWkb();
00186   if ( !wkb )
00187     return 0.0;
00188 
00189   unsigned char* ptr;
00190   unsigned int wkbType;
00191   double res, resTotal = 0;
00192   int count, i;
00193 
00194   memcpy( &wkbType, ( wkb + 1 ), sizeof( wkbType ) );
00195 
00196   // measure distance or area based on what is the type of geometry
00197   bool hasZptr = false;
00198 
00199   switch ( wkbType )
00200   {
00201     case QGis::WKBLineString25D:
00202       hasZptr = true;
00203     case QGis::WKBLineString:
00204       measureLine( wkb, &res, hasZptr );
00205       QgsDebugMsg( "returning " + QString::number( res ) );
00206       return res;
00207 
00208     case QGis::WKBMultiLineString25D:
00209       hasZptr = true;
00210     case QGis::WKBMultiLineString:
00211       count = *(( int* )( wkb + 5 ) );
00212       ptr = wkb + 9;
00213       for ( i = 0; i < count; i++ )
00214       {
00215         ptr = measureLine( ptr, &res, hasZptr );
00216         resTotal += res;
00217       }
00218       QgsDebugMsg( "returning " + QString::number( resTotal ) );
00219       return resTotal;
00220 
00221     case QGis::WKBPolygon25D:
00222       hasZptr = true;
00223     case QGis::WKBPolygon:
00224       measurePolygon( wkb, &res, 0, hasZptr );
00225       QgsDebugMsg( "returning " + QString::number( res ) );
00226       return res;
00227 
00228     case QGis::WKBMultiPolygon25D:
00229       hasZptr = true;
00230     case QGis::WKBMultiPolygon:
00231       count = *(( int* )( wkb + 5 ) );
00232       ptr = wkb + 9;
00233       for ( i = 0; i < count; i++ )
00234       {
00235         ptr = measurePolygon( ptr, &res, 0, hasZptr );
00236         resTotal += res;
00237       }
00238       QgsDebugMsg( "returning " + QString::number( resTotal ) );
00239       return resTotal;
00240 
00241     default:
00242       QgsDebugMsg( QString( "measure: unexpected geometry type: %1" ).arg( wkbType ) );
00243       return 0;
00244   }
00245 }
00246 
00247 double QgsDistanceArea::measurePerimeter( QgsGeometry* geometry )
00248 {
00249   if ( !geometry )
00250     return 0.0;
00251 
00252   unsigned char* wkb = geometry->asWkb();
00253   if ( !wkb )
00254     return 0.0;
00255 
00256   unsigned char* ptr;
00257   unsigned int wkbType;
00258   double res, resTotal = 0;
00259   int count, i;
00260 
00261   memcpy( &wkbType, ( wkb + 1 ), sizeof( wkbType ) );
00262 
00263   // measure distance or area based on what is the type of geometry
00264   bool hasZptr = false;
00265 
00266   switch ( wkbType )
00267   {
00268     case QGis::WKBLineString25D:
00269     case QGis::WKBLineString:
00270     case QGis::WKBMultiLineString25D:
00271     case QGis::WKBMultiLineString:
00272       return 0.0;
00273 
00274     case QGis::WKBPolygon25D:
00275       hasZptr = true;
00276     case QGis::WKBPolygon:
00277       measurePolygon( wkb, 0, &res, hasZptr );
00278       QgsDebugMsg( "returning " + QString::number( res ) );
00279       return res;
00280 
00281     case QGis::WKBMultiPolygon25D:
00282       hasZptr = true;
00283     case QGis::WKBMultiPolygon:
00284       count = *(( int* )( wkb + 5 ) );
00285       ptr = wkb + 9;
00286       for ( i = 0; i < count; i++ )
00287       {
00288         ptr = measurePolygon( ptr, 0, &res, hasZptr );
00289         resTotal += res;
00290       }
00291       QgsDebugMsg( "returning " + QString::number( resTotal ) );
00292       return resTotal;
00293 
00294     default:
00295       QgsDebugMsg( QString( "measure: unexpected geometry type: %1" ).arg( wkbType ) );
00296       return 0;
00297   }
00298 }
00299 
00300 
00301 unsigned char* QgsDistanceArea::measureLine( unsigned char* feature, double* area, bool hasZptr )
00302 {
00303   unsigned char *ptr = feature + 5;
00304   unsigned int nPoints = *(( int* )ptr );
00305   ptr = feature + 9;
00306 
00307   QList<QgsPoint> points;
00308   double x, y;
00309 
00310   QgsDebugMsg( "This feature WKB has " + QString::number( nPoints ) + " points" );
00311   // Extract the points from the WKB format into the vector
00312   for ( unsigned int i = 0; i < nPoints; ++i )
00313   {
00314     x = *(( double * ) ptr );
00315     ptr += sizeof( double );
00316     y = *(( double * ) ptr );
00317     ptr += sizeof( double );
00318     if ( hasZptr )
00319     {
00320       // totally ignore Z value
00321       ptr += sizeof( double );
00322     }
00323 
00324     points.append( QgsPoint( x, y ) );
00325   }
00326 
00327   *area = measureLine( points );
00328   return ptr;
00329 }
00330 
00331 double QgsDistanceArea::measureLine( const QList<QgsPoint>& points )
00332 {
00333   if ( points.size() < 2 )
00334     return 0;
00335 
00336   double total = 0;
00337   QgsPoint p1, p2;
00338 
00339   try
00340   {
00341     if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
00342       p1 = mCoordTransform->transform( points[0] );
00343     else
00344       p1 = points[0];
00345 
00346     for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
00347     {
00348       if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
00349       {
00350         p2 = mCoordTransform->transform( *i );
00351         total += computeDistanceBearing( p1, p2 );
00352       }
00353       else
00354       {
00355         p2 = *i;
00356         total += measureLine( p1, p2 );
00357       }
00358 
00359       p1 = p2;
00360     }
00361 
00362     return total;
00363   }
00364   catch ( QgsCsException &cse )
00365   {
00366     Q_UNUSED( cse );
00367     QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
00368     return 0.0;
00369   }
00370 
00371 }
00372 
00373 double QgsDistanceArea::measureLine( const QgsPoint& p1, const QgsPoint& p2 )
00374 {
00375   try
00376   {
00377     QgsPoint pp1 = p1, pp2 = p2;
00378     if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
00379     {
00380       pp1 = mCoordTransform->transform( p1 );
00381       pp2 = mCoordTransform->transform( p2 );
00382       return computeDistanceBearing( pp1, pp2 );
00383     }
00384     else
00385     {
00386       return sqrt(( p2.x() - p1.x() )*( p2.x() - p1.x() ) + ( p2.y() - p1.y() )*( p2.y() - p1.y() ) );
00387     }
00388   }
00389   catch ( QgsCsException &cse )
00390   {
00391     Q_UNUSED( cse );
00392     QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
00393     return 0.0;
00394   }
00395 }
00396 
00397 
00398 unsigned char* QgsDistanceArea::measurePolygon( unsigned char* feature, double* area, double* perimeter, bool hasZptr )
00399 {
00400   // get number of rings in the polygon
00401   unsigned int numRings = *(( int* )( feature + 1 + sizeof( int ) ) );
00402 
00403   if ( numRings == 0 )
00404     return 0;
00405 
00406   // Set pointer to the first ring
00407   unsigned char* ptr = feature + 1 + 2 * sizeof( int );
00408 
00409   QList<QgsPoint> points;
00410   QgsPoint pnt;
00411   double x, y;
00412   if ( area )
00413     *area = 0;
00414   if ( perimeter )
00415     *perimeter = 0;
00416 
00417   try
00418   {
00419     for ( unsigned int idx = 0; idx < numRings; idx++ )
00420     {
00421       int nPoints = *(( int* )ptr );
00422       ptr += 4;
00423 
00424       // Extract the points from the WKB and store in a pair of
00425       // vectors.
00426       for ( int jdx = 0; jdx < nPoints; jdx++ )
00427       {
00428         x = *(( double * ) ptr );
00429         ptr += sizeof( double );
00430         y = *(( double * ) ptr );
00431         ptr += sizeof( double );
00432         if ( hasZptr )
00433         {
00434           // totally ignore Z value
00435           ptr += sizeof( double );
00436         }
00437 
00438         pnt = QgsPoint( x, y );
00439 
00440         if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
00441         {
00442           pnt = mCoordTransform->transform( pnt );
00443         }
00444         points.append( pnt );
00445       }
00446 
00447       if ( points.size() > 2 )
00448       {
00449         if ( area )
00450         {
00451           double areaTmp = computePolygonArea( points );
00452           if ( idx == 0 )
00453           {
00454             // exterior ring
00455             *area += areaTmp;
00456           }
00457           else
00458           {
00459             *area -= areaTmp; // interior rings
00460           }
00461         }
00462 
00463         if ( perimeter )
00464         {
00465           *perimeter += measureLine( points );
00466         }
00467       }
00468 
00469       points.clear();
00470 
00471       if ( !area )
00472       {
00473         break;
00474       }
00475     }
00476   }
00477   catch ( QgsCsException &cse )
00478   {
00479     Q_UNUSED( cse );
00480     QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area or perimeter." ) );
00481   }
00482 
00483   return ptr;
00484 }
00485 
00486 
00487 double QgsDistanceArea::measurePolygon( const QList<QgsPoint>& points )
00488 {
00489 
00490   try
00491   {
00492     if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
00493     {
00494       QList<QgsPoint> pts;
00495       for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
00496       {
00497         pts.append( mCoordTransform->transform( *i ) );
00498       }
00499       return computePolygonArea( pts );
00500     }
00501     else
00502     {
00503       return computePolygonArea( points );
00504     }
00505   }
00506   catch ( QgsCsException &cse )
00507   {
00508     Q_UNUSED( cse );
00509     QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
00510     return 0.0;
00511   }
00512 }
00513 
00514 
00515 double QgsDistanceArea::bearing( const QgsPoint& p1, const QgsPoint& p2 )
00516 {
00517   QgsPoint pp1 = p1, pp2 = p2;
00518   double bearing;
00519 
00520   if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
00521   {
00522     pp1 = mCoordTransform->transform( p1 );
00523     pp2 = mCoordTransform->transform( p2 );
00524     computeDistanceBearing( pp1, pp2, &bearing );
00525   }
00526   else //compute simple planar azimuth
00527   {
00528     double dx = p2.x() - p1.x();
00529     double dy = p2.y() - p1.y();
00530     bearing = atan2( dx, dy );
00531   }
00532 
00533   return bearing;
00534 }
00535 
00536 
00538 // distance calculation
00539 
00540 double QgsDistanceArea::computeDistanceBearing(
00541   const QgsPoint& p1, const QgsPoint& p2,
00542   double* course1, double* course2 )
00543 {
00544   if ( p1.x() == p2.x() && p1.y() == p2.y() )
00545     return 0;
00546 
00547   // ellipsoid
00548   double a = mSemiMajor;
00549   double b = mSemiMinor;
00550   double f = 1 / mInvFlattening;
00551 
00552   double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
00553   double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
00554 
00555   double L = p2_lon - p1_lon;
00556   double U1 = atan(( 1 - f ) * tan( p1_lat ) );
00557   double U2 = atan(( 1 - f ) * tan( p2_lat ) );
00558   double sinU1 = sin( U1 ), cosU1 = cos( U1 );
00559   double sinU2 = sin( U2 ), cosU2 = cos( U2 );
00560   double lambda = L;
00561   double lambdaP = 2 * M_PI;
00562 
00563   double sinLambda = 0;
00564   double cosLambda = 0;
00565   double sinSigma = 0;
00566   double cosSigma = 0;
00567   double sigma = 0;
00568   double alpha = 0;
00569   double cosSqAlpha = 0;
00570   double cos2SigmaM = 0;
00571   double C = 0;
00572   double tu1 = 0;
00573   double tu2 = 0;
00574 
00575   int iterLimit = 20;
00576   while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
00577   {
00578     sinLambda = sin( lambda );
00579     cosLambda = cos( lambda );
00580     tu1 = ( cosU2 * sinLambda );
00581     tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
00582     sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
00583     cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
00584     sigma = atan2( sinSigma, cosSigma );
00585     alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
00586     cosSqAlpha = cos( alpha ) * cos( alpha );
00587     cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
00588     C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
00589     lambdaP = lambda;
00590     lambda = L + ( 1 - C ) * f * sin( alpha ) *
00591              ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
00592   }
00593 
00594   if ( iterLimit == 0 )
00595     return -1;  // formula failed to converge
00596 
00597   double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
00598   double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
00599   double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
00600   double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
00601                                        B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
00602   double s = b * A * ( sigma - deltaSigma );
00603 
00604   if ( course1 )
00605   {
00606     *course1 = atan2( tu1, tu2 );
00607   }
00608   if ( course2 )
00609   {
00610     // PI is added to return azimuth from P2 to P1
00611     *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
00612   }
00613 
00614   return s;
00615 }
00616 
00617 
00618 
00620 // stuff for measuring areas - copied from GRASS
00621 // don't know how does it work, but it's working .)
00622 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
00623 
00624 double QgsDistanceArea::getQ( double x )
00625 {
00626   double sinx, sinx2;
00627 
00628   sinx = sin( x );
00629   sinx2 = sinx * sinx;
00630 
00631   return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
00632 }
00633 
00634 
00635 double QgsDistanceArea::getQbar( double x )
00636 {
00637   double cosx, cosx2;
00638 
00639   cosx = cos( x );
00640   cosx2 = cosx * cosx;
00641 
00642   return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
00643 }
00644 
00645 
00646 void QgsDistanceArea::computeAreaInit()
00647 {
00648   double a2 = ( mSemiMajor * mSemiMajor );
00649   double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
00650   double e4, e6;
00651 
00652   m_TwoPI = M_PI + M_PI;
00653 
00654   e4 = e2 * e2;
00655   e6 = e4 * e2;
00656 
00657   m_AE = a2 * ( 1 - e2 );
00658 
00659   m_QA = ( 2.0 / 3.0 ) * e2;
00660   m_QB = ( 3.0 / 5.0 ) * e4;
00661   m_QC = ( 4.0 / 7.0 ) * e6;
00662 
00663   m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4  - ( 4.0 / 7.0 ) * e6;
00664   m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4  + ( 4.0 / 7.0 ) * e6;
00665   m_QbarC =                     - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
00666   m_QbarD = ( 4.0 / 49.0 ) * e6;
00667 
00668   m_Qp = getQ( M_PI / 2 );
00669   m_E  = 4 * M_PI * m_Qp * m_AE;
00670   if ( m_E < 0.0 )
00671     m_E = -m_E;
00672 }
00673 
00674 
00675 double QgsDistanceArea::computePolygonArea( const QList<QgsPoint>& points )
00676 {
00677   double x1, y1, x2, y2, dx, dy;
00678   double Qbar1, Qbar2;
00679   double area;
00680 
00681   QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
00682   if (( ! mProjectionsEnabled ) || ( mEllipsoid == "NONE" ) )
00683   {
00684     return computePolygonFlatArea( points );
00685   }
00686   int n = points.size();
00687   x2 = DEG2RAD( points[n-1].x() );
00688   y2 = DEG2RAD( points[n-1].y() );
00689   Qbar2 = getQbar( y2 );
00690 
00691   area = 0.0;
00692 
00693   for ( int i = 0; i < n; i++ )
00694   {
00695     x1 = x2;
00696     y1 = y2;
00697     Qbar1 = Qbar2;
00698 
00699     x2 = DEG2RAD( points[i].x() );
00700     y2 = DEG2RAD( points[i].y() );
00701     Qbar2 = getQbar( y2 );
00702 
00703     if ( x1 > x2 )
00704       while ( x1 - x2 > M_PI )
00705         x2 += m_TwoPI;
00706     else if ( x2 > x1 )
00707       while ( x2 - x1 > M_PI )
00708         x1 += m_TwoPI;
00709 
00710     dx = x2 - x1;
00711     area += dx * ( m_Qp - getQ( y2 ) );
00712 
00713     if (( dy = y2 - y1 ) != 0.0 )
00714       area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
00715   }
00716   if (( area *= m_AE ) < 0.0 )
00717     area = -area;
00718 
00719   /* kludge - if polygon circles the south pole the area will be
00720   * computed as if it cirlced the north pole. The correction is
00721   * the difference between total surface area of the earth and
00722   * the "north pole" area.
00723   */
00724   if ( area > m_E )
00725     area = m_E;
00726   if ( area > m_E / 2 )
00727     area = m_E - area;
00728 
00729   return area;
00730 }
00731 
00732 double QgsDistanceArea::computePolygonFlatArea( const QList<QgsPoint>& points )
00733 {
00734   // Normal plane area calculations.
00735   double area = 0.0;
00736   int i, size;
00737 
00738   size = points.size();
00739 
00740   // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
00741   for ( i = 0; i < size; i++ )
00742   {
00743     // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
00744     // Using '% size', so that we always end with the starting point
00745     // and thus close the polygon.
00746     area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
00747   }
00748   // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
00749   area = area / 2.0;
00750   return qAbs( area ); // All areas are positive!
00751 }
00752 
00753 QString QgsDistanceArea::textUnit( double value, int decimals, QGis::UnitType u, bool isArea, bool keepBaseUnit )
00754 {
00755   QString unitLabel;
00756 
00757   switch ( u )
00758   {
00759     case QGis::Meters:
00760       if ( isArea )
00761       {
00762         if ( keepBaseUnit )
00763         {
00764           unitLabel = QObject::trUtf8( " m²" );
00765         }
00766         else if ( qAbs( value ) > 1000000.0 )
00767         {
00768           unitLabel = QObject::trUtf8( " km²" );
00769           value = value / 1000000.0;
00770         }
00771         else if ( qAbs( value ) > 10000.0 )
00772         {
00773           unitLabel = QObject::tr( " ha" );
00774           value = value / 10000.0;
00775         }
00776         else
00777         {
00778           unitLabel = QObject::trUtf8( " m²" );
00779         }
00780       }
00781       else
00782       {
00783         if ( keepBaseUnit || qAbs( value ) == 0.0 )
00784         {
00785           unitLabel = QObject::tr( " m" );
00786         }
00787         else if ( qAbs( value ) > 1000.0 )
00788         {
00789           unitLabel = QObject::tr( " km" );
00790           value = value / 1000;
00791         }
00792         else if ( qAbs( value ) < 0.01 )
00793         {
00794           unitLabel = QObject::tr( " mm" );
00795           value = value * 1000;
00796         }
00797         else if ( qAbs( value ) < 0.1 )
00798         {
00799           unitLabel = QObject::tr( " cm" );
00800           value = value * 100;
00801         }
00802         else
00803         {
00804           unitLabel = QObject::tr( " m" );
00805         }
00806       }
00807       break;
00808     case QGis::Feet:
00809       if ( isArea )
00810       {
00811         if ( keepBaseUnit  || qAbs( value ) <= 0.5*43560.0 )
00812         {
00813           // < 0.5 acre show sq ft
00814           unitLabel = QObject::tr( " sq ft" );
00815         }
00816         else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
00817         {
00818           // < 0.5 sq mile show acre
00819           unitLabel = QObject::tr( " acres" );
00820           value /= 43560.0;
00821         }
00822         else
00823         {
00824           // above 0.5 acre show sq mi
00825           unitLabel = QObject::tr( " sq mile" );
00826           value /= 5280.0 * 5280.0;
00827         }
00828       }
00829       else
00830       {
00831         if ( qAbs( value ) <= 528.0 || keepBaseUnit )
00832         {
00833           if ( qAbs( value ) == 1.0 )
00834           {
00835             unitLabel = QObject::tr( " foot" );
00836           }
00837           else
00838           {
00839             unitLabel = QObject::tr( " feet" );
00840           }
00841         }
00842         else
00843         {
00844           unitLabel = QObject::tr( " mile" );
00845           value /= 5280.0;
00846         }
00847       }
00848       break;
00849     case QGis::Degrees:
00850       if ( isArea )
00851       {
00852         unitLabel = QObject::tr( " sq.deg." );
00853       }
00854       else
00855       {
00856         if ( qAbs( value ) == 1.0 )
00857           unitLabel = QObject::tr( " degree" );
00858         else
00859           unitLabel = QObject::tr( " degrees" );
00860       }
00861       break;
00862     case QGis::UnknownUnit:
00863       unitLabel = QObject::tr( " unknown" );
00864     default:
00865       QgsDebugMsg( QString( "Error: not picked up map units - actual value = %1" ).arg( u ) );
00866   };
00867 
00868 
00869   return QLocale::system().toString( value, 'f', decimals ) + unitLabel;
00870 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines