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