Quantum GIS API Documentation
1.8
|
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 }