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