QGIS API Documentation  2.0.1-Dufour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
qgsdistancearea.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsdistancearea.cpp - Distance and area calculations on the ellipsoid
3  ---------------------------------------------------------------------------
4  Date : September 2005
5  Copyright : (C) 2005 by Martin Dobias
6  email : won.der at centrum.sk
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include <cmath>
17 #include <sqlite3.h>
18 #include <QDir>
19 #include <QString>
20 #include <QLocale>
21 #include <QObject>
22 
23 #include "qgis.h"
24 #include "qgspoint.h"
25 #include "qgscoordinatetransform.h"
27 #include "qgsgeometry.h"
28 #include "qgsdistancearea.h"
29 #include "qgsapplication.h"
30 #include "qgslogger.h"
31 #include "qgsmessagelog.h"
32 
33 // MSVC compiler doesn't have defined M_PI in math.h
34 #ifndef M_PI
35 #define M_PI 3.14159265358979323846
36 #endif
37 
38 #define DEG2RAD(x) ((x)*M_PI/180)
39 
40 
42 {
43  // init with default settings
44  mEllipsoidalMode = false;
46  setSourceCrs( GEOCRS_ID ); // WGS 84
48 }
49 
50 
53 {
54  _copy( origDA );
55 }
56 
58 {
59  delete mCoordTransform;
60 }
61 
64 {
65  if ( this == & origDA )
66  {
67  // Do not copy unto self
68  return *this;
69  }
70  _copy( origDA );
71  return *this;
72 }
73 
76 {
78  mEllipsoid = origDA.mEllipsoid;
79  mSemiMajor = origDA.mSemiMajor;
80  mSemiMinor = origDA.mSemiMinor;
82  // Some calculations and trig. Should not be TOO time consuming.
83  // Alternatively we could copy the temp vars?
87 }
88 
90 {
91  mEllipsoidalMode = flag;
92 }
93 
95 {
97  srcCRS.createFromSrsId( srsid );
98  mCoordTransform->setSourceCrs( srcCRS );
99 }
100 
101 void QgsDistanceArea::setSourceAuthId( QString authId )
102 {
104  srcCRS.createFromOgcWmsCrs( authId );
105  mCoordTransform->setSourceCrs( srcCRS );
106 }
107 
108 bool QgsDistanceArea::setEllipsoid( const QString& ellipsoid )
109 {
110  QString radius, parameter2;
111  //
112  // SQLITE3 stuff - get parameters for selected ellipsoid
113  //
114  sqlite3 *myDatabase;
115  const char *myTail;
116  sqlite3_stmt *myPreparedStatement;
117  int myResult;
118 
119  // Shortcut if ellipsoid is none.
120  if ( ellipsoid == GEO_NONE )
121  {
123  return true;
124  }
125 
126  // Check if we have a custom projection, and set from text string.
127  // Format is "PARAMETER:<semi-major axis>:<semi minor axis>
128  // Numbers must be with (optional) decimal point and no other separators (C locale)
129  // Distances in meters. Flattening is calculated.
130  if ( ellipsoid.startsWith( "PARAMETER" ) )
131  {
132  QStringList paramList = ellipsoid.split( ":" );
133  bool semiMajorOk, semiMinorOk;
134  double semiMajor = paramList[1].toDouble( & semiMajorOk );
135  double semiMinor = paramList[2].toDouble( & semiMinorOk );
136  if ( semiMajorOk && semiMinorOk )
137  {
138  return setEllipsoid( semiMajor, semiMinor );
139  }
140  else
141  {
142  return false;
143  }
144  }
145 
146  // Continue with PROJ.4 list of ellipsoids.
147 
148  //check the db is available
149  myResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().data(), &myDatabase, SQLITE_OPEN_READONLY, NULL );
150  if ( myResult )
151  {
152  QgsMessageLog::logMessage( QObject::tr( "Can't open database: %1" ).arg( sqlite3_errmsg( myDatabase ) ) );
153  // XXX This will likely never happen since on open, sqlite creates the
154  // database if it does not exist.
155  return false;
156  }
157  // Set up the query to retrieve the projection information needed to populate the ELLIPSOID list
158  QString mySql = "select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid + "'";
159  myResult = sqlite3_prepare( myDatabase, mySql.toUtf8(), mySql.toUtf8().length(), &myPreparedStatement, &myTail );
160  // XXX Need to free memory from the error msg if one is set
161  if ( myResult == SQLITE_OK )
162  {
163  if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
164  {
165  radius = QString(( char * )sqlite3_column_text( myPreparedStatement, 0 ) );
166  parameter2 = QString(( char * )sqlite3_column_text( myPreparedStatement, 1 ) );
167  }
168  }
169  // close the sqlite3 statement
170  sqlite3_finalize( myPreparedStatement );
171  sqlite3_close( myDatabase );
172 
173  // row for this ellipsoid wasn't found?
174  if ( radius.isEmpty() || parameter2.isEmpty() )
175  {
176  QgsDebugMsg( QString( "setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
177  return false;
178  }
179 
180  // get major semiaxis
181  if ( radius.left( 2 ) == "a=" )
182  mSemiMajor = radius.mid( 2 ).toDouble();
183  else
184  {
185  QgsDebugMsg( QString( "setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
186  return false;
187  }
188 
189  // get second parameter
190  // one of values 'b' or 'f' is in field parameter2
191  // second one must be computed using formula: invf = a/(a-b)
192  if ( parameter2.left( 2 ) == "b=" )
193  {
194  mSemiMinor = parameter2.mid( 2 ).toDouble();
196  }
197  else if ( parameter2.left( 3 ) == "rf=" )
198  {
199  mInvFlattening = parameter2.mid( 3 ).toDouble();
201  }
202  else
203  {
204  QgsDebugMsg( QString( "setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
205  return false;
206  }
207 
208  QgsDebugMsg( QString( "setEllipsoid: a=%1, b=%2, 1/f=%3" ).arg( mSemiMajor ).arg( mSemiMinor ).arg( mInvFlattening ) );
209 
210 
211  // get spatial ref system for ellipsoid
212  QString proj4 = "+proj=longlat +ellps=" + ellipsoid + " +no_defs";
214  destCRS.createFromProj4( proj4 );
215  //TODO: createFromProj4 used to save to the user database any new CRS
216  // this behavior was changed in order to separate creation and saving.
217  // Not sure if it necessary to save it here, should be checked by someone
218  // familiar with the code (should also give a more descriptive name to the generated CRS)
219  if ( destCRS.srsid() == 0 )
220  {
221  QString myName = QString( " * %1 (%2)" )
222  .arg( QObject::tr( "Generated CRS", "A CRS automatically generated from layer info get this prefix for description" ) )
223  .arg( destCRS.toProj4() );
224  destCRS.saveAsUserCRS( myName );
225  }
226  //
227 
228  // set transformation from project CRS to ellipsoid coordinates
229  mCoordTransform->setDestCRS( destCRS );
230 
231  // precalculate some values for area calculations
232  computeAreaInit();
233 
235  return true;
236 }
237 
239 // Inverse flattening is calculated with invf = a/(a-b)
240 // Also, b = a-(a/invf)
241 bool QgsDistanceArea::setEllipsoid( double semiMajor, double semiMinor )
242 {
243  mEllipsoid = QString( "PARAMETER:%1:%2" ).arg( semiMajor ).arg( semiMinor );
244  mSemiMajor = semiMajor;
245  mSemiMinor = semiMinor;
247 
248  computeAreaInit();
249 
250  return true;
251 }
252 
253 
254 
256 {
257  if ( !geometry )
258  return 0.0;
259 
260  const unsigned char* wkb = geometry->asWkb();
261  if ( !wkb )
262  return 0.0;
263 
264  const unsigned char* ptr;
265  unsigned int wkbType;
266  double res, resTotal = 0;
267  int count, i;
268 
269  memcpy( &wkbType, ( wkb + 1 ), sizeof( wkbType ) );
270 
271  // measure distance or area based on what is the type of geometry
272  bool hasZptr = false;
273 
274  switch ( wkbType )
275  {
277  hasZptr = true;
278  case QGis::WKBLineString:
279  measureLine( wkb, &res, hasZptr );
280  QgsDebugMsg( "returning " + QString::number( res ) );
281  return res;
282 
284  hasZptr = true;
286  count = *(( int* )( wkb + 5 ) );
287  ptr = wkb + 9;
288  for ( i = 0; i < count; i++ )
289  {
290  ptr = measureLine( ptr, &res, hasZptr );
291  resTotal += res;
292  }
293  QgsDebugMsg( "returning " + QString::number( resTotal ) );
294  return resTotal;
295 
296  case QGis::WKBPolygon25D:
297  hasZptr = true;
298  case QGis::WKBPolygon:
299  measurePolygon( wkb, &res, 0, hasZptr );
300  QgsDebugMsg( "returning " + QString::number( res ) );
301  return res;
302 
304  hasZptr = true;
306  count = *(( int* )( wkb + 5 ) );
307  ptr = wkb + 9;
308  for ( i = 0; i < count; i++ )
309  {
310  ptr = measurePolygon( ptr, &res, 0, hasZptr );
311  if ( !ptr )
312  {
313  QgsDebugMsg( "measurePolygon returned 0" );
314  break;
315  }
316  resTotal += res;
317  }
318  QgsDebugMsg( "returning " + QString::number( resTotal ) );
319  return resTotal;
320 
321  default:
322  QgsDebugMsg( QString( "measure: unexpected geometry type: %1" ).arg( wkbType ) );
323  return 0;
324  }
325 }
326 
328 {
329  if ( !geometry )
330  return 0.0;
331 
332  const unsigned char* wkb = geometry->asWkb();
333  if ( !wkb )
334  return 0.0;
335 
336  const unsigned char* ptr;
337  unsigned int wkbType;
338  double res = 0.0, resTotal = 0.0;
339  int count, i;
340 
341  memcpy( &wkbType, ( wkb + 1 ), sizeof( wkbType ) );
342 
343  // measure distance or area based on what is the type of geometry
344  bool hasZptr = false;
345 
346  switch ( wkbType )
347  {
349  case QGis::WKBLineString:
352  return 0.0;
353 
354  case QGis::WKBPolygon25D:
355  hasZptr = true;
356  case QGis::WKBPolygon:
357  measurePolygon( wkb, 0, &res, hasZptr );
358  QgsDebugMsg( "returning " + QString::number( res ) );
359  return res;
360 
362  hasZptr = true;
364  count = *(( int* )( wkb + 5 ) );
365  ptr = wkb + 9;
366  for ( i = 0; i < count; i++ )
367  {
368  ptr = measurePolygon( ptr, 0, &res, hasZptr );
369  if ( !ptr )
370  {
371  QgsDebugMsg( "measurePolygon returned 0" );
372  break;
373  }
374  resTotal += res;
375  }
376  QgsDebugMsg( "returning " + QString::number( resTotal ) );
377  return resTotal;
378 
379  default:
380  QgsDebugMsg( QString( "measure: unexpected geometry type: %1" ).arg( wkbType ) );
381  return 0;
382  }
383 }
384 
385 
386 const unsigned char* QgsDistanceArea::measureLine( const unsigned char* feature, double* area, bool hasZptr )
387 {
388  const unsigned char *ptr = feature + 5;
389  unsigned int nPoints = *(( int* )ptr );
390  ptr = feature + 9;
391 
392  QList<QgsPoint> points;
393  double x, y;
394 
395  QgsDebugMsg( "This feature WKB has " + QString::number( nPoints ) + " points" );
396  // Extract the points from the WKB format into the vector
397  for ( unsigned int i = 0; i < nPoints; ++i )
398  {
399  x = *(( double * ) ptr );
400  ptr += sizeof( double );
401  y = *(( double * ) ptr );
402  ptr += sizeof( double );
403  if ( hasZptr )
404  {
405  // totally ignore Z value
406  ptr += sizeof( double );
407  }
408 
409  points.append( QgsPoint( x, y ) );
410  }
411 
412  *area = measureLine( points );
413  return ptr;
414 }
415 
416 double QgsDistanceArea::measureLine( const QList<QgsPoint>& points )
417 {
418  if ( points.size() < 2 )
419  return 0;
420 
421  double total = 0;
422  QgsPoint p1, p2;
423 
424  try
425  {
426  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
427  p1 = mCoordTransform->transform( points[0] );
428  else
429  p1 = points[0];
430 
431  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
432  {
433  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
434  {
435  p2 = mCoordTransform->transform( *i );
436  total += computeDistanceBearing( p1, p2 );
437  }
438  else
439  {
440  p2 = *i;
441  total += measureLine( p1, p2 );
442  }
443 
444  p1 = p2;
445  }
446 
447  return total;
448  }
449  catch ( QgsCsException &cse )
450  {
451  Q_UNUSED( cse );
452  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
453  return 0.0;
454  }
455 
456 }
457 
458 double QgsDistanceArea::measureLine( const QgsPoint& p1, const QgsPoint& p2 )
459 {
460  double result;
461 
462  try
463  {
464  QgsPoint pp1 = p1, pp2 = p2;
465 
466  QgsDebugMsgLevel( QString( "Measuring from %1 to %2" ).arg( p1.toString( 4 ) ).arg( p2.toString( 4 ) ), 3 );
467  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
468  {
469  QgsDebugMsgLevel( QString( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
470  QgsDebugMsgLevel( QString( "From proj4 : %1" ).arg( mCoordTransform->sourceCrs().toProj4() ), 4 );
471  QgsDebugMsgLevel( QString( "To proj4 : %1" ).arg( mCoordTransform->destCRS().toProj4() ), 4 );
472  pp1 = mCoordTransform->transform( p1 );
473  pp2 = mCoordTransform->transform( p2 );
474  QgsDebugMsgLevel( QString( "New points are %1 and %2, calculating..." ).arg( pp1.toString( 4 ) ).arg( pp2.toString( 4 ) ), 4 );
475  result = computeDistanceBearing( pp1, pp2 );
476  }
477  else
478  {
479  QgsDebugMsgLevel( "Cartesian calculation on canvas coordinates", 4 );
480  result = sqrt(( p2.x() - p1.x() ) * ( p2.x() - p1.x() ) + ( p2.y() - p1.y() ) * ( p2.y() - p1.y() ) );
481  }
482  }
483  catch ( QgsCsException &cse )
484  {
485  Q_UNUSED( cse );
486  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
487  result = 0.0;
488  }
489  QgsDebugMsgLevel( QString( "The result was %1" ).arg( result ), 3 );
490  return result;
491 }
492 
493 
494 const unsigned char* QgsDistanceArea::measurePolygon( const unsigned char* feature, double* area, double* perimeter, bool hasZptr )
495 {
496  if ( !feature )
497  {
498  QgsDebugMsg( "no feature to measure" );
499  return 0;
500  }
501 
502  // get number of rings in the polygon
503  unsigned int numRings = *(( int* )( feature + 1 + sizeof( int ) ) );
504 
505  if ( numRings == 0 )
506  {
507  QgsDebugMsg( "no rings to measure" );
508  return 0;
509  }
510 
511  // Set pointer to the first ring
512  const unsigned char* ptr = feature + 1 + 2 * sizeof( int );
513 
514  QList<QgsPoint> points;
515  QgsPoint pnt;
516  double x, y;
517  if ( area )
518  *area = 0;
519  if ( perimeter )
520  *perimeter = 0;
521 
522  try
523  {
524  for ( unsigned int idx = 0; idx < numRings; idx++ )
525  {
526  int nPoints = *(( int* )ptr );
527  ptr += 4;
528 
529  // Extract the points from the WKB and store in a pair of
530  // vectors.
531  for ( int jdx = 0; jdx < nPoints; jdx++ )
532  {
533  x = *(( double * ) ptr );
534  ptr += sizeof( double );
535  y = *(( double * ) ptr );
536  ptr += sizeof( double );
537  if ( hasZptr )
538  {
539  // totally ignore Z value
540  ptr += sizeof( double );
541  }
542 
543  pnt = QgsPoint( x, y );
544 
545  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
546  {
547  pnt = mCoordTransform->transform( pnt );
548  }
549  points.append( pnt );
550  }
551 
552  if ( points.size() > 2 )
553  {
554  if ( area )
555  {
556  double areaTmp = computePolygonArea( points );
557  if ( idx == 0 )
558  {
559  // exterior ring
560  *area += areaTmp;
561  }
562  else
563  {
564  *area -= areaTmp; // interior rings
565  }
566  }
567 
568  if ( perimeter )
569  {
570  if ( idx == 0 )
571  {
572  // exterior ring
573  *perimeter += measureLine( points );
574  }
575  }
576  }
577 
578  points.clear();
579  }
580  }
581  catch ( QgsCsException &cse )
582  {
583  Q_UNUSED( cse );
584  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area or perimeter." ) );
585  }
586 
587  return ptr;
588 }
589 
590 
591 double QgsDistanceArea::measurePolygon( const QList<QgsPoint>& points )
592 {
593  try
594  {
595  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
596  {
597  QList<QgsPoint> pts;
598  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
599  {
600  pts.append( mCoordTransform->transform( *i ) );
601  }
602  return computePolygonArea( pts );
603  }
604  else
605  {
606  return computePolygonArea( points );
607  }
608  }
609  catch ( QgsCsException &cse )
610  {
611  Q_UNUSED( cse );
612  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
613  return 0.0;
614  }
615 }
616 
617 
618 double QgsDistanceArea::bearing( const QgsPoint& p1, const QgsPoint& p2 )
619 {
620  QgsPoint pp1 = p1, pp2 = p2;
621  double bearing;
622 
623  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
624  {
625  pp1 = mCoordTransform->transform( p1 );
626  pp2 = mCoordTransform->transform( p2 );
627  computeDistanceBearing( pp1, pp2, &bearing );
628  }
629  else //compute simple planar azimuth
630  {
631  double dx = p2.x() - p1.x();
632  double dy = p2.y() - p1.y();
633  bearing = atan2( dx, dy );
634  }
635 
636  return bearing;
637 }
638 
639 
641 // distance calculation
642 
644  const QgsPoint& p1, const QgsPoint& p2,
645  double* course1, double* course2 )
646 {
647  if ( p1.x() == p2.x() && p1.y() == p2.y() )
648  return 0;
649 
650  // ellipsoid
651  double a = mSemiMajor;
652  double b = mSemiMinor;
653  double f = 1 / mInvFlattening;
654 
655  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
656  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
657 
658  double L = p2_lon - p1_lon;
659  double U1 = atan(( 1 - f ) * tan( p1_lat ) );
660  double U2 = atan(( 1 - f ) * tan( p2_lat ) );
661  double sinU1 = sin( U1 ), cosU1 = cos( U1 );
662  double sinU2 = sin( U2 ), cosU2 = cos( U2 );
663  double lambda = L;
664  double lambdaP = 2 * M_PI;
665 
666  double sinLambda = 0;
667  double cosLambda = 0;
668  double sinSigma = 0;
669  double cosSigma = 0;
670  double sigma = 0;
671  double alpha = 0;
672  double cosSqAlpha = 0;
673  double cos2SigmaM = 0;
674  double C = 0;
675  double tu1 = 0;
676  double tu2 = 0;
677 
678  int iterLimit = 20;
679  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
680  {
681  sinLambda = sin( lambda );
682  cosLambda = cos( lambda );
683  tu1 = ( cosU2 * sinLambda );
684  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
685  sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
686  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
687  sigma = atan2( sinSigma, cosSigma );
688  alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
689  cosSqAlpha = cos( alpha ) * cos( alpha );
690  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
691  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
692  lambdaP = lambda;
693  lambda = L + ( 1 - C ) * f * sin( alpha ) *
694  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
695  }
696 
697  if ( iterLimit == 0 )
698  return -1; // formula failed to converge
699 
700  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
701  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
702  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
703  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
704  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
705  double s = b * A * ( sigma - deltaSigma );
706 
707  if ( course1 )
708  {
709  *course1 = atan2( tu1, tu2 );
710  }
711  if ( course2 )
712  {
713  // PI is added to return azimuth from P2 to P1
714  *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
715  }
716 
717  return s;
718 }
719 
720 
721 
723 // stuff for measuring areas - copied from GRASS
724 // don't know how does it work, but it's working .)
725 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
726 
727 double QgsDistanceArea::getQ( double x )
728 {
729  double sinx, sinx2;
730 
731  sinx = sin( x );
732  sinx2 = sinx * sinx;
733 
734  return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
735 }
736 
737 
738 double QgsDistanceArea::getQbar( double x )
739 {
740  double cosx, cosx2;
741 
742  cosx = cos( x );
743  cosx2 = cosx * cosx;
744 
745  return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
746 }
747 
748 
750 {
751  double a2 = ( mSemiMajor * mSemiMajor );
752  double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
753  double e4, e6;
754 
755  m_TwoPI = M_PI + M_PI;
756 
757  e4 = e2 * e2;
758  e6 = e4 * e2;
759 
760  m_AE = a2 * ( 1 - e2 );
761 
762  m_QA = ( 2.0 / 3.0 ) * e2;
763  m_QB = ( 3.0 / 5.0 ) * e4;
764  m_QC = ( 4.0 / 7.0 ) * e6;
765 
766  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
767  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
768  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
769  m_QbarD = ( 4.0 / 49.0 ) * e6;
770 
771  m_Qp = getQ( M_PI / 2 );
772  m_E = 4 * M_PI * m_Qp * m_AE;
773  if ( m_E < 0.0 )
774  m_E = -m_E;
775 }
776 
777 
778 double QgsDistanceArea::computePolygonArea( const QList<QgsPoint>& points )
779 {
780  double x1, y1, x2, y2, dx, dy;
781  double Qbar1, Qbar2;
782  double area;
783 
784  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
785  if (( ! mEllipsoidalMode ) || ( mEllipsoid == GEO_NONE ) )
786  {
787  return computePolygonFlatArea( points );
788  }
789  int n = points.size();
790  x2 = DEG2RAD( points[n-1].x() );
791  y2 = DEG2RAD( points[n-1].y() );
792  Qbar2 = getQbar( y2 );
793 
794  area = 0.0;
795 
796  for ( int i = 0; i < n; i++ )
797  {
798  x1 = x2;
799  y1 = y2;
800  Qbar1 = Qbar2;
801 
802  x2 = DEG2RAD( points[i].x() );
803  y2 = DEG2RAD( points[i].y() );
804  Qbar2 = getQbar( y2 );
805 
806  if ( x1 > x2 )
807  while ( x1 - x2 > M_PI )
808  x2 += m_TwoPI;
809  else if ( x2 > x1 )
810  while ( x2 - x1 > M_PI )
811  x1 += m_TwoPI;
812 
813  dx = x2 - x1;
814  area += dx * ( m_Qp - getQ( y2 ) );
815 
816  if (( dy = y2 - y1 ) != 0.0 )
817  area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
818  }
819  if (( area *= m_AE ) < 0.0 )
820  area = -area;
821 
822  /* kludge - if polygon circles the south pole the area will be
823  * computed as if it cirlced the north pole. The correction is
824  * the difference between total surface area of the earth and
825  * the "north pole" area.
826  */
827  if ( area > m_E )
828  area = m_E;
829  if ( area > m_E / 2 )
830  area = m_E - area;
831 
832  return area;
833 }
834 
835 double QgsDistanceArea::computePolygonFlatArea( const QList<QgsPoint>& points )
836 {
837  // Normal plane area calculations.
838  double area = 0.0;
839  int i, size;
840 
841  size = points.size();
842 
843  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
844  for ( i = 0; i < size; i++ )
845  {
846  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
847  // Using '% size', so that we always end with the starting point
848  // and thus close the polygon.
849  area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
850  }
851  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
852  area = area / 2.0;
853  return qAbs( area ); // All areas are positive!
854 }
855 
856 QString QgsDistanceArea::textUnit( double value, int decimals, QGis::UnitType u, bool isArea, bool keepBaseUnit )
857 {
858  QString unitLabel;
859 
860  switch ( u )
861  {
862  case QGis::Meters:
863  if ( isArea )
864  {
865  if ( keepBaseUnit )
866  {
867  unitLabel = QObject::trUtf8( " m²" );
868  }
869  else if ( qAbs( value ) > 1000000.0 )
870  {
871  unitLabel = QObject::trUtf8( " km²" );
872  value = value / 1000000.0;
873  }
874  else if ( qAbs( value ) > 10000.0 )
875  {
876  unitLabel = QObject::tr( " ha" );
877  value = value / 10000.0;
878  }
879  else
880  {
881  unitLabel = QObject::trUtf8( " m²" );
882  }
883  }
884  else
885  {
886  if ( keepBaseUnit || qAbs( value ) == 0.0 )
887  {
888  unitLabel = QObject::tr( " m" );
889  }
890  else if ( qAbs( value ) > 1000.0 )
891  {
892  unitLabel = QObject::tr( " km" );
893  value = value / 1000;
894  }
895  else if ( qAbs( value ) < 0.01 )
896  {
897  unitLabel = QObject::tr( " mm" );
898  value = value * 1000;
899  }
900  else if ( qAbs( value ) < 0.1 )
901  {
902  unitLabel = QObject::tr( " cm" );
903  value = value * 100;
904  }
905  else
906  {
907  unitLabel = QObject::tr( " m" );
908  }
909  }
910  break;
911  case QGis::Feet:
912  if ( isArea )
913  {
914  if ( keepBaseUnit || qAbs( value ) <= 0.5*43560.0 )
915  {
916  // < 0.5 acre show sq ft
917  unitLabel = QObject::tr( " sq ft" );
918  }
919  else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
920  {
921  // < 0.5 sq mile show acre
922  unitLabel = QObject::tr( " acres" );
923  value /= 43560.0;
924  }
925  else
926  {
927  // above 0.5 acre show sq mi
928  unitLabel = QObject::tr( " sq mile" );
929  value /= 5280.0 * 5280.0;
930  }
931  }
932  else
933  {
934  if ( qAbs( value ) <= 528.0 || keepBaseUnit )
935  {
936  if ( qAbs( value ) == 1.0 )
937  {
938  unitLabel = QObject::tr( " foot" );
939  }
940  else
941  {
942  unitLabel = QObject::tr( " feet" );
943  }
944  }
945  else
946  {
947  unitLabel = QObject::tr( " mile" );
948  value /= 5280.0;
949  }
950  }
951  break;
952  case QGis::Degrees:
953  if ( isArea )
954  {
955  unitLabel = QObject::tr( " sq.deg." );
956  }
957  else
958  {
959  if ( qAbs( value ) == 1.0 )
960  unitLabel = QObject::tr( " degree" );
961  else
962  unitLabel = QObject::tr( " degrees" );
963  }
964  break;
965  case QGis::UnknownUnit:
966  unitLabel = QObject::tr( " unknown" );
967  default:
968  QgsDebugMsg( QString( "Error: not picked up map units - actual value = %1" ).arg( u ) );
969  };
970 
971 
972  return QLocale::system().toString( value, 'f', decimals ) + unitLabel;
973 }
974 
975 void QgsDistanceArea::convertMeasurement( double &measure, QGis::UnitType &measureUnits, QGis::UnitType displayUnits, bool isArea )
976 {
977  // Helper for converting between meters and feet
978  // The parameters measure and measureUnits are in/out
979 
980  if (( measureUnits == QGis::Degrees || measureUnits == QGis::Feet ) &&
981  mEllipsoid != GEO_NONE &&
983  {
984  // Measuring on an ellipsoid returned meters. Force!
985  measureUnits = QGis::Meters;
986  QgsDebugMsg( "We're measuring on an ellipsoid or using projections, the system is returning meters" );
987  }
988 
989  // Only convert between meters and feet
990  if ( measureUnits == QGis::Meters && displayUnits == QGis::Feet )
991  {
992  QgsDebugMsg( QString( "Converting %1 meters" ).arg( QString::number( measure ) ) );
993  measure /= 0.3048;
994  if ( isArea )
995  {
996  measure /= 0.3048;
997  }
998  QgsDebugMsg( QString( "to %1 feet" ).arg( QString::number( measure ) ) );
999  measureUnits = QGis::Feet;
1000  }
1001  if ( measureUnits == QGis::Feet && displayUnits == QGis::Meters )
1002  {
1003  QgsDebugMsg( QString( "Converting %1 feet" ).arg( QString::number( measure ) ) );
1004  measure *= 0.3048;
1005  if ( isArea )
1006  {
1007  measure *= 0.3048;
1008  }
1009  QgsDebugMsg( QString( "to %1 meters" ).arg( QString::number( measure ) ) );
1010  measureUnits = QGis::Meters;
1011  }
1012 }